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.You are not allowed to attach a file to this page.