Safe dynamic casting from root dictionary fascility in GetGenerator and GetAliGenHBTp...
[u/mrichter/AliRoot.git] / HBTAN / AliHBTAnalysis.cxx
CommitLineData
1b446896 1
2#include "AliHBTAnalysis.h"
3
4#include <iostream.h>
5
6#include "AliHBTRun.h"
7#include "AliHBTReader.h"
8#include "AliHBTParticle.h"
9#include "AliHBTParticleCut.h"
10#include "AliHBTPair.h"
11#include "AliHBTPairCut.h"
12#include "AliHBTFunction.h"
13
14#include <TList.h>
15
16
17
18ClassImp(AliHBTAnalysis)
19
20const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100;
21const Int_t AliHBTAnalysis::fgkHbtAnalyzeAll = 0;
22
23AliHBTAnalysis::AliHBTAnalysis()
24 {
25 fReader = 0x0;
26
27 fTrackFunctions = new AliHBTTwoPartFctn* [fgkFctnArraySize];
28 fParticleFunctions = new AliHBTTwoPartFctn* [fgkFctnArraySize];
29 fParticleAndTrackFunctions = new AliHBTFourPartFctn* [fgkFctnArraySize];
30
31 fNTrackFunctions = 0;
32 fNParticleFunctions = 0;
33 fNParticleAndTrackFunctions = 0;
34
35 fPairCut = new AliHBTEmptyPairCut();//empty cut - accepts all particles
36
37 }
38
39AliHBTAnalysis::~AliHBTAnalysis()
40 {
41 //destructor
42 //note that we do not delete functions itself
43 // they should be deleted by whom where created
44 //we only store pointers, and we use "new" only for pointers array
45 delete [] fTrackFunctions;
46 delete [] fParticleFunctions;
47 delete [] fParticleAndTrackFunctions;
48
49 delete fPairCut; // always have an copy of an object - we create we dstroy
50 }
51
52/*************************************************************************************/
53void AliHBTAnalysis::Process(Option_t* option)
54{
55 //default option = "TracksAndParticles"
56 //Main method of the HBT Analysis Package
57 //It triggers reading with the global cut (default is an empty cut)
58 //Than it checks options and data which are read
59 //if everything is OK, then it calls one of the looping methods
60 //depending on tfReaderhe option
61 //These methods differs on what thay are looping on
62 //
63 // METHOD OPTION
64 //--------------------------------------------------------------------
65 //ProcessTracksAndParticles - "TracksAndParticles"
66 // DEFAULT
67 // it loops over both, tracks(reconstructed) and particles(simulated)
68 // all function gethered in all 3 lists are called for each (double)pair
69 //
70 //ProcessTracks - "Tracks"
71 // it loops only on tracks(reconstructed),
72 // functions ONLY from fTrackFunctions list are called
73 //
74 //ProcessParticles - "Particles"
75 // it loops only on particles(simulated),
76 // functions ONLY from fParticleAndTrackFunctions list are called
77 //
78 //
79 if (!fReader)
80 {
81 Error("AliHBTAnalysis::Process","The reader is not set");
82 return;
83 }
84
85
86 const char *oT = strstr(option,"Tracks");
87 const char *oP = strstr(option,"Particles");
88
89 if(oT && oP)
90 {
91 if (fReader->GetNumberOfPartEvents() <1)
92 {
93 Error("AliHBTAnalysis::Process","There is no Particles. Maybe change the option?");
94 return;
95 }
96 if (fReader->GetNumberOfTrackEvents() <1)
97 {
98 Error("AliHBTAnalysis::Process","There is no Tracks. Maybe change the option?");
99 return;
100 }
101
102 if ( RunCoherencyCheck() )
103 {
104 Error("AliHBTAnalysis::Process",
105 "Coherency check not passed. Maybe change the option?\n");
106 return;
107 }
108 ProcessTracksAndParticles();
109 return;
110 }
111
112 if(oT)
113 {
114 ProcessTracks();
115 return;
116 }
117
118 if(oP)
119 {
120 ProcessParticles();
121 return;
122 }
123
124}
125
126/*************************************************************************************/
127void AliHBTAnalysis::ProcessTracksAndParticles()
128{
129
130//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
131//the loops are splited
132
133
134 AliHBTParticle * part1, * part2;
135 AliHBTParticle * track1, * track2;
136
137 AliHBTEvent * trackEvent, *partEvent;
138 AliHBTEvent * trackEvent2,*partEvent2;
139
140// Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
141
142 Int_t Nev = fReader->GetNumberOfTrackEvents();
143
144 /***************************************/
145 /****** Looping same events ********/
146 /****** filling numerators ********/
147 /***************************************/
148 AliHBTPair * trackpair = new AliHBTPair();
149 AliHBTPair * partpair = new AliHBTPair();
150
151 for (Int_t i = 0;i<Nev;i++)
152 {
153 partEvent= fReader->GetParticleEvent(i);
154 trackEvent = fReader->GetTrackEvent(i);
155
156 if (!partEvent) continue;
157
158 //N = 0;
159
160 for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
161 {
162 if ( (j%10) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
163
164 part1= partEvent->GetParticle(j);
165 track1= trackEvent->GetParticle(j);
166
167 for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
168 {
169 part2= partEvent->GetParticle(k);
170 partpair->SetParticles(part1,part2);
171
172 track2= trackEvent->GetParticle(k);
173 trackpair->SetParticles(track1,track2);
174
175 if(fPairCut->Pass(partpair) || (fPairCut->Pass(trackpair))) //check pair cut
176 { //do not meets crietria of the pair cut, try with swaped pairs
177 if( ( fPairCut->Pass(partpair->GetSwapedPair()) ) || ( fPairCut->Pass(trackpair->GetSwapedPair()) ) )
178 continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
179 else
180 { //swaped pair meets all the criteria
181 partpair = partpair->GetSwapedPair();
182 trackpair = trackpair->GetSwapedPair();
183 }
184 }
185 UInt_t ii;
186 for(ii = 0;ii<fNParticleFunctions;ii++)
187 fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
188
189 for(ii = 0;ii<fNTrackFunctions;ii++)
190 fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
191
192 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
193 fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(trackpair,partpair);
194 }
195 }
196 }
197
198 /***************************************/
199 /***** Filling denominators *********/
200 /***************************************/
201 for (Int_t i = 0;i<Nev;i++) //In each event ....
202 {
203
204 partEvent= fReader->GetParticleEvent(i);
205 if (!partEvent) continue;
206
207 trackEvent = fReader->GetTrackEvent(i);
208
209// N=0;
210
211 for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
212 {
213// if (N>MAXCOMB) break;
214
215 part1= partEvent->GetParticle(j);
216
217 track1= trackEvent->GetParticle(j);
218
219// for (Int_t k = i+1; k<Nev;k++) // ... Loop over all proceeding events ...
220 Int_t NNN;
221
222 if ( (i+2) < Nev) NNN = i+2;
223 else NNN = Nev;
224
225 for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
226 {
227
228 partEvent2= fReader->GetParticleEvent(k);
229 if (!partEvent2) continue;
230
231 trackEvent2 = fReader->GetTrackEvent(k);
232
233 if ( (j%10) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
234
235 for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
236 {
237
238 // if (N>MAXCOMB) break;
239
240 part2= partEvent2->GetParticle(l);
241 partpair->SetParticles(part1,part2);
242
243 track2= trackEvent2->GetParticle(l);
244 trackpair->SetParticles(track1,track2);
245
246 if( (fPairCut->Pass(partpair)) || (fPairCut->Pass(trackpair)) ) //check pair cut
247 { //do not meets crietria of the
248 if( ( fPairCut->Pass(partpair->GetSwapedPair()) ) || ( fPairCut->Pass(trackpair->GetSwapedPair()) ) )
249 continue;
250 else
251 {
252 partpair = partpair->GetSwapedPair();
253 trackpair = trackpair->GetSwapedPair();
254 }
255 }
256 UInt_t ii;
257 for(ii = 0;ii<fNParticleFunctions;ii++)
258 fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
259
260 for(ii = 0;ii<fNTrackFunctions;ii++)
261 fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
262
263 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
264 fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(trackpair,partpair);
265
266
267 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
268 }//for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
269 }
270
271 }
272
273 /***************************************/
274
275
276}
277/*************************************************************************************/
278
279void AliHBTAnalysis::ProcessTracks()
280{
281 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
282//the loops are splited
283 AliHBTParticle * track1, * track2;
284
285 AliHBTEvent * trackEvent;
286 AliHBTEvent * trackEvent2;
287
288// Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
289
290 Int_t Nev = fReader->GetNumberOfTrackEvents();
291
292 /***************************************/
293 /****** Looping same events ********/
294 /****** filling numerators ********/
295 /***************************************/
296 AliHBTPair * trackpair = new AliHBTPair();
297
298 for (Int_t i = 0;i<Nev;i++)
299 {
300 trackEvent = fReader->GetTrackEvent(i);
301 if (!trackEvent) continue;
302 //N = 0;
303
304 for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
305 {
306 if ( (j%10) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
307
308 track1= trackEvent->GetParticle(j);
309
310 for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
311 {
312 track2= trackEvent->GetParticle(k);
313 trackpair->SetParticles(track1,track2);
314 if(fPairCut->Pass(trackpair)) //check pair cut
315 { //do not meets crietria of the
316 if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
317 else trackpair = trackpair->GetSwapedPair();
318 }
319
320 UInt_t ii;
321
322 for(ii = 0;ii<fNTrackFunctions;ii++)
323 fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
324
325
326 }
327 }
328 }
329
330 /***************************************/
331 /***** Filling diff histogram *********/
332 /***************************************/
333 for (Int_t i = 0;i<Nev;i++) //In each event ....
334 {
335 trackEvent = fReader->GetTrackEvent(i);
336 if (!trackEvent) continue;
337// N=0;
338
339 for (Int_t j = 0; j< trackEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
340 {
341// if (N>MAXCOMB) break;
342
343 track1= trackEvent->GetParticle(j);
344
345// for (Int_t k = i+1; k<Nev;k++) // ... Loop over all proceeding events ...
346 Int_t NNN;
347
348 if ( (i+2) < Nev) NNN = i+2;
349 else NNN = Nev;
350
351 for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
352 {
353
354 trackEvent2 = fReader->GetTrackEvent(k);
355 if (!trackEvent2) continue;
356
357 if ( (j%10) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
358
359 for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++) // ... on all particles
360 {
361
362 // if (N>MAXCOMB) break;
363 track2= trackEvent2->GetParticle(l);
364 trackpair->SetParticles(track1,track2);
365
366 if(fPairCut->Pass(trackpair)) //check pair cut
367 { //do not meets crietria of the
368 if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
369 else trackpair = trackpair->GetSwapedPair();
370 }
371 UInt_t ii;
372 for(ii = 0;ii<fNTrackFunctions;ii++)
373 fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
374
375 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
376 }//for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
377 }
378
379 }
380
381 /***************************************/
382
383
384}
385
386/*************************************************************************************/
387void AliHBTAnalysis::ProcessParticles()
388{
389 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
390//the loops are splited
391 AliHBTParticle * part1, * part2;
392
393 AliHBTEvent * partEvent;
394 AliHBTEvent * partEvent2;
395
396// Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
397
398 Int_t Nev = fReader->GetNumberOfPartEvents();
399
400 /***************************************/
401 /****** Looping same events ********/
402 /****** filling numerators ********/
403 /***************************************/
404 AliHBTPair * partpair = new AliHBTPair();
405
406 for (Int_t i = 0;i<Nev;i++)
407 {
408 partEvent= fReader->GetParticleEvent(i);
409 if (!partEvent) continue;
410 //N = 0;
411
412 for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
413 {
414 if ( (j%10) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
415
416 part1= partEvent->GetParticle(j);
417
418 for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
419 {
420 part2= partEvent->GetParticle(k);
421 partpair->SetParticles(part1,part2);
422
423 if( fPairCut->Pass(partpair) ) //check pair cut
424 { //do not meets crietria of the pair cut, try with swaped pairs
425 if( fPairCut->Pass(partpair->GetSwapedPair() ) )
426 continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
427 else
428 { //swaped pair meets all the criteria
429 partpair = partpair->GetSwapedPair();
430 }
431 }
432
433 UInt_t ii;
434
435 for(ii = 0;ii<fNParticleFunctions;ii++)
436 fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
437 }
438 }
439 }
440
441 /***************************************/
442 /***** Filling diff histogram *********/
443 /***************************************/
444 for (Int_t i = 0;i<Nev;i++) //In each event ....
445 {
446 partEvent= fReader->GetParticleEvent(i);
447 if (!partEvent) continue;
448
449// N=0;
450
451 for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
452 {
453// if (N>MAXCOMB) break;
454
455 part1= partEvent->GetParticle(j);
456
457// for (Int_t k = i+1; k<Nev;k++) // ... Loop over all proceeding events ...
458 Int_t NNN;
459
460 if ( (i+2) < Nev) NNN = i+2; //loop over next event
461 else NNN = Nev;
462
463 for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
464 {
465
466 partEvent2= fReader->GetParticleEvent(k);
467 if (!partEvent2) continue;
468
469 if ( (j%10) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
470
471 for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
472 {
473
474 // if (N>MAXCOMB) break;
475 part2= partEvent2->GetParticle(l);
476 partpair->SetParticles(part1,part2);
477
478 if(fPairCut->Pass(partpair)) //check pair cut
479 { //do not meets crietria of the
480 if( fPairCut->Pass(partpair->GetSwapedPair()) ) continue;
481 else partpair = partpair->GetSwapedPair();
482 }
483 UInt_t ii;
484 for(ii = 0;ii<fNParticleFunctions;ii++)
485 fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
486
487 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
488 }//for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
489 }
490
491 }
492
493 /***************************************/
494
495
496}
497
498/*************************************************************************************/
499
500
501void AliHBTAnalysis::Write()
502{
503 UInt_t ii;
504 for(ii = 0;ii<fNParticleFunctions;ii++)
505 fParticleFunctions[ii]->Write();
506
507 for(ii = 0;ii<fNTrackFunctions;ii++)
508 fTrackFunctions[ii]->Write();
509
510 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
511 fParticleAndTrackFunctions[ii]->Write();
512}
513/*************************************************************************************/
514void AliHBTAnalysis::SetGlobalPairCut(AliHBTPairCut* cut)
515{
516 if (cut == 0x0)
517 {
518 Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
519 }
520 delete fPairCut;
521 fPairCut = (AliHBTPairCut*)cut->Clone();
522}
523
524/*************************************************************************************/
525void AliHBTAnalysis::AddTrackFunction(AliHBTTwoPartFctn* f)
526{
527 if (f == 0x0) return;
528 if (fNTrackFunctions == fgkFctnArraySize)
529 {
530 Error("AliHBTAnalysis::AddTrackFunction","Can not add this function, not enough place in the array.");
531 }
532 fTrackFunctions[fNTrackFunctions] = f;
533 fNTrackFunctions++;
534}
535void AliHBTAnalysis::AddParticleFunction(AliHBTTwoPartFctn* f)
536{
537 if (f == 0x0) return;
538
539 if (fNParticleFunctions == fgkFctnArraySize)
540 {
541 Error("AliHBTAnalysis::AddParticleFunction","Can not add this function, not enough place in the array.");
542 }
543 fParticleFunctions[fNParticleFunctions] = f;
544 fNParticleFunctions++;
545
546
547}
548void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTFourPartFctn* f)
549{
550 if (f == 0x0) return;
551 if (fNParticleAndTrackFunctions == fgkFctnArraySize)
552 {
553 Error("AliHBTAnalysis::AddParticleAndTrackFunction","Can not add this function, not enough place in the array.");
554 }
555 fParticleAndTrackFunctions[fNParticleAndTrackFunctions] = f;
556 fNParticleAndTrackFunctions++;
557}
558
559
560/*************************************************************************************/
561
562
563/*************************************************************************************/
564Bool_t AliHBTAnalysis::RunCoherencyCheck()
565{
566 //Checks if both HBTRuns are similar
567 //return true if error found
568 //if they seem to be OK return false
569 Int_t i;
570 cout<<"Checking HBT Runs Coherency"<<endl;
571
572 cout<<"Number of events ..."; fflush(stdout);
573
574 if (fReader->GetNumberOfPartEvents() == fReader->GetNumberOfTrackEvents() ) //check whether there is the same number of events
575 {
576 cout<<"OK. "<<fReader->GetNumberOfTrackEvents()<<" found"<<endl;
577 }
578 else
579 { //if not the same - ERROR
580 Error("AliHBTAnalysis::RunCoherencyCheck()",
581 "Number of simulated events (%d) is not equal to number of reconstructed events(%d)",
582 fReader->GetNumberOfPartEvents(),fReader->GetNumberOfTrackEvents());
583 return kTRUE;
584 }
585
586 cout<<"Checking number of Particles AND Particles Types in each event ...";fflush(stdout);
587
588 AliHBTEvent *partEvent;
589 AliHBTEvent *trackEvent;
590 for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
591 {
592 partEvent= fReader->GetParticleEvent(i); //gets the "ith" event
593 trackEvent = fReader->GetTrackEvent(i);
594
595 if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
596 if ( (partEvent == 0x0) || (partEvent == 0x0) )
597 {
598 Error("AliHBTAnalysis::RunCoherencyCheck()",
599 "One event is NULL and the other one not. Event Number %d",i);
600 return kTRUE;
601 }
602
603 if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
604 {
605 Error("AliHBTAnalysis::RunCoherencyCheck()",
606 "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
607 i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
608 return kTRUE;
609 }
610 else
611 for (Int_t j = 0; j<partEvent->GetNumberOfParticles(); j++)
612 {
613 if( partEvent->GetParticle(j)->GetPdgCode() != trackEvent->GetParticle(j)->GetPdgCode() )
614 {
615 Error("AliHBTAnalysis::RunCoherencyCheck()",
616 "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
617 i,j, partEvent->GetParticle(j)->GetPdgCode(),trackEvent->GetParticle(j)->GetPdgCode() );
618 return kTRUE;
619
620 }
621 }
622 }
623 cout<<" Done"<<endl;
624 cout<<" Everything looks OK"<<endl;
625 return kFALSE;
626}
627
628
629/*************************************************************************************/
630
631
632/*************************************************************************************/
633
634