Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:14:32

0001 {
0002 //  setTDRStyle();

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 //  ofstream fLorentzFit( "lorentzFit.txt", ios::trunc );

0021 //  fLorentzFit.precision( 4 );

0022 //  fLorentzFit << "module" << "\t" << "layer" << "\t" << "offset" << "\t" << "error" << "\t" << "slope" << "\t" << "error" << "\t" "rel.err" << "\t" "pull" << "\t" << "chi2" << "\t" << "prob" << endl;

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 // fill the histrograms with the ntpl

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_); // charge is given in 1000 e

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 //loop over bins in depth (z-local-coordinate) (in order to fit slices)

0129 for( int i = 1; i <= hist_depth_; i++){
0130   //                findMean(i, (i_module + (i_layer - 1) * 8));

0131   double npix = 0;
0132 
0133   h_drift_depth_adc_slice_->Reset("ICE");
0134         
0135   // determine sigma and sigma^2 of the adc counts and average adc counts

0136   //loop over bins in drift width

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   } // end loop over bins in drift width

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 }// end loop over bins in depth 

0158     
0159 TCanvas * c1 = new TCanvas("c1", "c1", 1200, 600);
0160 c1->Divide(2,1);
0161 c1->cd(1);
0162 //  h_drift_depth_noadc->

0163 h_drift_depth_noadc->Draw("colz");
0164 //  c1->cd(2);

0165 //  h_drift_depth_adc->Draw();

0166 c1->cd(2);
0167 h_mean->Draw();
0168 //  c1->cd(4);

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 //  delete h_mean;

0179 //  delete h_drift_depth_adc_slice_;

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 }