Attachment 'FindMergeAngle.C'

Download

   1 #include "GlobalVariables.C"
   2 
   3 const double PI = TMath::Pi();
   4 
   5 void FindMergeAngle(double e_min, double e_max, double costheta_min, double costheta_max, int tracking_cut){
   6 
   7   char filename[256];
   8   char line[128];
   9   sprintf(filename,"../results/MergeAngle_Emin%.1f_Emax%.1f_Cosmin%.3f_Cosmax%.1f_cut%i.root",
  10 	  e_min, e_max, costheta_min, costheta_min, tracking_cut);
  11   TFile* f_out = new TFile(filename,"RECREATE");
  12 
  13   double z_cut,r_cut;
  14   if (tracking_cut == 1){
  15     z_cut = 2246.0; // cut at tpc dimensions
  16     r_cut = 1739.0;
  17   }
  18   if (tracking_cut == 2){
  19     z_cut = 2450.0;
  20     r_cut = 1843.3;
  21   }
  22 
  23 
  24   double MatchCounter[18][100][4] = {0.0};  // For all 18 background files and all 100 cone angles count
  25                                             // FoundTrueMatches, TrueMatches, and FoundMatches.
  26                                             // Example: for file 1 field [0][2][0] holds the value of the third
  27                                             // tested cone value and fields [0][2][1],[0][2][2] and [0][2][3]
  28                                             // hold the corresponding values of FoundTrueMatches, TrueMatches, and FoundMatches.
  29 
  30 
  31   double GlobalMatchCounter[3][100][4] = {0.0};  // similar to Matchcounter, however only processes
  32                                                  // with one, two or three photons are concerned (a,aa,aaa).
  33 
  34 
  35   TTree* matchtree = new TTree("MatchTree","MatchTree");
  36   matchtree->Branch("MatchCounter",&MatchCounter,"MatchCounter[18][100][4]/D");
  37   matchtree->Branch("GlobalMatchCounter",&GlobalMatchCounter,"GlobalMatchCounter[3][100][4]/D");
  38 
  39 
  40 
  41   // out txt files with efficiencies and purities
  42   ofstream txtfile(Form("../results/Matching_Purity_Efficiency.txt"));
  43   ofstream globaltxtfile(Form("../results/Total_Matching_Purity_Efficiency.txt"));  
  44 
  45  
  46   FILE* file = fopen("../steeringfiles/ILD00_Processed_nxnxa_RootFiles.txt","r");
  47   TFile* f_in;
  48   TTree* tree;
  49 
  50   int n = 0;
  51   while(fscanf(file,"%s\n",filename) != EOF) {
  52     f_in = new TFile(filename,"OPEN");
  53     tree = (TTree*)f_in->Get("Tree");
  54     LoadTree(tree);
  55 
  56     // Fill MatchCounter array
  57     tree->GetEntry(0);
  58     cout << Process << endl;
  59 
  60     for ( int i(0); i < 2000; i++ ){ // for all entries in tree, if set by hand should not exeed 2100
  61       // which is the number of events in n3n3aaa files
  62       tree->GetEntry(i);
  63       double coneangle = 0.0;
  64       for ( int l(0); l < 100; l++){     // for all cone angles to evaluate
  65 
  66 	// write cone angle
  67 	MatchCounter[n][l][0] = coneangle;
  68 
  69 	for ( int j(0); j < nMCParticles; j++ ){
  70 	  if (MC_isReal[j] && MC_type[j] == 22 && MC_energy[j] > e_min && MC_energy[j] < e_max &&
  71 	      TMath::Abs(MC_cosTheta[j]) > costheta_min && TMath::Abs(MC_cosTheta[j]) < costheta_max){
  72 	    for ( int k(0); k < nRecoParticles; k++ ){
  73 	      if (Reco_type[k] == 22){
  74 		double angle = TMath::ACos((MC_px[j]*Reco_px[k]+MC_py[j]*Reco_py[k]
  75 					    +MC_pz[j]*Reco_pz[k])/(MC_absP[j]*Reco_absP[k]));
  76 		if (angle < coneangle){
  77 		  MatchCounter[n][l][3]++;     // increase FoundMatches
  78 		  if ( RecoToMCId[k] == MC_id[j] ){
  79 		    MatchCounter[n][l][1]++;    // increase FoundTrueMatches
  80 		  }
  81 		}
  82 		if ( RecoToMCId[k] == MC_id[j] ){
  83 		  MatchCounter[n][l][2]++;      // increase TrueMatches
  84 		}
  85 	      }
  86 	    }
  87 	  }
  88 	}
  89 	coneangle += 0.002;
  90       }
  91     }
  92     delete tree;
  93     n++;
  94   }
  95   fclose(file);
  96       
  97   
  98   // Calculate efficiencies for all subprocesses and polarisations and fill GlobalMatchCounter
  99   file = fopen("../steeringfiles/ILD00_Processed_nxnxa_RootFiles.txt","r");
 100   n = 0;
 101   while(fscanf(file,"%s\n",filename) != EOF) {
 102     f_in = new TFile(filename,"OPEN");
 103     tree  = (TTree*)f_in->Get("Tree");
 104     LoadTree(tree);
 105     tree->GetEntry(0);
 106 
 107     sprintf(line,"##################\n\n");
 108     txtfile.write(line,strlen(line));
 109     sprintf(line,"Process:                      %s\n",Process);
 110     txtfile.write(line,strlen(line));
 111     sprintf(line,"P_e:                          %i\nP_p:                          %i\n\n",P_e,P_p);
 112     txtfile.write(line,strlen(line));
 113 
 114     int s;
 115     if (strlen(Process) == 5) s = 0;  // for nxnxa
 116     if (strlen(Process) == 6) s = 1;  // for nxnxaa
 117     if (strlen(Process) == 7) s = 2;  // for nxnxaaa
 118 
 119     for ( int l(0); l < 100; l++){     // for all cone angles to evaluate
 120       double coneangle = MatchCounter[n][l][0];
 121       double FoundTrueMatches = MatchCounter[n][l][1];
 122       double TrueMatches =  MatchCounter[n][l][2];
 123       double FoundMatches = MatchCounter[n][l][3];
 124 
 125       GlobalMatchCounter[s][l][0] = coneangle;
 126       GlobalMatchCounter[s][l][1] += FoundTrueMatches;
 127       GlobalMatchCounter[s][l][2] += TrueMatches;
 128       GlobalMatchCounter[s][l][3] += FoundMatches;
 129  
 130     }
 131     n++;
 132   }
 133   
 134   matchtree->Fill();
 135 
 136   f_out->cd();
 137   matchtree->Write();
 138   f_out->Close();
 139 
 140   MakePlots( e_min,  e_max, costheta_min, costheta_max, tracking_cut);
 141 
 142   return;
 143 }
 144     
 145 
 146 void MakePlots(double e_min, double e_max, double costheta_min, double costheta_max, int tracking_cut){
 147 
 148 
 149   double MatchCounter[18][100][4] = {0.0};  // For all 18 background files and all 100 cone angles count
 150                                             // FoundTrueMatches, TrueMatches, and FoundMatches.
 151                                             // Example: for file 1 field [0][2][0] holds the value of the third
 152                                             // tested cone value and fields [0][2][1],[0][2][2] and [0][2][3]
 153                                             // hold the corresponding values of FoundTrueMatches, TrueMatches, and FoundMatches.
 154 
 155 
 156   double GlobalMatchCounter[3][100][4] = {0.0};  // similar to Matchcounter, however only processes
 157                                                  // with one, two or three photons are concerned (a,aa,aaa)
 158 
 159   char filename[256];
 160   sprintf(filename,"../results/MergeAngle_Emin%.1f_Emax%.1f_Cosmin%.3f_Cosmax%.1f_cut%i.root",
 161 	  e_min, e_max, costheta_min, costheta_min, tracking_cut);
 162   TFile* f_in = new TFile(filename,"OPEN");
 163   TTree* matchtree = (TTree*) f_in->Get("MatchTree");
 164 
 165   matchtree->SetBranchAddress("MatchCounter",&MatchCounter);
 166   matchtree->SetBranchAddress("GlobalMatchCounter",&GlobalMatchCounter);
 167   matchtree->GetEntry(0);
 168 
 169 
 170   TCanvas* c1 = new TCanvas("c1","n1n1 P_e:-1.0 P_p: 1.0");
 171   TCanvas* c2 = new TCanvas("c2","n1n1 P_e: 1.0 P_p:-1.0");
 172   TCanvas* c3 = new TCanvas("c3","n2n2 P_e:-1.0 P_p: 1.0");
 173   TCanvas* c4 = new TCanvas("c4","n2n2 P_e: 1.0 P_p:-1.0");
 174   TCanvas* c5 = new TCanvas("c5","n3n3 P_e:-1.0 P_p: 1.0");
 175   TCanvas* c6 = new TCanvas("c6","n3n3 P_e: 1.0 P_p:-1.0");
 176   TCanvas* c7 = new TCanvas("c7","All Processes");
 177   
 178   TH2D* h1 = new TH2D("h1","n1n1 P_e:-1.0 P_p: 1.0;Efficiency;Purity",1000,0.9,1.02,1000,0.9,1.02);
 179   TH2D* h2 = new TH2D("h2","n1n1 P_e: 1.0 P_p:-1.0;Efficiency;Purity",1000,0.9,1.02,1000,0.9,1.02);
 180   TH2D* h3 = new TH2D("h3","n2n2 P_e:-1.0 P_p: 1.0;Efficiency;Purity",1000,0.9,1.02,1000,0.9,1.02);
 181   TH2D* h4 = new TH2D("h4","n2n2 P_e: 1.0 P_p:-1.0;Efficiency;Purity",1000,0.9,1.02,1000,0.9,1.02);
 182   TH2D* h5 = new TH2D("h5","n3n3 P_e:-1.0 P_p: 1.0;Efficiency;Purity",1000,0.9,1.02,1000,0.9,1.02);
 183   TH2D* h6 = new TH2D("h6","n3n3 P_e: 1.0 P_p:-1.0;Efficiency;Purity",1000,0.9,1.02,1000,0.9,1.02);
 184   TH2D* h7 = new TH2D("h7","All Processes;Efficiency;Purity",1000,0.9,1.02,1000,0.9,1.02);
 185 
 186   h7->SetLineWidth(2);
 187 
 188   c7->cd();
 189   h7->Draw();
 190 
 191   for ( int n(0); n < 18; n++ ){         // for all 18 infiles
 192     for ( int l(0); l < 100; l++){     // for all cone angles to evaluate
 193   
 194       double coneangle = MatchCounter[n][l][0];
 195       double FoundTrueMatches = MatchCounter[n][l][1];
 196       double TrueMatches =  MatchCounter[n][l][2];
 197       double FoundMatches = MatchCounter[n][l][3];
 198 
 199       if (TrueMatches != 0 && FoundMatches != 0){
 200 	double efficiency = FoundTrueMatches/TrueMatches;
 201 	double purity = FoundTrueMatches/FoundMatches;
 202 	
 203 	// fill histo
 204 	if ( n == 0 || n == 2 || n == 4)
 205 	  h1->Fill(efficiency,purity);
 206 	if ( n == 1 || n == 3 || n == 5)
 207 	  h2->Fill(efficiency,purity);
 208 
 209 	if ( n == 6 || n == 8 || n == 10)
 210 	  h3->Fill(efficiency,purity);
 211 	if ( n == 7 || n == 9 || n == 11)
 212 	  h4->Fill(efficiency,purity);
 213 	
 214 	if ( n == 12 || n == 14 || n == 16)
 215 	  h5->Fill(efficiency,purity);
 216 	if ( n == 13 || n == 15 || n == 17)
 217 	  h6->Fill(efficiency,purity);
 218       }
 219     }
 220   }
 221   
 222   // Calculate total efficiencies and make plot beautiful with markers
 223   
 224   TMarker* marker[3];
 225   for ( int n(0); n < 3; n++ ){         // for all 3 Global processes
 226     marker[n] = new TMarker(0.0,0.0,20);
 227     marker[n]->SetMarkerColor(n+1);marker[n]->SetMarkerSize(1);
 228     
 229     int m = 0;  // counter for every tenth cone angle, needed later.
 230     
 231     if ( n==0 ) sprintf(Process,"nna");
 232     if ( n==1 ) sprintf(Process,"nnaa");
 233     if ( n==2 ) sprintf(Process,"nnaaa");
 234     
 235     
 236     for ( int l(0); l < 100; l++){     // for all cone angles to evaluate
 237       double coneangle = GlobalMatchCounter[n][l][0];
 238       double FoundTrueMatches = GlobalMatchCounter[n][l][1];
 239       double TrueMatches = GlobalMatchCounter[n][l][2];
 240       double FoundMatches = GlobalMatchCounter[n][l][3];
 241       
 242       if (TrueMatches != 0 && FoundMatches != 0){
 243   	double efficiency = FoundTrueMatches/TrueMatches;
 244   	double purity = FoundTrueMatches/FoundMatches;
 245   	
 246   	h7->Fill(efficiency,purity);
 247 	
 248   	// for every tenth opening angle draw marker in histo
 249 	
 250   	if( l%10 == 0 && coneangle > 0.03 && coneangle < 0.16 ){
 251   	  c7->cd();
 252   	  marker[n]->SetX(efficiency);marker[n]->SetY(purity);
 253   	  marker[n]->SetMarkerStyle(20+m);
 254   	  marker[n]->DrawClone();
 255 
 256   	  if (n == 2){
 257 	  cout <<"asd\n";
 258   	    TLatex latex;
 259   	    latex.SetTextSize(0.04);
 260   	    latex.SetTextAlign(13);
 261   	    latex.DrawLatex(efficiency-0.02,purity-0.002,Form("%.2f",coneangle));
 262   	  }
 263   	  m++;
 264   	}
 265       }
 266       //coneangle += 0.002; 
 267     }
 268     marker[n]->SetMarkerStyle(20);
 269   }
 270 
 271   //
 272   //  h1->GetXaxis()->SetLabelSize(0.04);
 273   //  h1->SetMarkerStyle(6);
 274   //  h2->GetXaxis()->SetLabelSize(0.04);
 275   //  h2->SetMarkerStyle(6);
 276   //  h3->GetXaxis()->SetLabelSize(0.04);
 277   //  h3->SetMarkerStyle(6);
 278   //  h4->GetXaxis()->SetLabelSize(0.04);
 279   //  h4->SetMarkerStyle(6);
 280   //  h5->GetXaxis()->SetLabelSize(0.04);
 281   //  h5->SetMarkerStyle(6);
 282   //  h6->GetXaxis()->SetLabelSize(0.04);
 283   //  h6->SetMarkerStyle(6);
 284   //  h7->GetXaxis()->SetLabelSize(0.04);
 285   //  h7->SetMarkerStyle(6);
 286 
 287 
 288 
 289 
 290   // Draw
 291 
 292   c1->cd();
 293   h1->DrawCopy();
 294 
 295   c2->cd();
 296   h2->DrawCopy();
 297   
 298   c3->cd();
 299   h3->DrawCopy();
 300   
 301   c4->cd();
 302   h4->DrawCopy();
 303   
 304   c5->cd();
 305   h5->DrawCopy();
 306   
 307   c6->cd();
 308   h6->DrawCopy();
 309   
 310   c7->cd();
 311   // h7->DrawCopy();
 312 
 313 
 314   // Draw legend 
 315 
 316    TLegend* leg =new TLegend(0.2,0.2,0.6,0.6);
 317    c7->cd();
 318    leg->AddEntry(marker[0],"#nu#nu#gamma","P");
 319    leg->AddEntry(marker[1],"#nu#nu#gamma#gamma","P");
 320    leg->AddEntry(marker[2],"#nu#nu#gamma#gamma#gamma","P");
 321    leg->Draw();
 322 
 323 
 324 
 325    //  // Print to file
 326    char plotname[256];
 327    
 328 
 329    sprintf(plotname,"../plots/MergeAngle_n1n1_Pe-1_Pp1_Emin%.1f_Emax%.1f_Cosmin%.3f_Cosmax%.1f_cut%i.eps",
 330 	   e_min, e_max, costheta_min, costheta_min, tracking_cut);
 331    c1->Print(plotname);
 332    sprintf(plotname,"../plots/MergeAngle_n1n1_Pe1_Pp-1_Emin%.1f_Emax%.1f_Cosmin%.3f_Cosmax%.1f_cut%i.eps",
 333 	   e_min, e_max, costheta_min, costheta_min, tracking_cut);   
 334    c2->Print(plotname);
 335    sprintf(plotname,"../plots/MergeAngle_n2n2_Pe-1_Pp1_Emin%.1f_Emax%.1f_Cosmin%.3f_Cosmax%.1f_cut%i.eps",
 336 	   e_min, e_max, costheta_min, costheta_min, tracking_cut);
 337    c3->Print(plotname);
 338    sprintf(plotname,"../plots/MergeAngle_n2n2_Pe1_Pp-1_Emin%.1f_Emax%.1f_Cosmin%.3f_Cosmax%.1f_cut%i.eps",
 339 	   e_min, e_max, costheta_min, costheta_min, tracking_cut);
 340    c4->Print(plotname);
 341    sprintf(plotname,"../plots/MergeAngle_n3n3_Pe-1_Pp1_Emin%.1f_Emax%.1f_Cosmin%.3f_Cosmax%.1f_cut%i.eps",
 342 	   e_min, e_max, costheta_min, costheta_min, tracking_cut);
 343    c5->Print(plotname);
 344    sprintf(plotname,"../plots/MergeAngle_n3n3_Pe1_Pp-1_Emin%.1f_Emax%.1f_Cosmin%.3f_Cosmax%.1f_cut%i.eps",
 345 	   e_min, e_max, costheta_min, costheta_min, tracking_cut);
 346    c6->Print(plotname);
 347    sprintf(plotname,"../plots/MergeAngle_Total_Emin%.1f_Emax%.1f_Cosmin%.3f_Cosmax%.1f_cut%i.eps",
 348 	   e_min, e_max, costheta_min, costheta_min, tracking_cut);
 349    c7->Print(plotname);
 350   
 351   return;
 352 }

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.