File indexing completed on 2024-04-06 11:59:36
0001 {
0002
0003
0004 TFile *f = new TFile("lorentzangle.root");
0005 f->cd();
0006
0007 TF1 *f1 = new TF1("f1","[0] + [1]*x",50., 235.);
0008 f1->SetParName(0,"p0");
0009 f1->SetParName(1,"p1");
0010 f1->SetParameter(0,0);
0011 f1->SetParameter(1,0.4);
0012
0013 int hist_drift_ = 200;
0014 int hist_depth_ = 50;
0015 double min_drift_ = -1000;
0016 double max_drift_ = 1000;
0017 double min_depth_ = -100;
0018 double max_depth_ = 400;
0019 double width_ = 0.0285;
0020
0021
0022
0023 TH2F * h_drift_depth_adc = new TH2F("h_drift_depth_adc", "h_drift_depth_adc",hist_drift_ , min_drift_, max_drift_, hist_depth_, min_depth_, max_depth_);
0024 TH2F * h_drift_depth_adc2 = new TH2F("h_drift_depth_adc2","h_drift_depth_adc2",hist_drift_ , min_drift_, max_drift_, hist_depth_, min_depth_, max_depth_);
0025 TH2F * h_drift_depth_noadc = new TH2F("h_drift_depth_noadc","h_drift_depth_noadc;drift in #mum;depth in #mum",hist_drift_ , min_drift_, max_drift_, hist_depth_, min_depth_, max_depth_);
0026
0027 int run_;
0028 int event_;
0029 int module_;
0030 int ladder_;
0031 int layer_;
0032 int isflipped_;
0033 float pt_;
0034 float eta_;
0035 float phi_;
0036 double chi2_;
0037 double ndof_;
0038 int maxpix = 200;
0039 struct Pixinfo
0040 {
0041 int npix;
0042 float row[maxpix];
0043 float col[maxpix];
0044 float adc[maxpix];
0045 float x[maxpix];
0046 float y[maxpix];
0047 } pixinfo_;
0048
0049 struct Hit{
0050 float x;
0051 float y;
0052 double alpha;
0053 double beta;
0054 double gamma;
0055 };
0056 Hit simhit_, trackhit_;
0057
0058 struct Clust {
0059 float x;
0060 float y;
0061 float charge;
0062 int size_x;
0063 int size_y;
0064 int maxPixelCol;
0065 int maxPixelRow;
0066 int minPixelCol;
0067 int minPixelRow;
0068 } clust_;
0069
0070 struct Rechit {
0071 float x;
0072 float y;
0073 } rechit_;
0074
0075
0076 TTree * LATree = (TTree*)f->Get("SiPixelLorentzAngleTree_");
0077 int nentries = LATree->GetEntries();
0078 LATree->SetBranchAddress("run", &run_);
0079 LATree->SetBranchAddress("event", &event_);
0080 LATree->SetBranchAddress("module", &module_);
0081 LATree->SetBranchAddress("ladder", &ladder_);
0082 LATree->SetBranchAddress("layer", &layer_);
0083 LATree->SetBranchAddress("isflipped", &isflipped_);
0084 LATree->SetBranchAddress("pt", &pt_);
0085 LATree->SetBranchAddress("eta", &eta_);
0086 LATree->SetBranchAddress("phi", &phi_);
0087 LATree->SetBranchAddress("chi2", &chi2_);
0088 LATree->SetBranchAddress("ndof", &ndof_);
0089 LATree->SetBranchAddress("trackhit", &trackhit_);
0090 LATree->SetBranchAddress("simhit", &simhit_);
0091 LATree->SetBranchAddress("npix", &pixinfo_.npix);
0092 LATree->SetBranchAddress("rowpix", pixinfo_.row);
0093 LATree->SetBranchAddress("colpix", pixinfo_.col);
0094 LATree->SetBranchAddress("adc", pixinfo_.adc);
0095 LATree->SetBranchAddress("xpix", pixinfo_.x);
0096 LATree->SetBranchAddress("ypix", pixinfo_.y);
0097 LATree->SetBranchAddress("clust", &clust_);
0098 LATree->SetBranchAddress("rechit", &rechit_);
0099
0100 cout << "Running over " << nentries << " hits" << endl;
0101
0102 for(int ientrie = 0 ; ientrie < nentries ; ientrie++){
0103 LATree->GetEntry(ientrie);
0104 bool large_pix = false;
0105 for (int j = 0; j < pixinfo_.npix; j++){
0106 int colpos = static_cast<int>(pixinfo_.col[j]);
0107 if (pixinfo_.row[j] == 0 || pixinfo_.row[j] == 79 || pixinfo_.row[j] == 80 || pixinfo_.row[j] == 159 || colpos % 52 == 0 || colpos % 52 == 51 ){
0108 large_pix = true;
0109 }
0110 }
0111
0112 double residual = TMath::Sqrt( (trackhit_.x - rechit_.x) * (trackhit_.x - rechit_.x) + (trackhit_.y - rechit_.y) * (trackhit_.y - rechit_.y) );
0113 if( (clust_.size_y >= 4) && (chi2_/ndof_) < 2 && !large_pix && residual < 0.005 && clust_.charge < 120){
0114 for (int j = 0; j < pixinfo_.npix; j++){
0115 float dx = (pixinfo_.x[j] - (trackhit_.x - width_/2. / TMath::Tan(trackhit_.alpha))) * 10000.;
0116 float dy = (pixinfo_.y[j] - (trackhit_.y - width_/2. / TMath::Tan(trackhit_.beta))) * 10000.;
0117 float depth = dy * tan(trackhit_.beta);
0118 float drift = dx - dy * tan(trackhit_.gamma);
0119 h_drift_depth_adc->Fill(drift, depth, pixinfo_.adc[j]);
0120 h_drift_depth_adc2->Fill(drift, depth, pixinfo_.adc[j]*pixinfo_.adc[j]);
0121 h_drift_depth_noadc->Fill(drift, depth);
0122 }
0123 }
0124 }
0125
0126 TH1F * h_mean = new TH1F("h_mean","h_mean;depth in #mum;drift in #mum", hist_depth_, min_depth_, max_depth_);
0127 TH1F * h_drift_depth_adc_slice_ = new TH1F("h_slice","h_slice", hist_drift_, min_drift_, max_drift_);
0128
0129 for( int i = 1; i <= hist_depth_; i++){
0130
0131 double npix = 0;
0132
0133 h_drift_depth_adc_slice_->Reset("ICE");
0134
0135
0136
0137 for( int j = 1; j<= hist_drift_; j++){
0138 if(h_drift_depth_noadc->GetBinContent(j, i) >= 1){
0139 double adc_error2 = (h_drift_depth_adc2->GetBinContent(j,i) - h_drift_depth_adc->GetBinContent(j,i)*h_drift_depth_adc->GetBinContent(j, i) / h_drift_depth_noadc->GetBinContent(j, i)) / h_drift_depth_noadc->GetBinContent(j, i);
0140 h_drift_depth_adc_slice_->SetBinContent(j, h_drift_depth_adc->GetBinContent(j,i));
0141 h_drift_depth_adc_slice_->SetBinError(j, sqrt(adc_error2));
0142 npix += h_drift_depth_noadc->GetBinContent(j,i);
0143 }else{
0144 h_drift_depth_adc_slice_->SetBinContent(j, h_drift_depth_adc->GetBinContent(j,i));
0145 h_drift_depth_adc_slice_->SetBinError(j, 0);
0146 }
0147 }
0148
0149 double mean = h_drift_depth_adc_slice_->GetMean(1);
0150 double error = 0;
0151 if(npix != 0){
0152 error = h_drift_depth_adc_slice_->GetRMS(1) / sqrt(npix);
0153 }
0154
0155 h_mean->SetBinContent(i, mean);
0156 h_mean->SetBinError(i, error);
0157 }
0158
0159 TCanvas * c1 = new TCanvas("c1", "c1", 1200, 600);
0160 c1->Divide(2,1);
0161 c1->cd(1);
0162
0163 h_drift_depth_noadc->Draw("colz");
0164
0165
0166 c1->cd(2);
0167 h_mean->Draw();
0168
0169
0170 h_mean->Fit(f1,"ERQ");
0171 double p0 = f1->GetParameter(0);
0172 double e0 = f1->GetParError(0);
0173 double p1 = f1->GetParameter(1);
0174 double e1 = f1->GetParError(1);
0175 double chi2 = f1->GetChisquare();
0176 double prob = f1->GetProb();
0177
0178
0179
0180 cout << "offset" << "\t" << "error" << "\t" << "slope" << "\t" << "error" << "\t" "rel.err" << "\t" << "chi2" << "\t" << "prob" << endl;
0181 cout << p0 << "\t" << e0 << "\t" << p1 << "\t" << e1 << "\t" << e1 / p1 *100. << "\t" << chi2 << "\t" << prob << endl;
0182
0183 }