Attachment 'CalibrationFunctions.C'

Download

   1 //###################################################################
   2 // Determine Reconstruction efficiency and Energy calibration
   3 // functions to be applied to theoretical background expectation (RecEff)
   4 // and Data samples (Energy calibration)
   5 //
   6 // Reconstruction efficiency: Look if a ProcessedReco photon 
   7 // connected to MC photon via MCTruth is found within 0.1 rad (~5.7 deg)
   8 //
   9 // ErecVsEgen:  Ratio of highest Processed Photon energy to MC energy.
  10 //              Processed photon has to be connected to MC photon
  11 //              and within 0.1 rad. Only MC photons not converted before
  12 //              ECAL (or before end of tracking) are taken into account
  13 //              Due to merging there should be only one Processed photon
  14 //              connected to MC photon in 0.1 rad 
  15 //###################################################################
  16 //###################################################################
  17 
  18 # include "GlobalVariables.C"
  19 
  20 // ########### CUTS
  21 
  22   double energy_cut = 2.0;
  23   double cosTheta_cut[2] = {0.0,0.995};
  24   int conv_cut = 1;
  25 
  26 void CalibrationFunctions(double e_min, double cos_min, double cos_max, int conv){
  27   
  28   energy_cut = e_min;
  29   cosTheta_cut[0] = cos_min;
  30   cosTheta_cut[1] = cos_max;
  31   conv_cut = conv;
  32 
  33   GammaReconstructonEfficiency();
  34   ErecVsEgen();
  35 
  36   gROOT->ProcessLine(".q");
  37 
  38   // return;
  39 }
  40 
  41 
  42 void GammaReconstructonEfficiency(){
  43 
  44   printf("Reconstruction efficiencies\n");
  45 
  46   // ############  OPEN OUTPUTFILE
  47   TFile* f_out = new TFile("../data/Auxiliary/ReconstructionEfficiencies.root","RECREATE");
  48 
  49   // ############  RECONSTRUCTION EFFICIENCY HISTOGRAMS IN COS(THETA), THETA, E AND PT
  50 
  51 
  52   TCanvas* c1 = new TCanvas("c1","Reconstruction efficiency in cos(#Theta)");
  53   TCanvas* c2 = new TCanvas("c2","Reconstruction efficiency in #Theta");
  54   TCanvas* c3 = new TCanvas("c3","Reconstruction efficiency in Energy");
  55   TCanvas* c4 = new TCanvas("c4","Reconstruction efficiency in PT");
  56   TCanvas* c5 = new TCanvas("c5","Reconstruction efficiency in E and cos(#Theta)");
  57 
  58   TProfile* p1 = new TProfile("RecEff_Cos","Reconstruction efficiency;cos(#Theta_{gen,#gamma});Efficiency",100,cosTheta_cut[0],cosTheta_cut[1]);
  59   TProfile* p2 = new TProfile("RecEff_Theta","Reconstruction efficiency;#Theta_{gen,#gamma} [rad];Efficiency",100,0.0,3.2);
  60   TProfile* p3 = new TProfile("RecEff_E","Reconstruction efficiency;E_{gen,#gamma} [GeV];Efficiency",250,0,250);
  61   TProfile* p4 = new TProfile("RecEff_PT","Reconstruction efficiency;P_{T,#gamma} [GeV];Efficiency",250,0,250);
  62   TProfile2D* p5 = new TProfile2D("RecEff_E_Cos","Reconstruction efficiency;E_{gen,#gamma} [GeV]; cos(#Theta_{gen,#gamma});Efficiency",25,0,250,10,cosTheta_cut[0],cosTheta_cut[1]);
  63 
  64 
  65   p1->GetYaxis()->SetRangeUser(0.8,1.05);
  66   p2->GetYaxis()->SetRangeUser(0.8,1.05);
  67   p3->GetYaxis()->SetRangeUser(0.8,1.05);
  68   p4->GetYaxis()->SetRangeUser(0.8,1.05);
  69   p5->GetZaxis()->SetRangeUser(0.8,1.05);
  70 
  71 
  72   // ##########  OPEN INPUT FILES
  73 
  74   TChain* tree  = new TChain("Tree");
  75   tree->Add("../data/Preprocessed/ILD00_n1n1a_Pe-1.0_Pp1.0Tree_1.root");
  76   tree->Add("../data/Preprocessed/ILD00_n1n1a_Pe-1.0_Pp1.0Tree_2.root");
  77   
  78   LoadTree(tree);
  79 
  80   TChain* proctree  = new TChain("ProcTree");
  81   proctree->Add("../data/Preprocessed/ILD00_n1n1a_Pe-1.0_Pp1.0Tree_1.root");
  82   proctree->Add("../data/Preprocessed/ILD00_n1n1a_Pe-1.0_Pp1.0Tree_2.root");
  83   LoadPreProcTree(proctree);
  84 
  85 
  86   // ##########  FILL HISTOGRAMS
  87 
  88   int nEntries = proctree->GetEntries();
  89   // nEntries = 1e5;
  90 
  91   for ( int i(0); i < nEntries; i++ ){
  92     tree->GetEntry(i);
  93     proctree->GetEntry(i);
  94     for ( int j(0); j < nMCParticles; j++ ){
  95       // Check if MC is photon and apply cuts
  96       if ( MC_type[j] == 22 && MC_isReal[j] && MC_energy[j] > energy_cut &&
  97 	   TMath::Abs(MC_cosTheta[j]) > cosTheta_cut[0] && TMath::Abs(MC_cosTheta[j]) < cosTheta_cut[1] ){  
  98 	int counter(0);
  99 	for ( int k(0); k < nProcessedRecos; k++ ){
 100 	  if (  ProcessedReco_type[k] == 22 ){
 101 	    double angle = TMath::ACos((MC_px[j]*ProcessedReco_px[k]+MC_py[j]*ProcessedReco_py[k]
 102 					+MC_pz[j]*ProcessedReco_pz[k])/(MC_absP[j]*ProcessedReco_absP[k]));
 103 	    if (angle < 0.0) cout << "ANGLE < 0.0" << endl;
 104 	    
 105 	    if ( MC_id[j] == ProcessedRecoRelatedToId[k] && angle < 0.1){
 106 	      counter++;
 107 	      break;            // break as soon as one PFlow photon is found 
 108 	    }
 109 	  }
 110 	}
 111 	
 112 	double PT = MC_energy[j]*TMath::Sin(TMath::ACos(MC_cosTheta[j]));
 113 
 114 	// Fill for each "real" MC photon, counter is either 0 or 1
 115 	p1->Fill(TMath::Abs(MC_cosTheta[j]),counter);
 116 	p2->Fill(TMath::ACos(MC_cosTheta[j]),counter);
 117 	p3->Fill(MC_energy[j],counter);
 118 	p4->Fill(PT,counter);
 119 	p5->Fill(MC_energy[j],TMath::Abs(MC_cosTheta[j]),counter);
 120       }
 121     }
 122   }
 123 
 124 
 125   // ########## DRAW
 126 
 127 
 128 
 129   TLatex latex;
 130   latex.SetNDC();
 131   latex.SetTextSize(0.05);
 132   latex.SetTextAlign(13);
 133 
 134   c1->cd();
 135   p1->DrawCopy();
 136   latex.DrawLatex(0.3,0.3,Form("Total Efficiency: %.3f",p1->GetMean(2)));
 137 
 138   c2->cd();
 139   p2->DrawCopy();
 140   c2->cd();
 141   latex.DrawLatex(0.3,0.3,Form("Total Efficiency: %.3f",p2->GetMean(2)));
 142 
 143   c3->cd();
 144   p3->DrawCopy();
 145   latex.DrawLatex(0.3,0.3,Form("Total Efficiency: %.3f",p3->GetMean(2)));
 146 
 147   c4->cd();
 148   p4->DrawCopy();
 149   latex.DrawLatex(0.3,0.3,Form("Total Efficiency: %.3f",p4->GetMean(2)));
 150 
 151   c5->cd();
 152   p5->DrawCopy();
 153 
 154 
 155   // ############ PRINT TO FILE
 156 
 157   f_out->cd();
 158 
 159   p1->Write();
 160   p2->Write();
 161   p3->Write();
 162   p4->Write();
 163   p5->Write();
 164   
 165   c1->Print(Form("Reconstruction_Efficiency_cosTheta.eps"));
 166   c1->Print(Form("Reconstruction_Efficiency_cosTheta.png"));
 167 
 168   
 169   c2->Print(Form("Reconstruction_Efficiency_Theta.eps"));
 170   c2->Print(Form("Reconstruction_Efficiency_Theta.png"));
 171   
 172   c3->Print(Form("Reconstruction_Efficiency_Energy.eps"));
 173   c3->Print(Form("Reconstruction_Efficiency_Energy.png"));
 174   
 175   c4->Print(Form("Reconstruction_Efficiency_PT.eps"));
 176   c4->Print(Form("Reconstruction_Efficiency_PT.png"));
 177     
 178   c5->Print(Form("Reconstruction_Efficiency_Energy_cosTheta.eps"));
 179   c5->Print(Form("Reconstruction_Efficiency_Energy_cosTheta.png"));
 180 
 181   //  
 182   f_out->Close();
 183   return;
 184 }
 185 
 186 
 187 // #########################################################################################################
 188 // #########################################################################################################
 189 
 190 
 191 
 192 void ErecVsEgen(){
 193 
 194   printf("Energy calibration\n");
 195 
 196   // ########  OPEN OUTPUTFILE
 197 
 198   TFile* f_out = new TFile("../data/Auxiliary/Energycalibration.root","RECREATE");
 199 
 200 
 201   char str[128];
 202 
 203   // ######## CUTS
 204 
 205   double z_cut,r_cut;
 206   
 207   if (conv_cut == 1){
 208     sprintf(str,"InECAL"); 
 209     z_cut = 2450.0;
 210     r_cut = 1843.3;
 211   }
 212   if (conv_cut == 2){
 213     sprintf(str,"AfterTracking"); 
 214     z_cut = 2246.0; // cut at tpc dimensions
 215     r_cut = 1739.0;
 216   }
 217 
 218   // ############  ENERGY RATIO HISTOGRAMS IN COS(THETA), THETA, E AND PT
 219 
 220   TCanvas* c1 = new TCanvas("c2","Ratio E_rec/E_gen vs cos(Theta_gen)");
 221   TCanvas* c2 = new TCanvas("c3","Ratio E_rec/E_gen vs Theta_gen");
 222   TCanvas* c3 = new TCanvas("c1","Ratio E_rec/E_gen vs E_gen");
 223   TCanvas* c4 = new TCanvas("c4","Ratio E_rec/E_gen vs PT_gen");
 224   TCanvas* c5 = new TCanvas("c5","Ratio E_rec/E_gen vs PT_gen");
 225 
 226   TProfile* p1 = new TProfile("ErecEgen_Cos","Ratio E_rec/E_gen vs cos(Theta_gen);cos(#Theta_{gen});#bar{E_{rec}/E_{gen}}",150,0,1);
 227   TProfile* p2 = new TProfile("ErecEgen_Theta","Ratio E_rec/E_gen vs Theta_gen;#Theta_{gen};#bar{E_{rec}/E_{gen}}",150,0,3.2);
 228   TProfile* p3 = new TProfile("ErecEgen_E","Ratio E_rec/E_gen vs E_gen;E_{gen} [GeV];#bar{E_{rec}/E_{gen}}",250,0,250);
 229   TProfile* p4 = new TProfile("ErecEgen_PT","Ratio E_rec/E_gen vs PT_gen;P_{T} [GeV];#bar{E_{rec}/E_{gen}}",250,0,250);
 230   TProfile2D* p5 = new TProfile2D("ErecEgen_E_Cos","Ratio E_rec/E_gen vs E and Cos;E_{gen} [GeV];cos(#Theta_{gen});#bar{E_{rec}/E_{gen}}",25,0,250,10,0.0,1.0);
 231 
 232   p1->GetYaxis()->SetRangeUser(0.8,1.05);
 233   p2->GetYaxis()->SetRangeUser(0.8,1.05);
 234   p3->GetYaxis()->SetRangeUser(0.8,1.05);
 235   p4->GetYaxis()->SetRangeUser(0.8,1.05);
 236   p5->GetZaxis()->SetRangeUser(0.8,1.05); 
 237 
 238 
 239   // ########## OPEN INPUT FILES
 240 
 241   TChain* tree  = new TChain("Tree");
 242   tree->Add("../data/Preprocessed/ILD00_n1n1a_Pe-1.0_Pp1.0Tree_1.root");
 243   tree->Add("../data/Preprocessed/ILD00_n1n1a_Pe-1.0_Pp1.0Tree_2.root"); 
 244 
 245   LoadTree(tree);
 246 
 247   TChain* proctree  = new TChain("ProcTree");
 248   proctree->Add("../data/Preprocessed/ILD00_n1n1a_Pe-1.0_Pp1.0Tree_1.root");
 249   proctree->Add("../data/Preprocessed/ILD00_n1n1a_Pe-1.0_Pp1.0Tree_2.root");
 250   LoadPreProcTree(proctree);
 251 
 252  
 253 
 254   // ####### FILL HISTOGRAMS
 255 
 256   int nEntries = tree->GetEntries();
 257   //nEntries = 1e5;
 258 
 259   for ( int i(0); i < nEntries; i++ ){
 260     tree->GetEntry(i);
 261     proctree->GetEntry(i);
 262     for ( int j(0); j < nMCParticles; j++ ){
 263       if ( MC_type[j] == 22 && MC_isReal[j] && MC_energy[j] > energy_cut &&
 264 	   TMath::Abs(MC_cosTheta[j]) > cosTheta_cut[0] && TMath::Abs(MC_cosTheta[j]) < cosTheta_cut[1] ){ 
 265 
 266 	// cut conversions
 267 	
 268 	double z = TMath::Abs(MC_endpoint_z[j]);
 269 	double r = TMath::Sqrt(MC_endpoint_x[j]*MC_endpoint_x[j] + MC_endpoint_y[j]*MC_endpoint_y[j]);
 270 	
 271 	if ( conv_cut && (z < z_cut && r < r_cut) )
 272 	  continue;
 273 	
 274 	// find max energy reco gamma
 275 	double maxEnergy(0.0);
 276 	int Recoid(-1);
 277 	for ( int k(0); k < nProcessedRecos; k++ ){   // Find PFlow photon with highest energy for MC photon
 278 	  if (  ProcessedReco_type[k] == 22 && MC_id[j] == ProcessedRecoRelatedToId[k] && ProcessedReco_energy[k] > maxEnergy){
 279 	    maxEnergy =  ProcessedReco_energy[k];
 280 	    Recoid = k;
 281 	  }
 282 	}
 283 	if (Recoid != -1){
 284 	  double angle = TMath::ACos((MC_px[j]*ProcessedReco_px[Recoid]+MC_py[j]*ProcessedReco_py[Recoid]
 285 				      +MC_pz[j]*ProcessedReco_pz[Recoid])/(MC_absP[j]*ProcessedReco_absP[Recoid]));
 286 	  if (angle < 0.1){
 287 	    double ratio = ProcessedReco_energy[Recoid]/MC_energy[j];
 288 	    p1->Fill(TMath::Abs(MC_cosTheta[j]),ratio);
 289 	    p2->Fill(TMath::ACos(MC_cosTheta[j]),ratio);
 290 	    p3->Fill(MC_energy[j],ratio);
 291 	    p4->Fill(TMath::Sin(TMath::ACos(MC_cosTheta[j]))*MC_energy[j],ratio);
 292 	    p5->Fill(MC_energy[j],MC_cosTheta[j],ratio);
 293 	  }
 294 	}
 295       }
 296     }
 297   }
 298 
 299   // ######## FIT
 300 
 301   TF1* func = new TF1("EnergyCalibration_Func",fitfunc,0.0,1.0,13);
 302 
 303   func->SetParameter(0,1.0);
 304 
 305   func->SetParameter(1,-0.1);
 306   func->SetParameter(2,0.23);
 307   func->SetParameter(3,0.01);
 308 
 309   func->SetParameter(4,-0.1);
 310   func->SetParameter(5,0.58);
 311   func->SetParameter(6,0.01);
 312 
 313   func->SetParameter(7,-0.1);
 314   func->SetParameter(8,0.77);
 315   func->SetParameter(9,0.01);
 316 
 317   func->SetParameter(10,-0.1);
 318   func->SetParameter(11,0.98);
 319   func->SetParameter(12,0.01);
 320 
 321   // ##########  DRAW
 322 
 323   c1->cd();
 324   p1->Fit(func,"R","");
 325   p1->DrawCopy();
 326 
 327   c2->cd();
 328   p2->DrawCopy();
 329 
 330   c3->cd();
 331   p3->DrawCopy();
 332 
 333   c4->cd();
 334   p4->DrawCopy();
 335 
 336   c5->cd();
 337   p5->DrawCopy();
 338   // ########### PRINT
 339 
 340   f_out->cd();
 341 
 342   p1->Write();
 343   p2->Write();
 344   p3->Write();
 345   p4->Write();
 346   p5->Write();
 347   func->Write();
 348 
 349   f_out->Close();
 350 
 351  
 352   c1->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_cosTheta_%s.eps",str));
 353   c1->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_cosTheta_%s.png",str));
 354   //c1->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_cosTheta_%s.root",str));
 355   
 356   c2->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_Theta_%s.eps",str));
 357   c2->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_Theta_%s.png",str));
 358 
 359   c3->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_Egen_%s.eps",str));
 360   c3->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_Egen_%s.png",str));
 361 
 362   c4->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_PT_%s.eps",str));
 363   c4->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_PT_%s.png",str));
 364 
 365   c5->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_Egen_cosTheta_%s.eps",str));
 366   c5->Print(Form("Ratio_PFlow_ProcessedPhotonenergy_MC_Energy_Vs_Egen_cosTheta_%s.png",str)); 
 367   return;
 368 }
 369 
 370 
 371 // ########### FITFUNCTION
 372 
 373 double constant(double* x, double* par){
 374   return par[0];
 375 }
 376 
 377 double gaussdip(double* x, double* par){
 378   return par[0]/(sqrt(2*TMath::Pi())*par[2])*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/(2*par[2]*par[2]));
 379 }
 380 
 381 double fitfunc(double* x, double* par){
 382   return constant(x,par) + gaussdip(x,&par[1]) + gaussdip(x,&par[4]) + gaussdip(x,&par[7]) + gaussdip(x,&par[10]);
 383 }

Attached Files

To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.
  • [get | view] (2012-01-30 13:28:08, 4.3 KB) [[attachment:CalibrateData.C]]
  • [get | view] (2012-01-30 13:27:40, 12.0 KB) [[attachment:CalibrationFunctions.C]]
  • [get | view] (2014-06-20 09:25:03, 5.2 KB) [[attachment:CreateDataBackground.C]]
  • [get | view] (2014-06-20 09:21:28, 22.3 KB) [[attachment:CreateDataForSignalAndBackground.C]]
  • [get | view] (2014-06-20 09:24:20, 0.2 KB) [[attachment:CreateDataForSignalAndBackground.py]]
  • [get | view] (2012-01-30 12:32:43, 12.0 KB) [[attachment:FindMergeAngle.C]]
  • [get | view] (2012-01-30 13:22:33, 9.8 KB) [[attachment:PreprocessData.C]]
  • [get | view] (2012-01-30 13:51:59, 16.8 KB) [[attachment:Selection.C]]
  • [get | view] (2014-06-20 09:26:52, 13.2 KB) [[attachment:TSysLimit_searches.C]]
  • [get | view] (2014-06-20 09:28:25, 14.9 KB) [[attachment:TSysLimit_searches_onebin.C]]
  • [get | view] (2014-06-20 09:25:07, 12.9 KB) [[attachment:WIMP_functions.C]]
  • [get | view] (2014-06-20 12:52:10, 1.3 KB) [[attachment:runLambdaSearch.sh]]
  • [get | view] (2014-06-20 12:52:16, 0.2 KB) [[attachment:submitter.sh]]
 All files | Selected Files: delete move to page copy to page

You are not allowed to attach a file to this page.