Attachment 'PreprocessData.C'

Download

   1 // Preprocessing of SM background files.
   2 // Merges PFLow photon candidates depending on angle seen from IP between
   3 // PFlow photon candidates.
   4 // Adds new leaves for ProcessedRecoParticles to tree.
   5 //
   6 //
   7 // 
   8 // 1. MergeRecoPhoton:       
   9 //                              Input parameter: new file name. Routine scans RecoParticles for the
  10 //                              PFlow photon with the highest energy, and then adds up the energies of all
  11 //                              PFlow photons within the MergeConeAngle (of course the first photon is included)
  12 //                              Then it proceeds with the next highest photon energy that is not merged into 
  13 //                              the previous one. Finally all other particles (that are not photons)
  14 //                              are appended to ProcessedReco list.
  15 //                              The new photon momentum is the scaled momenta sum (sum):
  16 //                              p_x = p_{x,sum}*E_{new}/p_{sum}, thus p^2 = p_{sum}^2 * E_{new}^2/p_{sum}^2 = E_{new}^2
  17 //
  18 // 2. ClearMergedArrays:        Helper function to clear all processed reco arrays before every event
  19 //
  20 // 3. MergeRecoPhotonMethodB:   NOT YET IMPLEMENTED
  21 //                              Same as 1., however the new ProcessedReco momenta are calculated differently.
  22 //                              Instead of just scaling the momentum of the highest particle contribution to
  23 //                              the ProcessedReco photon, the new momentum is a weighted sum of the contributors i.
  24 //                              The momenta are either weighted with E or 1/sqrt{E}. 
  25 //                              In the first case:  p_x = sum{p_{x,i}*E_{i}/E_{new}}     
  26 //                              In the second case: p_x = sum{p_{x,i}*1/sqrt{/E_{i}}/sum{1/sqrt{E_{i}}}}
  27 //                              Set the method in routine with variable "method" !!
  28 
  29 
  30 # include "GlobalVariables.C"
  31 
  32 
  33 TFile* f(0);               // file
  34 TTree* tree(0);            // tree
  35 TTree* preproctree(0);     // preprocessed tree to be created
  36 
  37 
  38 double MergeConeAngle(0.04);   // set merge cone angle here !!
  39 
  40 
  41 
  42 void PreprocessData(double mergeangle){
  43   MergeConeAngle = mergeangle;
  44   MergeRecoPhotons();
  45   gROOT->ProcessLine(".q");
  46   return;
  47 }
  48 
  49 
  50 
  51 
  52 //#########################################################################
  53 //#########################################################################
  54 
  55 
  56 
  57 // 1. -------------------------------------------------------------------------
  58 
  59 int MergeRecoPhotons(){
  60 
  61   printf("Merge Cone Angle set to %.2f rad\n",MergeConeAngle);
  62 
  63   FILE* file = fopen("../steeringfiles/ILD00_Processed_RootFiles.txt","r"); 
  64   char filename[500];
  65 
  66   while(fscanf(file,"%s\n",filename) != EOF) {
  67 
  68     // __________ OPEN ROOTFILE, UPDATE PROCTREE __________________
  69 
  70     printf("Preprocessing file \"%s\"\n",filename);
  71 
  72     f = new TFile(filename,"UPDATE");
  73     tree = (TTree*)f->Get("Tree");
  74     LoadTree(tree);
  75 
  76     if (f->Get("ProcTree")){
  77       f->Get("ProcTree")->Delete("all");
  78     }
  79     preproctree = new TTree("ProcTree","ProcTree");
  80     SetPreProcTreeBranches(preproctree);
  81 
  82     // _______________ LOOP OVER EVENTS ______________________________-
  83 
  84     double angle(0.0);
  85 
  86     int nentries = tree->GetEntries();
  87 
  88     for(int i(0); i < nentries; i++){   // for each event
  89       if (i%10000 == 0) printf("Processing event: %i\n",i);
  90 
  91       tree->GetEntry(i);
  92       ClearMergedArrays();
  93 
  94       nProcessedRecos = 0;             // Number of Processed Recos in event, start with 0
  95 
  96       // Count all PFlow photons in event, and set their state to "not yet merged", RecoIsMerged == 0.
  97 
  98       int* RecoIsMerged = new int[nRecoParticles];
  99       int NumberOfPhotonsLeft = 0;
 100       for (int j(0); j < nRecoParticles; j++){
 101 	RecoIsMerged[j] = 0;
 102 	if (Reco_type[j] == 22){
 103 	  NumberOfPhotonsLeft++;
 104 	}
 105       }
 106 
 107       // while there are photons left which are not merged
 108 
 109       while (NumberOfPhotonsLeft > 0){     
 110 	int controlflag(0);             // set to one if a merged photon is created
 111 
 112 	double energysum = 0.0;
 113 	float weightsum = 0.0;
 114 	int fragmentcounter = 0;
 115 
 116 	double px = 0.0;
 117 	double py = 0.0;
 118 	double pz = 0.0;      
 119 
 120 	// find the most energetic not yet merged photon
 121 
 122 	double maxRecoEnergy = 0.0; 
 123 	int maxRecoEnergyId = 0;
 124 	for (int j(0); j < nRecoParticles; j++){  
 125 	  if ( Reco_type[j] == 22 && !RecoIsMerged[j] && Reco_energy[j] > maxRecoEnergy){
 126 	    maxRecoEnergy = Reco_energy[j];
 127 	    maxRecoEnergyId = j;
 128 	  }
 129 	}
 130 
 131 	// merge all not yet merged photons in cone to photon just found
 132 	// photon just found is of course merged with itself
 133 	int isBCal = 0;
 134 	for (int j(0); j < nRecoParticles; j++){   
 135 	  if (Reco_type[j] == 22 && !RecoIsMerged[j]){
 136 	    angle = TMath::ACos((Reco_px[j]*Reco_px[maxRecoEnergyId] + Reco_py[j]*Reco_py[maxRecoEnergyId] 
 137 				 + Reco_pz[j]*Reco_pz[maxRecoEnergyId])/(Reco_absP[j]*Reco_absP[maxRecoEnergyId]));
 138 	    if (angle < MergeConeAngle){
 139 	      energysum += Reco_energy[j];
 140 	      //weightsum += RelationWeight[j];
 141 	      weightsum = sqrt(weightsum * weightsum + RelationWeight[j] * RelationWeight[j]);
 142 	      px += Reco_px[j];
 143 	      py += Reco_py[j];
 144 	      pz += Reco_pz[j];
 145 
 146 	      if (Reco_isBCal[j]) isBCal = 1;
 147 
 148 	      RecoIsMerged[j] = 1;
 149 	      NumberOfPhotonsLeft--;
 150 	      fragmentcounter++;
 151 	      controlflag = 1;     // A merged photon has been created
 152 	    }
 153 	  }
 154 	}
 155  
 156 	// if a merged photon has been created fill ProcessedReco array
 157 
 158 	if (controlflag){
 159 
 160 	  ProcessedReco_energy[nProcessedRecos] = energysum;
 161 	  ProcessedReco_id[nProcessedRecos] = nProcessedRecos + 1;    // make sure that id starts from 1.
 162 	  ProcessedReco_type[nProcessedRecos] = 22;                   // new particle is photon
 163 	  ProcessedReco_isBCal[nProcessedRecos] = isBCal;
 164 
 165 	  ProcessedReco_px[nProcessedRecos] = px;
 166 	  ProcessedReco_py[nProcessedRecos] = py;
 167 	  ProcessedReco_pz[nProcessedRecos] = pz;
 168 
 169 	  ProcessedReco_mass[nProcessedRecos] = 0.0;
 170 	  ProcessedReco_charge[nProcessedRecos] = 0.0;
 171 	  ProcessedReco_position_x[nProcessedRecos] = Reco_position_x[maxRecoEnergyId];
 172 	  ProcessedReco_position_y[nProcessedRecos] = Reco_position_y[maxRecoEnergyId];
 173 	  ProcessedReco_position_z[nProcessedRecos] = Reco_position_z[maxRecoEnergyId];
 174 
 175 	  //normalise momentum
 176 
 177 	  double normalisation_factor = ProcessedReco_energy[nProcessedRecos]/
 178 	                    TMath::Sqrt(ProcessedReco_px[nProcessedRecos]*ProcessedReco_px[nProcessedRecos]+
 179 				        ProcessedReco_py[nProcessedRecos]*ProcessedReco_py[nProcessedRecos]+
 180 				        ProcessedReco_pz[nProcessedRecos]*ProcessedReco_pz[nProcessedRecos]);
 181 
 182 
 183 
 184 	  ProcessedReco_px[nProcessedRecos] = ProcessedReco_px[nProcessedRecos] * normalisation_factor;
 185 	  ProcessedReco_py[nProcessedRecos] = ProcessedReco_py[nProcessedRecos] * normalisation_factor;
 186 	  ProcessedReco_pz[nProcessedRecos] = ProcessedReco_pz[nProcessedRecos] * normalisation_factor;
 187 
 188 	  // calculate absP
 189 	  ProcessedReco_absP[nProcessedRecos] = TMath::Sqrt(ProcessedReco_px[nProcessedRecos]*ProcessedReco_px[nProcessedRecos]+
 190 							    ProcessedReco_py[nProcessedRecos]*ProcessedReco_py[nProcessedRecos]+
 191 							    ProcessedReco_pz[nProcessedRecos]*ProcessedReco_pz[nProcessedRecos]);
 192 
 193 	  ProcessedReco_cosTheta[nProcessedRecos] = ProcessedReco_pz[nProcessedRecos]/ProcessedReco_absP[nProcessedRecos];;
 194 	  ProcessedReco_phi[nProcessedRecos] = TMath::ATan2(ProcessedReco_py[nProcessedRecos],ProcessedReco_px[nProcessedRecos]);
 195 	
 196 	  ProcessedToRecoId[nProcessedRecos] = Reco_id[maxRecoEnergyId];
 197 	  ProcessedRecoRelatedToId[nProcessedRecos] = RecoToMCId[maxRecoEnergyId];
 198 	  ProcessedRecoRelationWeight[nProcessedRecos] = weightsum;
 199 	  NumberOfMerged[nProcessedRecos] = fragmentcounter;
 200 
 201 	  //printf("AbsP:  %.2f   ; Energy:  %2.f\n",ProcessedReco_absP[nProcessedRecos],ProcessedReco_energy[nProcessedRecos]);
 202 	  nProcessedRecos++;  // Add a new processed reco photon
 203 	}
 204       }
 205 
 206 
 207       for ( int k(0); k < nRecoParticles; k++){ //Add all particle left to new collection
 208 	if (RecoIsMerged[k] == 0){
 209 	  ProcessedReco_energy[nProcessedRecos] = Reco_energy[k];
 210 	  ProcessedReco_id[nProcessedRecos] = nProcessedRecos + 1;
 211 	  ProcessedReco_type[nProcessedRecos] = Reco_type [k];
 212 	  ProcessedReco_isBCal[nProcessedRecos] = Reco_isBCal [k];
 213 	  ProcessedReco_px[nProcessedRecos] = Reco_px[k];
 214 	  ProcessedReco_py[nProcessedRecos] = Reco_py[k];
 215 	  ProcessedReco_pz[nProcessedRecos] = Reco_pz[k];
 216 	  ProcessedReco_mass[nProcessedRecos] = Reco_mass[k];
 217 	  ProcessedReco_charge[nProcessedRecos] = Reco_charge[k];
 218 	  ProcessedReco_position_x[nProcessedRecos] = Reco_position_x[k];
 219 	  ProcessedReco_position_y[nProcessedRecos] = Reco_position_y [k];
 220 	  ProcessedReco_position_z[nProcessedRecos] = Reco_position_z [k];
 221 	  ProcessedReco_absP[nProcessedRecos] = Reco_absP[k];
 222 	  ProcessedReco_cosTheta[nProcessedRecos] =Reco_cosTheta[k];
 223 	  ProcessedReco_phi[nProcessedRecos] = Reco_phi[k];
 224 	
 225 	  ProcessedToRecoId[nProcessedRecos] = Reco_id[k];
 226 	  ProcessedRecoRelatedToId[nProcessedRecos] = RecoToMCId[k];
 227 	  ProcessedRecoRelationWeight[nProcessedRecos] = RelationWeight[k];
 228 	  NumberOfMerged[nProcessedRecos] = 1;
 229 	  nProcessedRecos++;     // Add a new processed reco particle     
 230 	}
 231       }
 232       delete[] RecoIsMerged;
 233       preproctree->Fill();
 234     }
 235     preproctree->Write();
 236     //f->Close();
 237     delete preproctree;
 238   }
 239   return 0;
 240 }
 241 
 242 
 243 int ClearMergedArrays(){
 244 
 245   for (int i(0); i < 1000; i++){
 246     ProcessedReco_id[i] = 0;
 247     ProcessedReco_type[i] = 0;
 248     ProcessedReco_isBCal[i] = 0;
 249     ProcessedReco_px[i] = 0.0;
 250     ProcessedReco_py[i] = 0.0;
 251     ProcessedReco_pz[i] = 0.0;
 252     ProcessedReco_energy[i] = 0.0;
 253     ProcessedReco_mass[i] = 0.0;
 254     ProcessedReco_charge[i] = 0.0;
 255     ProcessedReco_position_x[i] = 0.0;
 256     ProcessedReco_position_y[i] = 0.0;
 257     ProcessedReco_position_z[i] = 0.0;
 258     
 259     ProcessedReco_absP[i] = 0.0;
 260     ProcessedReco_cosTheta[i] = 0.0;
 261     ProcessedReco_phi[i] = 0.0;
 262     
 263     ProcessedToRecoId[i] = 0;
 264 
 265     ProcessedRecoRelatedToId[i] = 0;
 266     ProcessedRecoRelationWeight[i] = 0.0;
 267     NumberOfMerged[i] = 0;
 268   }
 269   return 0;
 270 }

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.