source: trunk/MagicSoft/Mars/mtemp/mifae/macros/mergeClean.C@ 5999

Last change on this file since 5999 was 4822, checked in by rico, 21 years ago
*** empty log message ***
File size: 14.8 KB
Line 
1/* ======================================================================== *\
2! Take different MC input files and combine them into single one
3! Author(s): Oscar Blanch & Javier Rico
4\* ======================================================================== */
5
6void mergeClean()
7{
8 const TString fileout = "/local_disk/jrico/mc/prueba.root";
9
10 const Int_t nfiles = 11;
11 const Int_t nlevels= 5;
12
13 const Char_t* levelname[nlevels]={"noclean","3015","4025","5035","6045"};
14
15 // input files -- no noise
16 const Char_t* NonoiseOutFilename[nfiles] =
17 {"/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0_cleaned_nonoise.root",
18 "/local_disk/jrico/mc/Gamma_zbin0_0_7_1130to1139_w0_cleaned_nonoise.root",
19 "/local_disk/jrico/mc/Gamma_zbin0_90_7_1260to1269_w0_cleaned_nonoise.root",
20 "/local_disk/jrico/mc/Gamma_zbin1_0_7_1010to1019_w0_cleaned_nonoise.root",
21 "/local_disk/jrico/mc/Gamma_zbin1_0_7_1140to1149_w0_cleaned_nonoise.root",
22 "/local_disk/jrico/mc/Gamma_zbin1_90_7_1270to1279_w0_cleaned_nonoise.root",
23 "/local_disk/jrico/mc/Gamma_zbin2_0_7_1020to1029_w0_cleaned_nonoise.root",
24 "/local_disk/jrico/mc/Gamma_zbin2_0_7_1150to1159_w0_cleaned_nonoise.root",
25 "/local_disk/jrico/mc/Gamma_zbin3_0_7_1030to1039_w0_cleaned_nonoise.root",
26 "/local_disk/jrico/mc/Gamma_zbin3_0_7_1160to1169_w0_cleaned_nonoise.root",
27 "/local_disk/jrico/mc/Gamma_zbin3_90_7_1290to1299_w0_cleaned_nonoise.root"
28 };
29
30 // input files -- different cleaning levels
31 const Char_t* NoiseOutFilename[nlevels][nfiles] =
32 {{"/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0_cleaned_noclean.root",
33 "/local_disk/jrico/mc/Gamma_zbin0_0_7_1130to1139_w0_cleaned_noclean.root",
34 "/local_disk/jrico/mc/Gamma_zbin0_90_7_1260to1269_w0_cleaned_noclean.root",
35 "/local_disk/jrico/mc/Gamma_zbin1_0_7_1010to1019_w0_cleaned_noclean.root",
36 "/local_disk/jrico/mc/Gamma_zbin1_0_7_1140to1149_w0_cleaned_noclean.root",
37 "/local_disk/jrico/mc/Gamma_zbin1_90_7_1270to1279_w0_cleaned_noclean.root",
38 "/local_disk/jrico/mc/Gamma_zbin2_0_7_1020to1029_w0_cleaned_noclean.root",
39 "/local_disk/jrico/mc/Gamma_zbin2_0_7_1150to1159_w0_cleaned_noclean.root",
40 "/local_disk/jrico/mc/Gamma_zbin3_0_7_1030to1039_w0_cleaned_noclean.root",
41 "/local_disk/jrico/mc/Gamma_zbin3_0_7_1160to1169_w0_cleaned_noclean.root",
42 "/local_disk/jrico/mc/Gamma_zbin3_90_7_1290to1299_w0_cleaned_noclean.root"},
43 {"/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0_cleaned_3015.root",
44 "/local_disk/jrico/mc/Gamma_zbin0_0_7_1130to1139_w0_cleaned_3015.root",
45 "/local_disk/jrico/mc/Gamma_zbin0_90_7_1260to1269_w0_cleaned_3015.root",
46 "/local_disk/jrico/mc/Gamma_zbin1_0_7_1010to1019_w0_cleaned_3015.root",
47 "/local_disk/jrico/mc/Gamma_zbin1_0_7_1140to1149_w0_cleaned_3015.root",
48 "/local_disk/jrico/mc/Gamma_zbin1_90_7_1270to1279_w0_cleaned_3015.root",
49 "/local_disk/jrico/mc/Gamma_zbin2_0_7_1020to1029_w0_cleaned_3015.root",
50 "/local_disk/jrico/mc/Gamma_zbin2_0_7_1150to1159_w0_cleaned_3015.root",
51 "/local_disk/jrico/mc/Gamma_zbin3_0_7_1030to1039_w0_cleaned_3015.root",
52 "/local_disk/jrico/mc/Gamma_zbin3_0_7_1160to1169_w0_cleaned_3015.root",
53 "/local_disk/jrico/mc/Gamma_zbin3_90_7_1290to1299_w0_cleaned_3015.root"},
54 {"/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0_cleaned_4025.root",
55 "/local_disk/jrico/mc/Gamma_zbin0_0_7_1130to1139_w0_cleaned_4025.root",
56 "/local_disk/jrico/mc/Gamma_zbin0_90_7_1260to1269_w0_cleaned_4025.root",
57 "/local_disk/jrico/mc/Gamma_zbin1_0_7_1010to1019_w0_cleaned_4025.root",
58 "/local_disk/jrico/mc/Gamma_zbin1_0_7_1140to1149_w0_cleaned_4025.root",
59 "/local_disk/jrico/mc/Gamma_zbin1_90_7_1270to1279_w0_cleaned_4025.root",
60 "/local_disk/jrico/mc/Gamma_zbin2_0_7_1020to1029_w0_cleaned_4025.root",
61 "/local_disk/jrico/mc/Gamma_zbin2_0_7_1150to1159_w0_cleaned_4025.root",
62 "/local_disk/jrico/mc/Gamma_zbin3_0_7_1030to1039_w0_cleaned_4025.root",
63 "/local_disk/jrico/mc/Gamma_zbin3_0_7_1160to1169_w0_cleaned_4025.root",
64 "/local_disk/jrico/mc/Gamma_zbin3_90_7_1290to1299_w0_cleaned_4025.root"},
65 {"/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0_cleaned_5035.root",
66 "/local_disk/jrico/mc/Gamma_zbin0_0_7_1130to1139_w0_cleaned_5035.root",
67 "/local_disk/jrico/mc/Gamma_zbin0_90_7_1260to1269_w0_cleaned_5035.root",
68 "/local_disk/jrico/mc/Gamma_zbin1_0_7_1010to1019_w0_cleaned_5035.root",
69 "/local_disk/jrico/mc/Gamma_zbin1_0_7_1140to1149_w0_cleaned_5035.root",
70 "/local_disk/jrico/mc/Gamma_zbin1_90_7_1270to1279_w0_cleaned_5035.root",
71 "/local_disk/jrico/mc/Gamma_zbin2_0_7_1020to1029_w0_cleaned_5035.root",
72 "/local_disk/jrico/mc/Gamma_zbin2_0_7_1150to1159_w0_cleaned_5035.root",
73 "/local_disk/jrico/mc/Gamma_zbin3_0_7_1030to1039_w0_cleaned_5035.root",
74 "/local_disk/jrico/mc/Gamma_zbin3_0_7_1160to1169_w0_cleaned_5035.root",
75 "/local_disk/jrico/mc/Gamma_zbin3_90_7_1290to1299_w0_cleaned_5035.root"},
76 {"/local_disk/jrico/mc/Gamma_zbin0_0_7_1000to1009_w0_cleaned_6045.root",
77 "/local_disk/jrico/mc/Gamma_zbin0_0_7_1130to1139_w0_cleaned_6045.root",
78 "/local_disk/jrico/mc/Gamma_zbin0_90_7_1260to1269_w0_cleaned_6045.root",
79 "/local_disk/jrico/mc/Gamma_zbin1_0_7_1010to1019_w0_cleaned_6045.root",
80 "/local_disk/jrico/mc/Gamma_zbin1_0_7_1140to1149_w0_cleaned_6045.root",
81 "/local_disk/jrico/mc/Gamma_zbin1_90_7_1270to1279_w0_cleaned_6045.root",
82 "/local_disk/jrico/mc/Gamma_zbin2_0_7_1020to1029_w0_cleaned_6045.root",
83 "/local_disk/jrico/mc/Gamma_zbin2_0_7_1150to1159_w0_cleaned_6045.root",
84 "/local_disk/jrico/mc/Gamma_zbin3_0_7_1030to1039_w0_cleaned_6045.root",
85 "/local_disk/jrico/mc/Gamma_zbin3_0_7_1160to1169_w0_cleaned_6045.root",
86 "/local_disk/jrico/mc/Gamma_zbin3_90_7_1290to1299_w0_cleaned_6045.root"}
87 };
88
89
90 // Getting input braches
91 TChain* tno = new TChain("Events");
92 TChain* tyes[nlevels];
93 for(Int_t i=0;i<nlevels;i++)
94 tyes[i] = new TChain("Events");
95
96 for(Int_t i=0;i<nfiles;i++)
97 {
98 tno->Add(NonoiseOutFilename[i]);
99 for(Int_t j=0;j<nlevels;j++)
100 tyes[j]->Add(NoiseOutFilename[j][i]);
101 }
102
103 // input braches for no noise file
104 MCerPhotEvt* nphotNo;
105 MHillas* hillasNo;
106 MHillasSrc* hillassrcNo;
107 MMcEvt* mcevtNo;
108 MRawEvtHeader* evtheaderNo;
109 MRawRunHeader* runheaderNo;
110
111 tno->SetBranchAddress("MCerPhotEvt.",&nphotNo);
112 tno->SetBranchAddress("MHillas.",&hillasNo);
113 tno->SetBranchAddress("MHillasSrc.",&hillassrcNo);
114 tno->SetBranchAddress("MMcEvt.",&mcevtNo);
115 tno->SetBranchAddress("MRawEvtHeader.",&evtheaderNo);
116 tno->SetBranchAddress("MRawRunHeader.",&runheaderNo);
117
118
119 // input branches for noise files
120 MCerPhotEvt* nphotYes[nlevels];
121 MHillas* hillasYes[nlevels];
122 MHillasSrc* hillassrcYes[nlevels];
123 MMcEvt* mcevtYes[nlevels];
124 MRawEvtHeader* evtheaderYes[nlevels];
125 MRawRunHeader* runheaderYes[nlevels];
126
127 for(Int_t i=0;i<nlevels;i++)
128 {
129 tyes[i]->SetBranchAddress("MCerPhotEvt.",&nphotYes[i]);
130 tyes[i]->SetBranchAddress("MHillas.",&hillasYes[i]);
131 tyes[i]->SetBranchAddress("MHillasSrc.",&hillassrcYes[i]);
132 tyes[i]->SetBranchAddress("MMcEvt.",&mcevtYes[i]);
133 tyes[i]->SetBranchAddress("MRawEvtHeader.",&evtheaderYes[i]);
134 tyes[i]->SetBranchAddress("MRawRunHeader.",&runheaderYes[i]);
135 }
136
137 // Setting ouput branches
138 TFile* fout = new TFile(fileout,"RECREATE");
139 TTree* to = new TTree("Events","Events");
140
141 // common containers
142 TBranch* bmcevt = to->Branch("MMcEvt.", "MMcEvt", &mcevtNo);
143 TBranch* brawevtheader = to->Branch("MRawEvtHeader.","MRawEvtHeader",&evtheaderNo);
144 TBranch* brawrunheader = to->Branch("MRawRunHeader.","MRawRunHeader",&runheaderNo);
145
146 // noise/cleaning dependent containers
147 to->Branch("MCerPhotEvt0.", "MCerPhotEvt", &nphotNo);
148 to->Branch("MHillas0.", "MHillas", &hillasNo);
149 to->Branch("MHillasSrc0.", "MHillasSrc", &hillassrcNo);
150
151 // cout << "llego1" << endl;
152 for(Int_t i=0;i<nlevels;i++)
153 {
154 // cout << "i: " << i << endl;
155
156 char mcerphot[100];
157 char mhillas[100];
158 char mhillassrc[100];
159
160
161 sprintf(mcerphot, "MCerPhotEvt%s.",levelname[i]);
162 sprintf(mhillas, "MHillas%s." ,levelname[i]);
163 sprintf(mhillassrc,"MHillasSrc%s." ,levelname[i]);
164
165 to->Branch(mcerphot, "MCerPhotEvt", &nphotYes[i]);
166 to->Branch(mhillas, "MHillas", &hillasYes[i]);
167 to->Branch(mhillassrc,"MHillasSrc", &hillassrcYes[i]);
168 }
169 // cout << "llego2" << endl;
170
171
172 Int_t nentries1 = (Int_t)tno->GetEntries();
173 cout << "Entries: No noise file: " << nentries1 << endl;
174 Int_t nentries2[nlevels];
175 Int_t latestk[nlevels];
176 Int_t latestevt[nlevels];
177 Int_t nfound[nlevels];
178 for(Int_t i=0;i<nlevels;i++)
179 {
180 nfound[i]=0;
181 latestk[i]=0;
182 latestevt[i]=99999999;
183 nentries2[i] = (Int_t)tyes[i]->GetEntries();
184 cout << "Entries: cleaning level [" << levelname[i] << "] :" << nentries2[i] << endl;
185 }
186
187 // loop over all events in no-noise file and look for them in the other files
188 for(Int_t i=0;i<nentries1;i++)
189 // for(Int_t i=0;i<100;i++)
190 {
191 tno->GetEntry(i);
192 Int_t runno = runheaderNo->GetRunNumber();
193 Int_t evtno = evtheaderNo->GetDAQEvtNumber();
194
195 cout << "Checking run " << runno<< ", event " << evtno << endl;
196
197 //check if the no-noise event is present in all of the other files
198 for(Int_t j=0;j<nlevels;j++)
199 {
200 Int_t runyes;
201 Int_t evtyes;
202 cout << "Level " << levelname[j] << endl;
203 Int_t k=latestk[j];
204 while(k<nentries2[j])
205 {
206 // cout << k << endl;
207 tyes[j]->GetEntry(k++);
208 runyes = runheaderYes[j]->GetRunNumber();
209 evtyes = evtheaderYes[j]->GetDAQEvtNumber();
210 if(runyes==runno && evtyes==evtno)
211 break;
212 // next condition speeds up selection when events are grouped by runs and
213 // in event increasing order
214 if(runyes==runno && latestevt[j]<evtno && evtno<evtyes)
215 break;
216 // next condition speeds up selection when events are grouped by runs
217 // and in run increasing order
218 // if(runyes>runno)
219 // break;
220 }
221
222 if(k>=nentries2[j])
223 {
224 k=0;
225 while(k<latestk[j])
226 {
227 // cout << k << endl;
228 tyes[j]->GetEntry(k++);
229 runyes = runheaderYes[j]->GetRunNumber();
230 evtyes = evtheaderYes[j]->GetDAQEvtNumber();
231 if(runyes==runno && evtyes==evtno)
232 break;
233 // next condition speeds up selection when events are grouped by runs and
234 // in event increasing order
235 if(runyes==runno && latestevt[j]<evtno && evtno<evtyes)
236 break;
237 // next condition speeds up selection when events are grouped by runs
238 // and in run increasing order
239 // if(runyes>runno)
240 // break;
241 }
242 }
243
244 // the event has not been found, assign dummy containers
245 if(runyes=!runno || evtyes!=evtno)
246 {
247 cout << "Not Found!" << endl;
248 for(Int_t l=j;l<nlevels;l++)
249 {
250 nphotYes[l]->Reset();
251 hillasYes[l]->Reset();
252 hillassrcYes[l]->Reset();
253 }
254 break;
255 }
256 else
257 {
258 cout << "Found!" << endl;
259 latestevt[j] = evtyes;
260 latestk[j]=k;
261 nfound[j]++;
262 }
263 }
264 to->Fill();
265 }
266
267 // loop over all events in first noise file and look for them in the no-noise file in order to save also purely noisy events
268 cout << "***************************************" << endl;
269 cout << "*** Looking for purely noisy events ***" << endl;
270 cout << "***************************************" << endl;
271 Int_t nfound2[nlevels];
272 for(Int_t i=0;i<nlevels;i++)
273 {
274 nfound2[i]=0;
275 latestk[i]=0;
276 latestevt[i]=99999999;
277 }
278
279 // readdress the common container addresses
280 bmcevt->SetAddress(&mcevtYes[0]);
281 brawevtheader->SetAddress(&evtheaderYes[0]);
282 brawrunheader->SetAddress(&runheaderYes[0]);
283
284 for(Int_t i=0;i<nentries2[0];i++)
285 // for(Int_t i=0;i<100;i++)
286 {
287 tyes[0]->GetEntry(i);
288 Int_t run = runheaderYes[0]->GetRunNumber();
289 Int_t evt = evtheaderYes[0]->GetDAQEvtNumber();
290 Int_t runno;
291 Int_t evtno;
292
293 cout << "Checking run " << run << ", event " << evt << endl;
294
295 //check if the event is present in the no-noise file
296 Int_t kk=latestk[0];
297 while(kk<nentries1)
298 {
299 // cout << kk << endl;
300 tno->GetEntry(kk++);
301 runno = runheaderNo->GetRunNumber();
302 evtno = evtheaderNo->GetDAQEvtNumber();
303 if(run==runno && evt==evtno)
304 break;
305 if(run==runno && latestevt[0]<evt && evt<evtno)
306 break;
307 // if(runno>run)
308 // break;
309 }
310
311 if(kk>=nentries1)
312 {
313 kk=0;
314 while(kk<latestk[0])
315 {
316 // cout << kk << endl;
317 tno->GetEntry(kk++);
318 runno = runheaderNo->GetRunNumber();
319 evtno = evtheaderNo->GetDAQEvtNumber();
320 if(run==runno && evt==evtno)
321 break;
322 if(run==runno && latestevt[0]<evt && evt<evtno)
323 break;
324 // if(runno>run)
325 // break;
326 }
327 }
328
329 // the event is already included because it is also a no-noise event, dismiss it
330 if(run==runno && evt==evtno)
331 {
332 cout << "Found already in the no-noise sample, continuing... " << endl;
333 latestevt[0] = evtno;
334 latestk[0]=kk;
335 continue;
336 }
337
338 nfound2[0]++;
339 nphotNo->Reset();
340 hillasNo->Reset();
341 hillassrcNo->Reset();
342
343 // look for the event in the rest of the noisy samples
344 for(Int_t j=1;j<nlevels;j++)
345 {
346 Int_t runyes;
347 Int_t evtyes;
348 cout << "Level " << levelname[j] << endl;
349 Int_t k=latestk[j];
350 while(k<nentries2[j])
351 {
352 // cout << k << endl;
353 tyes[j]->GetEntry(k++);
354 runyes = runheaderYes[j]->GetRunNumber();
355 evtyes = evtheaderYes[j]->GetDAQEvtNumber();
356 if(runyes==run && evtyes==evt)
357 break;
358 if(runyes==run && latestevt[j]<evt && evt<evtyes)
359 break;
360 // if(runyes>run)
361 // break;
362 }
363
364 if(k>=nentries2[j])
365 {
366 k=0;
367 while(k<latestk[j])
368 {
369 // cout << k << endl;
370 tyes[j]->GetEntry(k++);
371 runyes = runheaderYes[j]->GetRunNumber();
372 evtyes = evtheaderYes[j]->GetDAQEvtNumber();
373 if(runyes==run && evtyes==evt)
374 break;
375 if(runyes==run && latestevt[j]<evt && evt<evtyes)
376 break;
377 // if(runyes>run)
378 // break;
379 }
380 }
381
382 // the event has not been found, assign dummy containers
383 if(runyes=!run || evtyes!=evt)
384 {
385 cout << "Not Found!" << endl;
386 for(Int_t l=j;l<nlevels;l++)
387 {
388 nphotYes[l]->Reset();
389 hillasYes[l]->Reset();
390 hillassrcYes[l]->Reset();
391 }
392 break;
393 }
394 else
395 {
396 cout << "Found!" << endl;
397 latestevt[j] = evtyes;
398 latestk[j]=k;
399 nfound2[j]++;
400 }
401 }
402 to->Fill();
403 }
404
405 fout->Write();
406 fout->Close();
407
408 for(Int_t i=0;i<nlevels;i++)
409 cout << "Found " << nfound[i] << "+" << nfound2[i] << " for cleaning level " << levelname[i] << endl;
410
411}
Note: See TracBrowser for help on using the repository browser.