SiPixelLorentzAngle/test/calculateLorentzAngleFromClusterSizeFpix.C

0001 double fitFunction(double *x, double *par)
0002 {
0003     double a;
0004     if(x[0] < par[0]) a = par[1] + (x[0] - par[0]) * par[2]  + (x[0] - par[0]) *(x[0] - par[0]) * par[4] + (x[0] - par[0]) * (x[0] - par[0]) *(x[0] - par[0]) * par[6];
0005     else        a = par[1] + (x[0] - par[0]) * par[3] + (x[0] - par[0]) * (x[0] - par[0]) * par[5] + (x[0] - par[0]) * (x[0] - par[0]) * (x[0] - par[0]) * par[7];
0006     return a;
0007 }
0008
0009 double fitFunction2(double *x, double *par)
0010 {
0011     double a;
0012     if(x[0] < par[0]) a = par[1] + (x[0] - par[0]) * par[2];
0013     else        a = par[1] + (x[0] - par[0]) * par[3];
0014     return a;
0015 }
0016
0017 double fitFunction3(double *x, double *par)
0018 {
0019     return par[1] + TMath::Sqrt(par[2] + par[3] * (x[0] - par[0]) * (x[0] - par[0]) );
0020 }
0021
0022 double fitf(double *x, double *par)
0023 {
0024 double arg;
0025     if(x[0] < par[3]) {
0026     arg = par[1]*par[1]+par[2]*par[2]*(x[0]-par[3])*(x[0]-par[3]);
0027     } else {
0028     arg = par[1]*par[1]+par[4]*par[4]*(x[0]-par[3])*(x[0]-par[3]);
0029     }
0030     double fitval = par[0]+sqrt(arg);
0031     return fitval;
0032 }
0033 int calculateLorentzAngleFromClusterSizeFpix()
0034 {
0035 //  setTDRStyle();

0036     cout << "hallo" << endl;
0037
0038 //  TFile *f = new TFile("/nfs/data5/wilke/TrackerPointing_ALL_V9/lorentzangle.root");
0039     TFile *f = new TFile("/data1/Users/wilke/LorentzAngle/CMSSW_2_2_3/lorentzangle_0T.root");
0040
0041     f->cd();
0042 //

0043     TF1 *f1 = new TF1("f1",fitFunction,60, 140, 8);
0044     f1->SetParName(0,"p0");
0045     f1->SetParName(1,"p1");
0046     f1->SetParName(2,"p2");
0047     f1->SetParName(3,"p3");
0048     f1->SetParName(4,"p4");
0049     f1->SetParName(5,"p5");
0050     f1->SetParameters(114,1,-0.01,0.01);
0051
0052     TF1 *f2 = new TF1("f2",fitFunction2, -1.0, 1.0, 4);
0053     f2->SetParName(0,"p0");
0054     f2->SetParName(1,"p1");
0055     f2->SetParName(2,"p2");
0056     f2->SetParName(3,"p3");
0057     f2->SetParameters(-0.4,1,-1.0,1.0);
0058
0059     TF1 *func = new TF1("func", fitf, -1, 1,5);//3.8T

0060     //TF1 *func = new TF1("func", fitf, -1.5, 1.5,5);//0T

0061     func->SetParameters(1.0,0.1,1.6,-0.4,1.2);
0062     func->SetParNames("Offset","RMS Constant","SlopeL","cot(#alpha)_min","SlopeR");
0063
0064
0065     TF1 *func_beta = new TF1("func_beta", fitf, -1., 1, 5);
    func_beta->SetParameters(1.,0.1,1.6,-0.,1.2);
      func_beta->SetParNames("Offset","RMS Constant","SlopeL","cot(beta)_min","SlopeR");
0066
0067     int hist_drift_ = 200;
0068     int hist_depth_ = 50;
0069     double min_drift_ = -1000;
0070     double max_drift_ = 1000;
0071     double min_depth_ = -100;
0072     double max_depth_ = 400;
0073     double width_ = 0.0285;
0074     int anglebins = 90;
0075     int anglebinscotan = 60;
0076     //int anglebinscotan = 36; //for 0T

0077
0078     TH2F * h_sizex_alpha = new TH2F("h_sizex_alpha", "h_sizex_alpha", anglebins, 0, 180,10 , .5, 10.5);
0079     TH2F * h_sizex_alpha_cotan = new TH2F("h_sizex_cotanalpha", "h_sizex_cotanalpha", anglebinscotan, -3, 3,10 , .5, 10.5);
0080     TH2F * h_sizey_beta_cotan = new TH2F("h_sizey_cotanbeta", "h_sizey_cotanbeta", anglebinscotan, -3, 3,10 , .5, 10.5);
0081     TH2F * h_alpha_beta_cotan = new TH2F("h_cotanalpha_cotanbeta", "h_cotanalpha_cotanbeta", 100, -5, 5,100 , -5, 5);
0082     TH2F * h_alpha_beta = new TH2F("h_alpha_beta", "h_alpha_beta", anglebins, -180, 180,anglebins, -180, 180);
0083     TH1F * h_alpha_cotan = new TH1F("h_cotanalpha", "h_cotanalpha", anglebinscotan, -3, 3 );
0084     TH1F * h_beta_cotan = new TH1F("h_cotanbeta", "h_cotanbeta", anglebinscotan, -3, 3 );
0085     TH1F * h_alpha = new TH1F("h_alpha", "h_alpha", anglebins, -180, 180 );
0086     TH1F * h_beta = new TH1F("h_beta", "h_beta", anglebins, -180, 180 );
0087     TH1F * h_chi2 = new TH1F("h_chi2", "h_chi2", 200, 0, 10 );
0088
0089     int run_;
0090     int event_;
0091     int module_;
0092     int side_;
0093     int panel_;
0094     //int ladder_;

0095     //int layer_;

0096     //int isflipped_;

0097     float pt_;
0098     float eta_;
0099     float phi_;
0100     double chi2_;
0101     double ndof_;
0102     const int maxpix = 100;
0103     struct Pixinfo
0104     {
0105         int npix;
0106         float row[200];
0107         float col[200];
0108         float adc[200];
0109         float x[200];
0110         float y[200];
0111     } pixinfo_;
0112
0113     struct Hit{
0114         float x;
0115         float y;
0116         double alpha;
0117         double beta;
0118         double gamma;
0119     };
0120     Hit simhit_, trackhit_;
0121
0122     struct Clust {
0123         float x;
0124         float y;
0125         float charge;
0126         int size_x;
0127         int size_y;
0128         int maxPixelCol;
0129         int maxPixelRow;
0130         int minPixelCol;
0131         int minPixelRow;
0132     } clust_;
0133
0134     struct Rechit {
0135         float x;
0136         float y;
0137     } rechit_;
0138
0139     // fill the histrograms with the ntpl

0140     TTree * LATree = (TTree*)f->Get("SiPixelLorentzAngleTreeForward_");
0141     int nentries = LATree->GetEntries();
0142     LATree->SetBranchAddress("run", &run_);
0143     LATree->SetBranchAddress("event", &event_);
0144     LATree->SetBranchAddress("module", &module_);
0145     LATree->SetBranchAddress("side", &side_);
0146     LATree->SetBranchAddress("panel", &panel_);
0147     //LATree->SetBranchAddress("ladder", &ladder_);

0148     //LATree->SetBranchAddress("layer", &layer_);

0149     //LATree->SetBranchAddress("isflipped", &isflipped_);

0150     LATree->SetBranchAddress("pt", &pt_);
0151     LATree->SetBranchAddress("eta", &eta_);
0152     LATree->SetBranchAddress("phi", &phi_);
0153     LATree->SetBranchAddress("chi2", &chi2_);
0154     LATree->SetBranchAddress("ndof", &ndof_);
0155     LATree->SetBranchAddress("trackhit", &trackhit_);
0156     LATree->SetBranchAddress("simhit", &simhit_);
0157     LATree->SetBranchAddress("npix", &pixinfo_.npix);
0158     LATree->SetBranchAddress("rowpix", pixinfo_.row);
0159     LATree->SetBranchAddress("colpix", pixinfo_.col);
0160     LATree->SetBranchAddress("adc", pixinfo_.adc);
0161     LATree->SetBranchAddress("xpix", pixinfo_.x);
0162     LATree->SetBranchAddress("ypix", pixinfo_.y);
0163     LATree->SetBranchAddress("clust", &clust_);
0164     LATree->SetBranchAddress("rechit", &rechit_);
0165
0166     cout << "Running over " << nentries << " hits" << endl;
0167     ofstream fAngles( "cotanangles.txt", ios::trunc );
0168     int passcut = 0;
0169     for(int ientrie = 0 ; ientrie < nentries ; ientrie++){
0170         LATree->GetEntry(ientrie);
0171         bool large_pix = false;
0172
0173         // is it a large pixel (needs to be excluded)

0174         for (int j = 0; j <  pixinfo_.npix; j++){
0175             int colpos = static_cast<int>(pixinfo_.col[j]);
0176             if (pixinfo_.row[j] == 0 || pixinfo_.row[j] == 79 || pixinfo_.row[j] == 80 || pixinfo_.row[j] == 159 || colpos % 52 == 0 || colpos % 52 == 51 ){
0177                 large_pix = true;
0178             }
0179         }
0180
0181
0182         if(clust_.size_y < 2) continue;
0183         //double residual = TMath::Sqrt( (trackhit_.x - rechit_.x) * (trackhit_.x - rechit_.x) + (trackhit_.y - rechit_.y) * (trackhit_.y - rechit_.y) );
0184         h_chi2->Fill((chi2_/ndof_));
0185         //if( (chi2_/ndof_) < 10 && !large_pix){

0186         if( (chi2_/ndof_) < 2 && !large_pix ){
0187         //if( (chi2_/ndof_) < 10 && side_ == 1 ){ // for -Z side

0188         //if( (chi2_/ndof_) < 10 && side_ == 2 ){ // for +Z side

0189         passcut++;
0190             if(trackhit_.alpha > 0){
0191                 h_sizex_alpha->Fill(trackhit_.alpha*180. / TMath::Pi(),clust_.size_x);
0192                 h_sizex_alpha_cotan->Fill(TMath::Tan(TMath::Pi()/2. - trackhit_.alpha),clust_.size_x);
0193                 h_sizey_beta_cotan->Fill(TMath::Tan(TMath::Pi()/2. - trackhit_.beta),clust_.size_y);
0194                 h_alpha_cotan->Fill(TMath::Tan(TMath::Pi()/2. - trackhit_.alpha));
0195                 h_beta_cotan->Fill(TMath::Tan(TMath::Pi()/2. - trackhit_.beta));
0196                 fAngles << TMath::Tan(TMath::Pi()/2. - trackhit_.alpha) << "\t" << TMath::Tan(TMath::Pi()/2. - trackhit_.beta) << endl;
0197             }
0198             else{
0199                 h_sizex_alpha->Fill( (trackhit_.alpha + TMath::Pi())*180. / TMath::Pi(),clust_.size_x);
0200                 h_sizex_alpha_cotan->Fill(TMath::Tan(TMath::Pi()/2. - trackhit_.alpha),clust_.size_x);
0201                 h_sizey_beta_cotan->Fill(TMath::Tan(TMath::Pi()/2. - trackhit_.beta),clust_.size_y);
0202                 h_alpha_cotan->Fill(TMath::Tan(TMath::Pi()/2. - trackhit_.alpha));
0203                 h_beta_cotan->Fill(TMath::Tan(TMath::Pi()/2. - trackhit_.beta));
0204                 fAngles << TMath::Tan(TMath::Pi()/2. - trackhit_.alpha) << "\t" << TMath::Tan(TMath::Pi()/2. - trackhit_.beta) << endl;
0205             }
0206             h_alpha->Fill( trackhit_.alpha*180. / TMath::Pi());
0207             h_beta->Fill( trackhit_.beta*180. / TMath::Pi());
0208             h_alpha_beta->Fill( trackhit_.alpha*180. / TMath::Pi(), trackhit_.beta*180. / TMath::Pi());
0209             if(TMath::Tan(TMath::Pi()/2. - trackhit_.alpha) < -0.3 && TMath::Tan(TMath::Pi()/2. - trackhit_.alpha) > -0.5) h_alpha_beta_cotan->Fill(TMath::Tan(TMath::Pi()/2. - trackhit_.alpha), TMath::Tan(TMath::Pi()/2. - trackhit_.beta));
0210
0211         }
0212     }
0213
0214     cout << TMath::Pi()/2. << endl;
0215     cout << "Passing selection " << passcut << " hits" << endl;
0216     TH1F * h_mean = new TH1F("h_mean","h_mean", anglebins, 0, 180);
0217     TH1F * h_slice_ = new TH1F("h_slice","h_slice", 10, .5, 10.5);
0218     TH1F * h_mean_cotan = new TH1F("h_mean_cotan","h_mean_contan", anglebinscotan, -3, 3);
0219     TH1F * h_slice_cotan_ = new TH1F("h_slice_cotan","h_slice_cotan", 10, .5, 10.5);
0220     TH1F * h_mean_cotan_beta = new TH1F("h_mean_cotan_beta","h_mean_contan_beta", anglebinscotan, -3, 3);
0221     TH1F * h_slice_cotan_beta = new TH1F("h_slice_cotan_beta","h_slice_cotan_beta", 10, .5, 10.5);
0222     //loop over bins in depth (z-local-coordinate) (in order to fit slices)

0223     for( int i = 1; i <= anglebins; i++){
0224 //              findMean(i, (i_module + (i_layer - 1) * 8));

0225         h_slice_->Reset("ICE");
0226
0227         // determine sigma and sigma^2 of the adc counts and average adc counts

0228             //loop over bins in drift width

0229         for( int j = 1; j<= 10; j++){
0230             h_slice_->SetBinContent(j,h_sizex_alpha->GetBinContent(i,j));
0231             h_slice_cotan_beta->SetBinContent(j,h_sizey_beta_cotan->GetBinContent(i,j));
0232         } // end loop over bins in drift width

0233
0234         double mean = h_slice_->GetMean(1);
0235       double error = h_slice_->GetMeanError(1);
0236
0237         h_mean->SetBinContent(i, mean);
0238         h_mean->SetBinError(i, error);
0239
0240         double mean_beta = h_slice_cotan_beta->GetMean(1);
0241         double error_beta = h_slice_cotan_beta->GetMeanError(1);
0242
0243         h_mean_cotan_beta->SetBinContent(i, mean_beta);
0244         h_mean_cotan_beta->SetBinError(i, error_beta);
0245     }// end loop over bins in depth

0246
0247     for( int i = 1; i <= anglebinscotan; i++){
0248         h_slice_cotan_->Reset("ICE");
0249
0250         //loop over bins in drift width

0251         for( int j = 1; j<= 10; j++){
0252             h_slice_cotan_->SetBinContent(j,h_sizex_alpha_cotan->GetBinContent(i,j));
0253         } // end loop over bins in drift width

0254
0255         double mean = h_slice_cotan_->GetMean(1);
0256         double error = h_slice_cotan_->GetMeanError(1);
0257
0258         h_mean_cotan->SetBinContent(i, mean);
0259         h_mean_cotan->SetBinError(i, error);
0260     }// end loop over bins in depth

0261
0262
0263     gStyle->SetOptStat(0);
0264     gStyle->SetOptTitle(0);
0265     TCanvas * c1 = new TCanvas("c1", "c1", 1200, 600);
0266     c1->Divide(2,1);
0267     c1->cd(1);
0268     h_sizex_alpha->GetXaxis()->SetTitle("#alpha [^{o}]");
0269     h_sizex_alpha->GetYaxis()->SetTitle("cluster size [pixels]");
0270     h_sizex_alpha->Draw("colz");
0271     c1->cd(2);
0272     h_mean->GetXaxis()->SetTitle("#alpha [^{o}]");
0273     h_mean->GetYaxis()->SetTitle("average cluster size [pixels]");
0274     h_mean->Draw();
0275     //h_mean->Fit(f1,"ERQ");

0276
0277     TCanvas * c2 = new TCanvas("c2", "c2", 1200, 600);
0278     c2->Divide(2,1);
0279     c2->cd(1);
0280     h_sizex_alpha_cotan->GetXaxis()->SetTitle("cotan(#alpha)");
0281     h_sizex_alpha_cotan->GetYaxis()->SetTitle("cluster size [pixels]");
0282     h_sizex_alpha_cotan->Draw("colz");
0283     c2->cd(2);
0284     h_mean_cotan->GetXaxis()->SetTitle("cotan(#alpha)");
0285     h_mean_cotan->GetYaxis()->SetTitle("average cluster size [pixels]");
0286     h_mean_cotan->Draw();
0287     //h_mean_cotan->Fit(f2,"ERQ");

0288     h_mean_cotan->Fit(func,"R");
0289
0290
0291     TCanvas * c3 = new TCanvas("c3", "c3", 600, 600);
0292     c3->Divide(1,1);
0293     c3->cd(1);
0294     //h_mean_cotan->GetXaxis()->SetTitle("cotan(#alpha)");

0295     //h_mean_cotan->GetYaxis()->SetTitle("average cluster size [pixels]");

0296     h_mean_cotan->Draw();
0297     //h_mean_cotan->Fit(func,"R");
0298
0299         TCanvas * c8 = new TCanvas("c8", "c8", 600, 600);
0300     h_mean_cotan_beta->GetXaxis()->SetTitle("cotan(#beta)");
0301     h_mean_cotan_beta->GetYaxis()->SetTitle("longitudinal cluster size [pixels]");
0302     h_mean_cotan_beta->Draw();
0303     h_mean_cotan_beta->Fit(func_beta,"ERQ");
0304
0305     TCanvas * c4 = new TCanvas("c4", "c4", 600, 600);
0306     h_alpha_cotan->Draw();
0307
0308     TCanvas * c5 = new TCanvas("c5", "c5", 600, 600);
0309     h_beta_cotan->Draw();
0310
0311     TCanvas * c6 = new TCanvas("c6", "c6", 600, 600);
0312     h_alpha->GetXaxis()->SetTitle("#alpha [^{o}]");
0313     h_alpha->GetYaxis()->SetTitle("number of hits");
0314     h_alpha->Draw();
0315
0316
0317     TCanvas * c7 = new TCanvas("c7", "c7", 600, 600);
0318     h_beta->GetXaxis()->SetTitle("#beta [^{o}]");
0319     h_beta->GetYaxis()->SetTitle("number of hits");
0320     h_beta->Draw();
0321
0322     TCanvas * c9 = new TCanvas("c9", "c9", 600, 600);
0323     h_alpha_beta->GetXaxis()->SetTitle("#alpha [^{o}]");
0324     h_alpha_beta->GetYaxis()->SetTitle("#beta [^{o}]");
0325     h_alpha_beta->Draw("col");
0326
0327     TCanvas * c10 = new TCanvas("c10", "c10", 600, 600);
0328     h_chi2->GetXaxis()->SetTitle("#chi^{2}/ndof");
0329     h_chi2->GetYaxis()->SetTitle("number of hits");
0330     h_chi2->Draw();
0331
0332     return 0;
0333 }