Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 09:35:13 

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 fitFunction1(double *x, double *par)
0010 {
0011     double a;
0012     if(x[0] < par[0]) a = par[1] + (x[0] - par[0]) * par[2]  + (x[0] - par[0]) *(x[0] - par[0]) * par[4];
0013     else        a = par[1] + (x[0] - par[0]) * par[3] + (x[0] - par[0]) * (x[0] - par[0]) * par[5];
0014     return a;
0015 }
0016 
0017 double fitFunction2(double *x, double *par)
0018 {
0019     double a;
0020     if(x[0] < par[0]) a = par[1] + (x[0] - par[0]) * par[2];
0021     else        a = par[1] + (x[0] - par[0]) * par[3];
0022     return a;
0023 }
0024 
0025 double fitFunction3(double *x, double *par)
0026 {
0027     return par[1] + TMath::Sqrt(par[2] + par[3] * (x[0] - par[0]) * (x[0] - par[0]) );
0028 }
0029 
0030 Double_t fitf(Double_t *x,Double_t *par)
{
 Double_t arg;
      if(x[0] < par[3]) {
        arg = par[1]*par[1] + par[2]*par[2]*(x[0]-par[3])*(x[0]-par[3]);
   } else {
       arg = par[1]*par[1] + par[4]*par[4]*(x[0]-par[3])*(x[0]-par[3]);
   }
      Double_t fitval = par[0]+sqrt(arg);
    return fitval;
}
0031  
0032 int calculateLorentzAngleFromClusterSize()
0033 {
0034     
0035     int nhits = 0;
0036 //  setTDRStyle();

0037     cout << "hallo" << endl;
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     f->cd();
0041 //  

0042     TF1 *f1 = new TF1("f1",fitFunction,80, 150, 8); 
0043     f1->SetParName(0,"p0");
0044     f1->SetParName(1,"p1"); 
0045     f1->SetParName(2,"p2");
0046     f1->SetParName(3,"p3");
0047     f1->SetParName(4,"p4");
0048     f1->SetParName(5,"p5");
0049     f1->SetParameters(110,1,-0.01,0.01);
0050     
0051     TF1 *f2 = new TF1("f2",fitFunction1, -1., 1.0, 6); 
0052     f2->SetParName(0,"p0");
0053     f2->SetParName(1,"p1"); 
0054     f2->SetParName(2,"p2");
0055     f2->SetParName(3,"p3");
0056     f2->SetParameters(-0.5, 1., -1.,1.0);
0057     
0058     TF1 *func = new TF1("func", fitf, -1., 1.0, 5);
    func->SetParameters(1.,0.1,1.6,-0.4,1.2);
      func->SetParNames("Offset","RMS Constant","SlopeL","cot(alpha)_min","SlopeR");
0059     
0060     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");
     
0061     int hist_drift_ = 200;
0062     int hist_depth_ = 50;
0063     double min_drift_ = -1000;
0064     double max_drift_ = 1000;
0065     double min_depth_ = -100;
0066     double max_depth_ = 400;
0067     double width_ = 0.0285;
0068     int anglebins = 90; 
0069     int anglebinscotan = 120;
0070 
0071     TH2F * h_sizex_alpha = new TH2F("h_sizex_alpha", "h_sizex_alpha", anglebins, 0, 180,10 , .5, 10.5);
0072     TH2F * h_sizex_alpha_cotan = new TH2F("h_sizex_cotanalpha", "h_sizex_cotanalpha", anglebinscotan, -3, 3,10 , .5, 10.5);
0073     TH2F * h_sizey_beta_cotan = new TH2F("h_sizey_cotanbeta", "h_sizey_cotanbeta", anglebinscotan, -3, 3,10 , .5, 10.5);
0074     TH2F * h_alpha_beta_cotan = new TH2F("h_cotanalpha_cotanbeta", "h_cotanalpha_cotanbeta", 100, -5, 5,100 , -5, 5);
0075     TH2F * h_alpha_beta = new TH2F("h_alpha_beta", "h_alpha_beta", anglebins, -180, 180,anglebins, -180, 180);
0076     TH1F * h_alpha_cotan = new TH1F("h_cotanalpha", "h_cotanalpha", anglebinscotan, -3, 3 );
0077     TH1F * h_beta_cotan = new TH1F("h_cotanbeta", "h_cotanbeta", anglebinscotan, -3, 3 );
0078     TH1F * h_alpha = new TH1F("h_alpha", "h_alpha", anglebins, -180, 180 );
0079     TH1F * h_beta = new TH1F("h_beta", "h_beta", anglebins, -180, 180 );
0080     TH1F * h_chi2 = new TH1F("h_chi2", "h_chi2", 200, 0, 10 );
0081     TH1F * h_charge = new TH1F("h_chi2", "h_chi2", 200, 0, 200 );
0082     
0083     int run_;
0084     int event_;
0085     int module_;
0086     int ladder_;
0087     int layer_;
0088     int isflipped_;
0089     float pt_;
0090     float eta_;
0091     float phi_;
0092     double chi2_;
0093     double ndof_;
0094     const int maxpix = 300;
0095     struct Pixinfo
0096     {
0097         int npix;
0098         float row[300];
0099         float col[300];
0100         float adc[300];
0101         float x[300];
0102         float y[300];
0103     } pixinfo_;
0104     
0105     struct Hit{
0106         float x;
0107         float y;
0108         double alpha;
0109         double beta;
0110         double gamma;
0111     }; 
0112     Hit simhit_, trackhit_;
0113     
0114     struct Clust {
0115         float x;
0116         float y;
0117         float charge;
0118         int size_x;
0119         int size_y;
0120         int maxPixelCol;
0121         int maxPixelRow;
0122         int minPixelCol;
0123         int minPixelRow;
0124     } clust_;
0125     
0126     struct Rechit {
0127         float x;
0128         float y;
0129     } rechit_;
0130     
0131     // fill the histrograms with the ntpl

0132     TTree * LATree = (TTree*)f->Get("SiPixelLorentzAngleTree_");
0133     int nentries = LATree->GetEntries();
0134     LATree->SetBranchAddress("run", &run_);
0135     LATree->SetBranchAddress("event", &event_);
0136     LATree->SetBranchAddress("module", &module_);
0137     LATree->SetBranchAddress("ladder", &ladder_);
0138     LATree->SetBranchAddress("layer", &layer_);
0139     LATree->SetBranchAddress("isflipped", &isflipped_);
0140     LATree->SetBranchAddress("pt", &pt_);
0141     LATree->SetBranchAddress("eta", &eta_);
0142     LATree->SetBranchAddress("phi", &phi_);
0143     LATree->SetBranchAddress("chi2", &chi2_);
0144     LATree->SetBranchAddress("ndof", &ndof_);
0145     LATree->SetBranchAddress("trackhit", &trackhit_);
0146     LATree->SetBranchAddress("simhit", &simhit_);
0147     LATree->SetBranchAddress("npix", &pixinfo_.npix);
0148     LATree->SetBranchAddress("rowpix", pixinfo_.row);
0149     LATree->SetBranchAddress("colpix", pixinfo_.col);
0150     LATree->SetBranchAddress("adc", pixinfo_.adc);
0151     LATree->SetBranchAddress("xpix", pixinfo_.x);
0152     LATree->SetBranchAddress("ypix", pixinfo_.y);
0153     LATree->SetBranchAddress("clust", &clust_);
0154     LATree->SetBranchAddress("rechit", &rechit_);
0155     
0156     cout << "Running over " << nentries << " hits" << endl;
0157     ofstream fAngles( "cotanangles.txt", ios::trunc );
0158     for(int ientrie = 0 ; ientrie < nentries ; ientrie++){
0159         LATree->GetEntry(ientrie);  
0160         bool large_pix = false;
0161         
0162         // is it a large pixel (needs to be excluded)

0163         for (int j = 0; j <  pixinfo_.npix; j++){
0164             int colpos = static_cast<int>(pixinfo_.col[j]);
0165             if (pixinfo_.row[j] == 0 || pixinfo_.row[j] == 79 || pixinfo_.row[j] == 80 || pixinfo_.row[j] == 159 || colpos % 52 == 0 || colpos % 52 == 51 ){
0166                 large_pix = true;   
0167             }
0168         }
0169         
0170         // is it one of the problematic half ladders? (needs to be excluded)

0171         //if( (layer_ == 1 && (ladder_ == 5 || ladder_ == 6 || ladder_ == 15 || ladder_ == 16)) ||(layer_ == 2 && (ladder_ == 8 || ladder_ == 9 || ladder_ == 24 || ladder_ == 25)) ||(layer_ == 3 && (ladder_ ==11 || ladder_ == 12 || ladder_ == 33 || ladder_ == 34)) ) {

0172         //  continue;

0173     //s }

0174         
0175         if(clust_.size_y < 2) continue;
0176         
0177 //      double residual = TMath::Sqrt( (trackhit_.x - rechit_.x) * (trackhit_.x - rechit_.x) + (trackhit_.y - rechit_.y) * (trackhit_.y - rechit_.y) );
0178         //if(!isflipped_) continue;
0179         if(!large_pix){
0180             h_chi2->Fill((chi2_/ndof_));
0181             h_charge->Fill(clust_.charge);
0182         }
0183     
0184         if( (chi2_/ndof_) < 2 && !large_pix ){
0185             nhits++;
0186             if(trackhit_.alpha > 0){    
0187                 h_sizex_alpha->Fill(trackhit_.alpha*180. / TMath::Pi(),clust_.size_x);
0188                 h_sizex_alpha_cotan->Fill(TMath::Tan(TMath::Pi()/2. - trackhit_.alpha),clust_.size_x);
0189                 h_sizey_beta_cotan->Fill(TMath::Tan(TMath::Pi()/2. - trackhit_.beta),clust_.size_y);
0190                 h_alpha_cotan->Fill(TMath::Tan(TMath::Pi()/2. - trackhit_.alpha));
0191                 h_beta_cotan->Fill(TMath::Tan(TMath::Pi()/2. - trackhit_.beta));
0192                 fAngles << TMath::Tan(TMath::Pi()/2. - trackhit_.alpha) << "\t" << TMath::Tan(TMath::Pi()/2. - trackhit_.beta) << endl;
0193             }
0194             else{ 
0195                 h_sizex_alpha->Fill( (trackhit_.alpha + TMath::Pi())*180. / TMath::Pi(),clust_.size_x);
0196                 h_sizex_alpha_cotan->Fill(TMath::Tan(TMath::Pi()/2. - trackhit_.alpha),clust_.size_x);
0197                 h_sizey_beta_cotan->Fill(TMath::Tan(TMath::Pi()/2. - trackhit_.beta),clust_.size_y);
0198                 h_alpha_cotan->Fill(TMath::Tan(TMath::Pi()/2. - trackhit_.alpha));
0199                 h_beta_cotan->Fill(TMath::Tan(TMath::Pi()/2. - trackhit_.beta));
0200                 fAngles << TMath::Tan(TMath::Pi()/2. - trackhit_.alpha) << "\t" << TMath::Tan(TMath::Pi()/2. - trackhit_.beta) << endl;
0201             }
0202             h_alpha->Fill( trackhit_.alpha*180. / TMath::Pi());
0203             h_beta->Fill( trackhit_.beta*180. / TMath::Pi());
0204             h_alpha_beta->Fill( trackhit_.alpha*180. / TMath::Pi(), trackhit_.beta*180. / TMath::Pi());
0205             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));
0206         } 
0207     }
0208     
0209     cout << TMath::Pi()/2. << endl;
0210     TH1F * h_mean = new TH1F("h_mean","h_mean", anglebins, 0, 180);
0211     TH1F * h_slice_ = new TH1F("h_slice","h_slice", 10, .5, 10.5);
0212     TH1F * h_mean_cotan_alpha = new TH1F("h_mean_cotan_alpha","h_mean_contan_alpha", anglebinscotan, -3, 3);
0213     TH1F * h_slice_cotan_alpha = new TH1F("h_slice_cotan_alpha","h_slice_cotan_alpha", 10, .5, 10.5);
0214     TH1F * h_mean_cotan_beta = new TH1F("h_mean_cotan_beta","h_mean_contan_beta", anglebinscotan, -3, 3);
0215     TH1F * h_slice_cotan_beta = new TH1F("h_slice_cotan_beta","h_slice_cotan_beta", 10, .5, 10.5);
0216     //loop over bins in depth (z-local-coordinate) (in order to fit slices)

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

0219         h_slice_->Reset("ICE");
0220         
0221         // determine sigma and sigma^2 of the adc counts and average adc counts

0222             //loop over bins in drift width

0223         for( int j = 1; j<= 10; j++){
0224             h_slice_->SetBinContent(j,h_sizex_alpha->GetBinContent(i,j));
0225         } // end loop over bins in drift width

0226             
0227         double mean = h_slice_->GetMean(1); 
0228       double error = h_slice_->GetMeanError(1);
0229             
0230         h_mean->SetBinContent(i, mean);
0231         h_mean->SetBinError(i, error);  
0232     }// end loop over bins in depth 

0233     
0234     for( int i = 1; i <= anglebinscotan; i++){
0235         h_slice_cotan_alpha->Reset("ICE");
0236         h_slice_cotan_beta->Reset("ICE");
0237         //loop over bins in drift width

0238         for( int j = 1; j<= 10; j++){
0239             h_slice_cotan_alpha->SetBinContent(j,h_sizex_alpha_cotan->GetBinContent(i,j));
0240             h_slice_cotan_beta->SetBinContent(j,h_sizey_beta_cotan->GetBinContent(i,j));
0241         } // end loop over bins in drift width

0242             
0243         double mean = h_slice_cotan_alpha->GetMean(1); 
0244         double error = h_slice_cotan_alpha->GetMeanError(1);
0245             
0246         h_mean_cotan_alpha->SetBinContent(i, mean);
0247         h_mean_cotan_alpha->SetBinError(i, error);  
0248         
0249         double mean_beta = h_slice_cotan_beta->GetMean(1); 
0250         double error_beta = h_slice_cotan_beta->GetMeanError(1);
0251             
0252         h_mean_cotan_beta->SetBinContent(i, mean_beta);
0253         h_mean_cotan_beta->SetBinError(i, error_beta);  
0254     }// end loop over bins in depth 

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

0273     //c2->cd(1);

0274     h_sizex_alpha_cotan->GetXaxis()->SetTitle("cotan(#alpha)");
0275     h_sizex_alpha_cotan->GetYaxis()->SetTitle("cluster size x[pixels]");
0276     h_sizex_alpha_cotan->Draw("col");
0277     //c2->cd(2);    

0278     h_mean_cotan_alpha->GetXaxis()->SetTitle("cotan(#alpha)");
0279     h_mean_cotan_alpha->GetYaxis()->SetTitle("transverse cluster size [pixels]");
0280     //h_mean_cotan_alpha->Draw();

0281 
0282     
0283     TCanvas * c3 = new TCanvas("c3", "c3", 600, 600);
0284     h_mean_cotan_alpha->Draw();
0285     h_mean_cotan_alpha->Fit(func,"ERQ");
0286     double correlation = h_alpha_beta_cotan->GetCorrelationFactor();
0287     cout << "corrlation between cotan(alpha) and cotan(beta): " << correlation << endl;
0288     cout << "number of hits used: " << nhits << endl;
0289     
0290     TCanvas * c8 = new TCanvas("c8", "c8", 600, 600);
0291     h_mean_cotan_beta->GetXaxis()->SetTitle("cotan(#beta)");
0292     h_mean_cotan_beta->GetYaxis()->SetTitle("longitudinal cluster size [pixels]");
0293     h_mean_cotan_beta->Draw();
0294     h_mean_cotan_beta->Fit(func_beta,"ERQ");
0295     
0296     TCanvas * c4 = new TCanvas("c4", "c4", 600, 600);
0297     h_alpha_cotan->Draw();
0298     
0299     TCanvas * c5 = new TCanvas("c5", "c5", 600, 600);
0300     h_beta_cotan->Draw();
0301     
0302     TCanvas * c6 = new TCanvas("c6", "c6", 600, 600);
0303     h_alpha->GetXaxis()->SetTitle("#alpha [^{o}]");
0304     h_alpha->GetYaxis()->SetTitle("number of hits");
0305     h_alpha->Draw();
0306     
0307     
0308     TCanvas * c7 = new TCanvas("c7", "c7", 600, 600);
0309     h_beta->GetXaxis()->SetTitle("#beta [^{o}]");
0310     h_beta->GetYaxis()->SetTitle("number of hits");
0311     h_beta->Draw();
0312     
0313     TCanvas * c9 = new TCanvas("c9", "c9", 600, 600);
0314     h_alpha_beta->GetXaxis()->SetTitle("#alpha [^{o}]");
0315     h_alpha_beta->GetYaxis()->SetTitle("#beta [^{o}]");
0316     h_alpha_beta->Draw("col");
0317     
0318     TCanvas * c10 = new TCanvas("c10", "c10", 600, 600);
0319     h_chi2->GetXaxis()->SetTitle("#chi^{2}/ndof");
0320     h_chi2->GetYaxis()->SetTitle("number of hits");
0321     h_chi2->Draw();
0322     
0323     TCanvas * c11 = new TCanvas("c11", "c11", 600, 600);
0324     h_charge->GetXaxis()->SetTitle("charge [1000 e^{-}]");
0325     h_charge->GetYaxis()->SetTitle("number of hits");
0326     h_charge->Draw();
0327     return 0;
0328 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!