git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
example analysis macro extended for monitor functions
[u/mrichter/AliRoot.git] / HBTAN / AliHBTAnalysis.cxx
CommitLineData
1b446896 1
2#include "AliHBTAnalysis.h"
3
4#include <iostream.h>
5
6#include "AliHBTRun.h"
dc2c3f36 7#include "AliHBTEvent.h"
1b446896 8#include "AliHBTReader.h"
9#include "AliHBTParticle.h"
10#include "AliHBTParticleCut.h"
11#include "AliHBTPair.h"
12#include "AliHBTPairCut.h"
13#include "AliHBTFunction.h"
14
dc2c3f36 15#include <TBenchmark.h>
1b446896 16#include <TList.h>
17
18
19
20ClassImp(AliHBTAnalysis)
21
22const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100;
23const Int_t AliHBTAnalysis::fgkHbtAnalyzeAll = 0;
24
25AliHBTAnalysis::AliHBTAnalysis()
26 {
27 fReader = 0x0;
28
27b3fe5d 29 fTrackFunctions = new AliHBTOnePairFctn* [fgkFctnArraySize];
30 fParticleFunctions = new AliHBTOnePairFctn* [fgkFctnArraySize];
31 fParticleAndTrackFunctions = new AliHBTTwoPairFctn* [fgkFctnArraySize];
1b446896 32
33 fNTrackFunctions = 0;
34 fNParticleFunctions = 0;
35 fNParticleAndTrackFunctions = 0;
7a2c8238 36
1b446896 37 fPairCut = new AliHBTEmptyPairCut();//empty cut - accepts all particles
38
7a2c8238 39 fBufferSize = 2;
1b446896 40 }
491d1b5d 41/*************************************************************************************/
1b446896 42
43AliHBTAnalysis::~AliHBTAnalysis()
44 {
45 //destructor
46 //note that we do not delete functions itself
47 // they should be deleted by whom where created
48 //we only store pointers, and we use "new" only for pointers array
49 delete [] fTrackFunctions;
50 delete [] fParticleFunctions;
51 delete [] fParticleAndTrackFunctions;
52
53 delete fPairCut; // always have an copy of an object - we create we dstroy
54 }
55
56/*************************************************************************************/
491d1b5d 57
1b446896 58void AliHBTAnalysis::Process(Option_t* option)
59{
60 //default option = "TracksAndParticles"
61 //Main method of the HBT Analysis Package
62 //It triggers reading with the global cut (default is an empty cut)
63 //Than it checks options and data which are read
64 //if everything is OK, then it calls one of the looping methods
65 //depending on tfReaderhe option
66 //These methods differs on what thay are looping on
67 //
68 // METHOD OPTION
69 //--------------------------------------------------------------------
70 //ProcessTracksAndParticles - "TracksAndParticles"
71 // DEFAULT
72 // it loops over both, tracks(reconstructed) and particles(simulated)
73 // all function gethered in all 3 lists are called for each (double)pair
74 //
75 //ProcessTracks - "Tracks"
76 // it loops only on tracks(reconstructed),
77 // functions ONLY from fTrackFunctions list are called
78 //
79 //ProcessParticles - "Particles"
80 // it loops only on particles(simulated),
81 // functions ONLY from fParticleAndTrackFunctions list are called
82 //
83 //
84 if (!fReader)
85 {
01725374 86 Error("Process","The reader is not set");
1b446896 87 return;
88 }
89
90
91 const char *oT = strstr(option,"Tracks");
92 const char *oP = strstr(option,"Particles");
93
dc2c3f36 94 Bool_t nonid = IsNonIdentAnalysis();
95
1b446896 96 if(oT && oP)
97 {
98 if (fReader->GetNumberOfPartEvents() <1)
99 {
01725374 100 Error("Process","There is no Particles. Maybe change the option?");
1b446896 101 return;
102 }
103 if (fReader->GetNumberOfTrackEvents() <1)
104 {
01725374 105 Error("Process","There is no Tracks. Maybe change the option?");
1b446896 106 return;
107 }
108
109 if ( RunCoherencyCheck() )
110 {
01725374 111 Error("Process",
1b446896 112 "Coherency check not passed. Maybe change the option?\n");
113 return;
114 }
dc2c3f36 115 if (nonid) ProcessTracksAndParticlesNonIdentAnal();
116 else ProcessTracksAndParticles();
1b446896 117 return;
118 }
119
120 if(oT)
121 {
7836ee94 122 if (fReader->GetNumberOfTrackEvents() <1)
123 {
124 Error("Process","There is no data to analyze.");
125 return;
126 }
dc2c3f36 127 if (nonid) ProcessTracksNonIdentAnal();
128 else ProcessTracks();
1b446896 129 return;
130 }
131
132 if(oP)
133 {
7836ee94 134 if (fReader->GetNumberOfPartEvents() <1)
135 {
136 Error("Process","There is no data to analyze.");
137 return;
138 }
dc2c3f36 139 if (nonid) ProcessParticlesNonIdentAnal();
140 else ProcessParticles();
141
142// cout<<"NON ID"<<endl;
143// ProcessParticlesNonIdentAnal();
144// cout<<"\n\n\n NORMAL"<<endl;
145// ProcessParticles();
146
1b446896 147 return;
148 }
149
150}
151
152/*************************************************************************************/
491d1b5d 153
1b446896 154void AliHBTAnalysis::ProcessTracksAndParticles()
155{
156
157//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
158//the loops are splited
159
160
161 AliHBTParticle * part1, * part2;
162 AliHBTParticle * track1, * track2;
163
164 AliHBTEvent * trackEvent, *partEvent;
165 AliHBTEvent * trackEvent2,*partEvent2;
166
167// Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
168
169 Int_t Nev = fReader->GetNumberOfTrackEvents();
170
171 /***************************************/
172 /****** Looping same events ********/
173 /****** filling numerators ********/
174 /***************************************/
175 AliHBTPair * trackpair = new AliHBTPair();
176 AliHBTPair * partpair = new AliHBTPair();
491d1b5d 177
178 AliHBTPair * tmptrackpair;//temprary pointers to pairs
179 AliHBTPair * tmppartpair;
1b446896 180
7a2c8238 181
182
1b446896 183 for (Int_t i = 0;i<Nev;i++)
184 {
185 partEvent= fReader->GetParticleEvent(i);
186 trackEvent = fReader->GetTrackEvent(i);
187
188 if (!partEvent) continue;
189
190 //N = 0;
191
192 for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
193 {
491d1b5d 194 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
1b446896 195
196 part1= partEvent->GetParticle(j);
197 track1= trackEvent->GetParticle(j);
198
199 for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
200 {
201 part2= partEvent->GetParticle(k);
202 partpair->SetParticles(part1,part2);
203
204 track2= trackEvent->GetParticle(k);
205 trackpair->SetParticles(track1,track2);
206
207 if(fPairCut->Pass(partpair) || (fPairCut->Pass(trackpair))) //check pair cut
208 { //do not meets crietria of the pair cut, try with swaped pairs
209 if( ( fPairCut->Pass(partpair->GetSwapedPair()) ) || ( fPairCut->Pass(trackpair->GetSwapedPair()) ) )
210 continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
211 else
212 { //swaped pair meets all the criteria
491d1b5d 213 tmppartpair = partpair->GetSwapedPair();
214 tmptrackpair = trackpair->GetSwapedPair();
215
1b446896 216 }
217 }
491d1b5d 218 else
219 {//meets criteria of the pair cut
220 tmptrackpair = trackpair;
221 tmppartpair = partpair;
222 }
1b446896 223 UInt_t ii;
224 for(ii = 0;ii<fNParticleFunctions;ii++)
491d1b5d 225 fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
1b446896 226
227 for(ii = 0;ii<fNTrackFunctions;ii++)
491d1b5d 228 fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
1b446896 229
230 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
491d1b5d 231 fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair,tmppartpair);
1b446896 232 }
233 }
234 }
235
dc2c3f36 236
1b446896 237 /***** Filling denominators *********/
238 /***************************************/
7a2c8238 239 for (Int_t i = 0;i<Nev-1;i++) //In each event (but last) ....
1b446896 240 {
241
242 partEvent= fReader->GetParticleEvent(i);
243 if (!partEvent) continue;
244
245 trackEvent = fReader->GetTrackEvent(i);
246
247// N=0;
248
249 for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
250 {
1b446896 251
252 part1= partEvent->GetParticle(j);
253
254 track1= trackEvent->GetParticle(j);
255
1b446896 256 Int_t NNN;
7a2c8238 257
258 if ( ((i+fBufferSize) >= Nev) ||( fBufferSize < 0) ) //if buffer size is negative
259 //or current event+buffersize is greater
260 //than max nuber of events
261 {
262 NNN = Nev; //set the max event number
263 }
264 else
265 {
266 NNN = i+fBufferSize; //set the current event number + fBufferSize
267 }
1b446896 268
269 for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
270 {
271
272 partEvent2= fReader->GetParticleEvent(k);
273 if (!partEvent2) continue;
274
275 trackEvent2 = fReader->GetTrackEvent(k);
276
491d1b5d 277 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
1b446896 278
279 for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
280 {
281
282 // if (N>MAXCOMB) break;
283
284 part2= partEvent2->GetParticle(l);
285 partpair->SetParticles(part1,part2);
286
287 track2= trackEvent2->GetParticle(l);
288 trackpair->SetParticles(track1,track2);
289
290 if( (fPairCut->Pass(partpair)) || (fPairCut->Pass(trackpair)) ) //check pair cut
291 { //do not meets crietria of the
292 if( ( fPairCut->Pass(partpair->GetSwapedPair()) ) || ( fPairCut->Pass(trackpair->GetSwapedPair()) ) )
293 continue;
294 else
295 {
491d1b5d 296 tmppartpair = partpair->GetSwapedPair();
297 tmptrackpair = trackpair->GetSwapedPair();
1b446896 298 }
299 }
491d1b5d 300 else
301 {//meets criteria of the pair cut
302 tmptrackpair = trackpair;
303 tmppartpair = partpair;
304 }
1b446896 305 UInt_t ii;
306 for(ii = 0;ii<fNParticleFunctions;ii++)
491d1b5d 307 fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
1b446896 308
309 for(ii = 0;ii<fNTrackFunctions;ii++)
491d1b5d 310 fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
1b446896 311
312 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
491d1b5d 313 fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair);
1b446896 314
315
316 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
317 }//for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
318 }
319
320 }
321
322 /***************************************/
323
324
325}
326/*************************************************************************************/
327
328void AliHBTAnalysis::ProcessTracks()
329{
330 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
331//the loops are splited
332 AliHBTParticle * track1, * track2;
333
334 AliHBTEvent * trackEvent;
335 AliHBTEvent * trackEvent2;
336
337// Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
338
339 Int_t Nev = fReader->GetNumberOfTrackEvents();
340
341 /***************************************/
342 /****** Looping same events ********/
343 /****** filling numerators ********/
344 /***************************************/
345 AliHBTPair * trackpair = new AliHBTPair();
491d1b5d 346 AliHBTPair * tmptrackpair; //temporary pointer
1b446896 347
348 for (Int_t i = 0;i<Nev;i++)
349 {
350 trackEvent = fReader->GetTrackEvent(i);
351 if (!trackEvent) continue;
352 //N = 0;
353
354 for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
355 {
491d1b5d 356 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
1b446896 357
358 track1= trackEvent->GetParticle(j);
359
360 for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
361 {
362 track2= trackEvent->GetParticle(k);
363 trackpair->SetParticles(track1,track2);
364 if(fPairCut->Pass(trackpair)) //check pair cut
365 { //do not meets crietria of the
366 if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
491d1b5d 367 else tmptrackpair = trackpair->GetSwapedPair();
368 }
369 else
370 {
371 tmptrackpair = trackpair;
1b446896 372 }
1b446896 373 UInt_t ii;
374
375 for(ii = 0;ii<fNTrackFunctions;ii++)
491d1b5d 376 fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
1b446896 377
378
379 }
380 }
381 }
382
383 /***************************************/
384 /***** Filling diff histogram *********/
385 /***************************************/
7a2c8238 386 for (Int_t i = 0;i<Nev-1;i++) //In each event (but last) ....
1b446896 387 {
388 trackEvent = fReader->GetTrackEvent(i);
389 if (!trackEvent) continue;
390// N=0;
391
392 for (Int_t j = 0; j< trackEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
393 {
394// if (N>MAXCOMB) break;
395
396 track1= trackEvent->GetParticle(j);
397
1b446896 398 Int_t NNN;
7a2c8238 399
400 if ( ((i+fBufferSize) >= Nev) ||( fBufferSize < 0) ) //if buffer size is negative
401 //or current event+buffersize is greater
402 //than max nuber of events
403 {
404 NNN = Nev; //set the max event number
405 }
406 else
407 {
408 NNN = i+fBufferSize; //set the current event number + fBufferSize
409 }
1b446896 410
411 for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
412 {
413
414 trackEvent2 = fReader->GetTrackEvent(k);
415 if (!trackEvent2) continue;
416
491d1b5d 417 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
1b446896 418
419 for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++) // ... on all particles
420 {
421
422 // if (N>MAXCOMB) break;
423 track2= trackEvent2->GetParticle(l);
424 trackpair->SetParticles(track1,track2);
425
426 if(fPairCut->Pass(trackpair)) //check pair cut
427 { //do not meets crietria of the
428 if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
491d1b5d 429 else tmptrackpair = trackpair->GetSwapedPair();
430 }
431 else
432 {
433 tmptrackpair = trackpair;
1b446896 434 }
435 UInt_t ii;
436 for(ii = 0;ii<fNTrackFunctions;ii++)
491d1b5d 437 fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
1b446896 438
439 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
440 }//for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
441 }
442
443 }
444
445 /***************************************/
446
447
448}
449
450/*************************************************************************************/
491d1b5d 451
1b446896 452void AliHBTAnalysis::ProcessParticles()
453{
454 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
455//the loops are splited
456 AliHBTParticle * part1, * part2;
457
458 AliHBTEvent * partEvent;
459 AliHBTEvent * partEvent2;
491d1b5d 460
461 AliHBTPair * partpair = new AliHBTPair();
462 AliHBTPair * tmppartpair; //temporary pointer to the pair
1b446896 463
464// Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
465
466 Int_t Nev = fReader->GetNumberOfPartEvents();
467
7a2c8238 468 // Nev = 1;
1b446896 469 /***************************************/
470 /****** Looping same events ********/
471 /****** filling numerators ********/
472 /***************************************/
dc2c3f36 473// gBenchmark->Start("id");
7a2c8238 474
1b446896 475 for (Int_t i = 0;i<Nev;i++)
476 {
477 partEvent= fReader->GetParticleEvent(i);
478 if (!partEvent) continue;
479 //N = 0;
480
481 for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
482 {
491d1b5d 483 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
1b446896 484
485 part1= partEvent->GetParticle(j);
486
487 for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
488 {
489 part2= partEvent->GetParticle(k);
490 partpair->SetParticles(part1,part2);
491
492 if( fPairCut->Pass(partpair) ) //check pair cut
493 { //do not meets crietria of the pair cut, try with swaped pairs
494 if( fPairCut->Pass(partpair->GetSwapedPair() ) )
495 continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
496 else
497 { //swaped pair meets all the criteria
491d1b5d 498 tmppartpair = partpair->GetSwapedPair();
1b446896 499 }
500 }
491d1b5d 501 else
502 {
503 tmppartpair = partpair;
504 }
1b446896 505
506 UInt_t ii;
507
508 for(ii = 0;ii<fNParticleFunctions;ii++)
491d1b5d 509 fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
1b446896 510 }
511 }
512 }
513
514 /***************************************/
515 /***** Filling diff histogram *********/
516 /***************************************/
7a2c8238 517 for (Int_t i = 0;i<Nev-1;i++) //In each event (but last)....
1b446896 518 {
519 partEvent= fReader->GetParticleEvent(i);
520 if (!partEvent) continue;
521
522// N=0;
523
524 for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
525 {
526// if (N>MAXCOMB) break;
527
528 part1= partEvent->GetParticle(j);
529
1b446896 530 Int_t NNN;
7a2c8238 531
532 if ( ((i+fBufferSize) >= Nev) ||( fBufferSize < 0) ) //if buffer size is negative
533 //or current event+buffersize is greater
534 //than max nuber of events
535 {
536 NNN = Nev; //set the max event number
537 }
538 else
539 {
540 NNN = i+fBufferSize; //set the current event number + fBufferSize
541 }
542
543// cout<<"NNN = "<<NNN<<endl;
1b446896 544 for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
545 {
546
547 partEvent2= fReader->GetParticleEvent(k);
548 if (!partEvent2) continue;
549
491d1b5d 550 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
1b446896 551
552 for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
553 {
554
555 // if (N>MAXCOMB) break;
556 part2= partEvent2->GetParticle(l);
557 partpair->SetParticles(part1,part2);
558
559 if(fPairCut->Pass(partpair)) //check pair cut
560 { //do not meets crietria of the
561 if( fPairCut->Pass(partpair->GetSwapedPair()) ) continue;
491d1b5d 562 else tmppartpair = partpair->GetSwapedPair();
563 }
564 else
565 {
566 tmppartpair = partpair;
1b446896 567 }
568 UInt_t ii;
569 for(ii = 0;ii<fNParticleFunctions;ii++)
491d1b5d 570 fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
1b446896 571
572 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
573 }//for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
574 }
575
576 }
577
578 /***************************************/
579
dc2c3f36 580// gBenchmark->Stop("id");
581// gBenchmark->Show("id");
1b446896 582
583}
584
585/*************************************************************************************/
586
587
491d1b5d 588void AliHBTAnalysis::WriteFunctions()
1b446896 589{
590 UInt_t ii;
591 for(ii = 0;ii<fNParticleFunctions;ii++)
592 fParticleFunctions[ii]->Write();
593
594 for(ii = 0;ii<fNTrackFunctions;ii++)
595 fTrackFunctions[ii]->Write();
596
597 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
598 fParticleAndTrackFunctions[ii]->Write();
599}
600/*************************************************************************************/
491d1b5d 601
1b446896 602void AliHBTAnalysis::SetGlobalPairCut(AliHBTPairCut* cut)
603{
604 if (cut == 0x0)
605 {
606 Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
607 }
608 delete fPairCut;
609 fPairCut = (AliHBTPairCut*)cut->Clone();
610}
611
612/*************************************************************************************/
491d1b5d 613
27b3fe5d 614void AliHBTAnalysis::AddTrackFunction(AliHBTOnePairFctn* f)
1b446896 615{
616 if (f == 0x0) return;
617 if (fNTrackFunctions == fgkFctnArraySize)
618 {
619 Error("AliHBTAnalysis::AddTrackFunction","Can not add this function, not enough place in the array.");
620 }
621 fTrackFunctions[fNTrackFunctions] = f;
622 fNTrackFunctions++;
623}
491d1b5d 624/*************************************************************************************/
625
27b3fe5d 626void AliHBTAnalysis::AddParticleFunction(AliHBTOnePairFctn* f)
1b446896 627{
628 if (f == 0x0) return;
629
630 if (fNParticleFunctions == fgkFctnArraySize)
631 {
632 Error("AliHBTAnalysis::AddParticleFunction","Can not add this function, not enough place in the array.");
633 }
634 fParticleFunctions[fNParticleFunctions] = f;
635 fNParticleFunctions++;
636
637
638}
27b3fe5d 639void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTTwoPairFctn* f)
1b446896 640{
641 if (f == 0x0) return;
642 if (fNParticleAndTrackFunctions == fgkFctnArraySize)
643 {
644 Error("AliHBTAnalysis::AddParticleAndTrackFunction","Can not add this function, not enough place in the array.");
645 }
646 fParticleAndTrackFunctions[fNParticleAndTrackFunctions] = f;
647 fNParticleAndTrackFunctions++;
648}
649
650
651/*************************************************************************************/
652
653
654/*************************************************************************************/
491d1b5d 655
1b446896 656Bool_t AliHBTAnalysis::RunCoherencyCheck()
657{
658 //Checks if both HBTRuns are similar
659 //return true if error found
660 //if they seem to be OK return false
661 Int_t i;
662 cout<<"Checking HBT Runs Coherency"<<endl;
663
664 cout<<"Number of events ..."; fflush(stdout);
665
666 if (fReader->GetNumberOfPartEvents() == fReader->GetNumberOfTrackEvents() ) //check whether there is the same number of events
667 {
668 cout<<"OK. "<<fReader->GetNumberOfTrackEvents()<<" found"<<endl;
669 }
670 else
671 { //if not the same - ERROR
672 Error("AliHBTAnalysis::RunCoherencyCheck()",
673 "Number of simulated events (%d) is not equal to number of reconstructed events(%d)",
674 fReader->GetNumberOfPartEvents(),fReader->GetNumberOfTrackEvents());
675 return kTRUE;
676 }
677
678 cout<<"Checking number of Particles AND Particles Types in each event ...";fflush(stdout);
679
680 AliHBTEvent *partEvent;
681 AliHBTEvent *trackEvent;
682 for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
683 {
684 partEvent= fReader->GetParticleEvent(i); //gets the "ith" event
685 trackEvent = fReader->GetTrackEvent(i);
686
687 if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
688 if ( (partEvent == 0x0) || (partEvent == 0x0) )
689 {
690 Error("AliHBTAnalysis::RunCoherencyCheck()",
691 "One event is NULL and the other one not. Event Number %d",i);
692 return kTRUE;
693 }
694
695 if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
696 {
697 Error("AliHBTAnalysis::RunCoherencyCheck()",
698 "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
699 i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
700 return kTRUE;
701 }
702 else
703 for (Int_t j = 0; j<partEvent->GetNumberOfParticles(); j++)
704 {
705 if( partEvent->GetParticle(j)->GetPdgCode() != trackEvent->GetParticle(j)->GetPdgCode() )
706 {
707 Error("AliHBTAnalysis::RunCoherencyCheck()",
708 "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
709 i,j, partEvent->GetParticle(j)->GetPdgCode(),trackEvent->GetParticle(j)->GetPdgCode() );
710 return kTRUE;
711
712 }
713 }
714 }
715 cout<<" Done"<<endl;
716 cout<<" Everything looks OK"<<endl;
717 return kFALSE;
718}
719
dc2c3f36 720/*************************************************************************************/
721
722void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
723{
724 AliHBTParticle * part1, * part2;
725 AliHBTParticle * track1, * track2;
726
727 AliHBTEvent * trackEvent1,*partEvent1;
728 AliHBTEvent * trackEvent2,*partEvent2;
729 AliHBTEvent * trackEvent3,*partEvent3;
730
731 AliHBTEvent * rawtrackEvent, *rawpartEvent;
732// AliHBTEvent * rawtrackEvent2,*rawpartEvent2;
733
734
735 Int_t Nev = fReader->GetNumberOfTrackEvents();
736
737 AliHBTPair * trackpair = new AliHBTPair();
738 AliHBTPair * partpair = new AliHBTPair();
739
740 TList tbuffer;
741 TList pbuffer;
742 Int_t ninbuffer = 0;
743
744 trackEvent1 = new AliHBTEvent();
745 partEvent1 = new AliHBTEvent();
746 trackEvent1->SetOwner(kFALSE);
747 partEvent1->SetOwner(kFALSE);;
748
749 cout<<"**************************************"<<endl;
750 cout<<"***** NON IDENT MODE ****************"<<endl;
751 cout<<"**************************************"<<endl;
752 for (Int_t i = 0;i<Nev;i++)
753 {
754 rawpartEvent = fReader->GetParticleEvent(i);
755 rawtrackEvent = fReader->GetTrackEvent(i);
756 if ( (rawpartEvent == 0x0) || (rawtrackEvent == 0x0) ) continue;//in case of any error
757
758 /********************************/
759 /* Filtering out */
760 /********************************/
761 if (ninbuffer > fBufferSize)
762 {
763 partEvent2 = (AliHBTEvent*)pbuffer.Remove(pbuffer.Last()); //remove from the end to be reset, filled and put on the beginning
764 trackEvent2 = (AliHBTEvent*)tbuffer.Remove(tbuffer.Last());
765 ninbuffer--;
766 }
767 else
768 {
769 partEvent2 = new AliHBTEvent();
770 trackEvent2 = new AliHBTEvent();
771 partEvent2->SetOwner(kFALSE);
772 trackEvent2->SetOwner(kFALSE);
773 }
774 FilterOut(partEvent1, partEvent2, rawpartEvent, trackEvent1, trackEvent2, rawtrackEvent);
775
776 for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
777 {
778 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
779
780 part1= partEvent1->GetParticle(j);
781 track1= trackEvent1->GetParticle(j);
782
783
784 /***************************************/
785 /****** filling numerators ********/
786 /****** (particles from event2) ********/
787 /***************************************/
788 for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++)
789 {
790 part2= partEvent2->GetParticle(k);
791 partpair->SetParticles(part1,part2);
792
793 track2= trackEvent2->GetParticle(k);
794 trackpair->SetParticles(track1,track2);
795
796
797 if( (fPairCut->PassPairProp(partpair)) || (fPairCut->PassPairProp(trackpair))) //check pair cut
798 { //do not meets crietria of the pair cut
799 continue;
800 }
801 else
802 {//meets criteria of the pair cut
803 UInt_t ii;
804 for(ii = 0;ii<fNParticleFunctions;ii++)
805 fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
806
807 for(ii = 0;ii<fNTrackFunctions;ii++)
808 fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
809
810 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
811 fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(trackpair,partpair);
812 }
813 }
814 /***************************************/
815 /***** Filling denominators *********/
816 /***************************************/
817 TIter piter(&pbuffer);
818 TIter titer(&tbuffer);
819 Int_t nmonitor = 0;
820
821 while ( (partEvent3 = (AliHBTEvent*)piter.Next()) != 0x0)
822 {
823 trackEvent3 = (AliHBTEvent*)titer.Next();
824 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i - (++nmonitor)<<endl;
825
826 for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
827 {
828 part2= partEvent3->GetParticle(k);
829 partpair->SetParticles(part1,part2);
1b446896 830
dc2c3f36 831 track2= trackEvent3->GetParticle(k);
832 trackpair->SetParticles(track1,track2);
833
834
835 if( (fPairCut->PassPairProp(partpair)) || (fPairCut->PassPairProp(trackpair))) //check pair cut
836 { //do not meets crietria of the pair cut
837 continue;
838 }
839 else
840 {//meets criteria of the pair cut
841 UInt_t ii;
842 for(ii = 0;ii<fNParticleFunctions;ii++)
843 fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
844
845 for(ii = 0;ii<fNTrackFunctions;ii++)
846 fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
847
848 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
849 fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(trackpair,partpair);
850 }
851 }// for particles event2
852 }//while event2
853 }//for over particles in event1
854
855 pbuffer.AddFirst(partEvent2);
856 tbuffer.AddFirst(trackEvent2);
857 ninbuffer++;
858
859 }//end of loop over events (1)
860
861 pbuffer.SetOwner(); //to clean stored events by the way of buffer destruction
862 tbuffer.SetOwner();
863}
1b446896 864/*************************************************************************************/
865
dc2c3f36 866void AliHBTAnalysis::ProcessTracksNonIdentAnal()
867{
868 AliHBTParticle * track1, * track2;
869
870 AliHBTEvent * trackEvent1;
871 AliHBTEvent * trackEvent2;
872 AliHBTEvent * trackEvent3;
873
874 AliHBTEvent * rawtrackEvent;
875// AliHBTEvent * rawtrackEvent2,*rawpartEvent2;
876
877
878 Int_t Nev = fReader->GetNumberOfTrackEvents();
879
880 AliHBTPair * trackpair = new AliHBTPair();
881
882 TList tbuffer;
883 Int_t ninbuffer = 0;
884
885 trackEvent1 = new AliHBTEvent();
886 trackEvent1->SetOwner(kFALSE);
887
888 cout<<"**************************************"<<endl;
889 cout<<"***** NON IDENT MODE ****************"<<endl;
890 cout<<"**************************************"<<endl;
891 for (Int_t i = 0;i<Nev;i++)
892 {
893 rawtrackEvent = fReader->GetTrackEvent(i);
894 if (rawtrackEvent == 0x0) continue;//in case of any error
895
896 /********************************/
897 /* Filtering out */
898 /********************************/
899 if (ninbuffer > fBufferSize)
900 {
901 trackEvent2 = (AliHBTEvent*)tbuffer.Remove(tbuffer.Last());
902 ninbuffer--;
903 }
904 else
905 {
906 trackEvent2 = new AliHBTEvent();
907 trackEvent2->SetOwner(kFALSE);
908 }
909 FilterOut(trackEvent1, trackEvent2, rawtrackEvent);
910
911 for (Int_t j = 0; j<trackEvent1->GetNumberOfParticles() ; j++)
912 {
913 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
914
915 track1= trackEvent1->GetParticle(j);
916
917
918 /***************************************/
919 /****** filling numerators ********/
920 /****** (particles from event2) ********/
921 /***************************************/
922 for (Int_t k = 0; k < trackEvent2->GetNumberOfParticles() ; k++)
923 {
924 track2= trackEvent2->GetParticle(k);
925 trackpair->SetParticles(track1,track2);
926
927
928 if( fPairCut->PassPairProp(trackpair)) //check pair cut
929 { //do not meets crietria of the pair cut
930 continue;
931 }
932 else
933 {//meets criteria of the pair cut
934 UInt_t ii;
935 for(ii = 0;ii<fNTrackFunctions;ii++)
936 fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
937 }
938 }
939 /***************************************/
940 /***** Filling denominators *********/
941 /***************************************/
942 TIter titer(&tbuffer);
943 Int_t nmonitor = 0;
944
945 while ( (trackEvent3 = (AliHBTEvent*)titer.Next()) != 0x0)
946 {
947
948 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i - (++nmonitor)<<endl;
949
950 for (Int_t k = 0; k < trackEvent3->GetNumberOfParticles() ; k++)
951 {
952
953 track2= trackEvent3->GetParticle(k);
954 trackpair->SetParticles(track1,track2);
955
956
957 if( fPairCut->PassPairProp(trackpair)) //check pair cut
958 { //do not meets crietria of the pair cut
959 continue;
960 }
961 else
962 {//meets criteria of the pair cut
963 UInt_t ii;
964 for(ii = 0;ii<fNTrackFunctions;ii++)
965 fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
966
967 }
968 }// for particles event2
969 }//while event2
970 }//for over particles in event1
971
972 tbuffer.AddFirst(trackEvent2);
973 ninbuffer++;
974
975 }//end of loop over events (1)
976
977 tbuffer.SetOwner();
978}
979/*************************************************************************************/
980
981void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
982{
983 AliHBTParticle * part1 = 0x0, * part2 = 0x0;
984
985 AliHBTEvent * partEvent1 = 0x0;
986 AliHBTEvent * partEvent2 = 0x0;
987 AliHBTEvent * partEvent3 = 0x0;
988
989 AliHBTEvent * rawpartEvent = 0x0;
990
991 Int_t Nev = fReader->GetNumberOfPartEvents();
992
993 AliHBTPair * partpair = new AliHBTPair();
994
995 TList pbuffer;
996 Int_t ninbuffer = 0;
997
998 partEvent1 = new AliHBTEvent();
999 partEvent1->SetOwner(kFALSE);;
1000
1001 cout<<"**************************************"<<endl;
1002 cout<<"***** PART NON IDENT MODE ******"<<endl;
1003 cout<<"**************************************"<<endl;
1004
1005// gBenchmark->Start("non_id");
1006 for (Int_t i = 0;i<Nev;i++)
1007 {
1008 rawpartEvent = fReader->GetParticleEvent(i);
1009 if ( rawpartEvent == 0x0 ) continue;//in case of any error
1010
1011 /********************************/
1012 /* Filtering out */
1013 /********************************/
1014 if (ninbuffer > fBufferSize)
1015 {
1016 partEvent2 = (AliHBTEvent*)pbuffer.Remove(pbuffer.Last()); //remove from the end to be reset, filled and put on the beginning
1017 ninbuffer--;
1018 }
1019 else
1020 {
1021 partEvent2 = new AliHBTEvent();
1022 partEvent2->SetOwner(kFALSE);
1023 }
1024 FilterOut(partEvent1, partEvent2, rawpartEvent);
1025
1026 for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
1027 {
1028 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
1029
1030 part1= partEvent1->GetParticle(j);
1031
1032
1033 /***************************************/
1034 /****** filling numerators ********/
1035 /****** (particles from event2) ********/
1036 /***************************************/
1037 for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++)
1038 {
1039 part2= partEvent2->GetParticle(k);
1040 partpair->SetParticles(part1,part2);
1041
1042 if(fPairCut->PassPairProp(partpair) ) //check pair cut
1043 { //do not meets crietria of the pair cut
1044 continue;
1045 }
1046 else
1047 {//meets criteria of the pair cut
1048 UInt_t ii;
1049 for(ii = 0;ii<fNParticleFunctions;ii++)
1050 fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
1051
1052 }
1053 }
1054 /***************************************/
1055 /***** Filling denominators *********/
1056 /***************************************/
1057 TIter piter(&pbuffer);
1058 Int_t nmonitor = 0;
1059
1060 while ( (partEvent3 = (AliHBTEvent*)piter.Next()) != 0x0)
1061 {
1062 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i - (++nmonitor)<<endl;
1063
1064 for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
1065 {
1066 part2= partEvent3->GetParticle(k);
1067 partpair->SetParticles(part1,part2);
1068
1069
1070 if(fPairCut->PassPairProp(partpair) ) //check pair cut
1071 { //do not meets crietria of the pair cut
1072 continue;
1073 }
1074 else
1075 {//meets criteria of the pair cut
1076 UInt_t ii;
1077 for(ii = 0;ii<fNParticleFunctions;ii++)
1078 fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
1079
1080 }
1081 }// for particles event2
1082 }//while event2
1083 }//for over particles in event1
1084
1085 pbuffer.AddFirst(partEvent2);
1086 ninbuffer++;
1087
1088 }//end of loop over events (1)
1089
1090// gBenchmark->Stop("non_id");
1091// gBenchmark->Show("non_id");
1092 pbuffer.SetOwner();//to delete stered events.
1093}
1094
1095/*************************************************************************************/
1096void AliHBTAnalysis::FilterOut(AliHBTEvent* outpart1, AliHBTEvent* outpart2, AliHBTEvent* inpart,
1097 AliHBTEvent* outtrack1, AliHBTEvent* outtrack2, AliHBTEvent* intrack)
1098{
1099 //Puts particles accepted as a first particle by global cut in out1
1100 //and as a second particle in out2
1101
1102 AliHBTParticle* part, *track;
1103
1104 outpart1->Reset();
1105 outpart2->Reset();
1106 outtrack1->Reset();
1107 outtrack2->Reset();
1108
1109 AliHBTParticleCut *cut1 = fPairCut->GetFirstPartCut();
1110 AliHBTParticleCut *cut2 = fPairCut->GetSecondPartCut();
1111
1112 Bool_t in1, in2;
1113
1114 for (Int_t i = 0; i < inpart->GetNumberOfParticles(); i++)
1115 {
1116 in1 = in2 = kTRUE;
1117 part = inpart->GetParticle(i);
1118 track = intrack->GetParticle(i);
1119
1120 if ( (cut1->Pass(part)) || ( cut1->Pass(track) ) ) in1 = kFALSE; //if any, part or track, is rejected by cut1, in1 is false
1121 if ( (cut2->Pass(part)) || ( cut2->Pass(track) ) ) in2 = kFALSE; //if any, part or track, is rejected by cut2, in2 is false
1122
1123 if (gDebug)//to be removed in real analysis
1124 if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1125 {
1126 //Particle accpted by both cuts
1127 Error("FilterOut","Particle accepted by both cuts");
1128 continue;
1129 }
1130
1131 if (in1)
1132 {
1133 outpart1->AddParticle(part);
1134 outtrack1->AddParticle(track);
1135 continue;
1136 }
1137
1138 if (in2)
1139 {
1140 outpart2->AddParticle(part);
1141 outtrack2->AddParticle(track);
1142 continue;
1143 }
1144 }
1145
1146}
1147
1b446896 1148
1149/*************************************************************************************/
dc2c3f36 1150void AliHBTAnalysis::FilterOut(AliHBTEvent* out1, AliHBTEvent* out2, AliHBTEvent* in)
1151{
1152 //Puts particles accepted as a first particle by global cut in out1
1153 //and as a second particle in out2
1154 AliHBTParticle* part;
1155
1156 out1->Reset();
1157 out2->Reset();
1158
1159 AliHBTParticleCut *cut1 = fPairCut->GetFirstPartCut();
1160 AliHBTParticleCut *cut2 = fPairCut->GetSecondPartCut();
1161
1162 Bool_t in1, in2;
1163
1164 for (Int_t i = 0; i < in->GetNumberOfParticles(); i++)
1165 {
1166 in1 = in2 = kTRUE;
1167 part = in->GetParticle(i);
1168
1169 if ( cut1->Pass(part) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
1170 if ( cut2->Pass(part) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
1171
1172 if (gDebug)//to be removed in real analysis
1173 if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1174 {
1175 //Particle accpted by both cuts
1176 Error("FilterOut","Particle accepted by both cuts");
1177 continue;
1178 }
1b446896 1179
dc2c3f36 1180 if (in1)
1181 {
1182 out1->AddParticle(part);
1183 continue;
1184 }
1185
1186 if (in2)
1187 {
1188 out2->AddParticle(part);
1189 continue;
1190 }
1191 }
1192}
1193/*************************************************************************************/
1b446896 1194
dc2c3f36 1195Bool_t AliHBTAnalysis::IsNonIdentAnalysis()
1196{
1197 //checks if it is possible to use special analysis for non identical particles
1198 //it means - in global pair cut first particle id is different than second one
1199 //and both are different from 0
1200 //in the future is possible to perform more sophisticated check
1201 //if cuts have excluding requirements
1202
1203 if (fPairCut->IsEmpty()) return kFALSE;
1204 if (fPairCut->GetFirstPartCut()->IsEmpty()) return kFALSE;
1205 if (fPairCut->GetSecondPartCut()->IsEmpty()) return kFALSE;
1206
1207 Int_t id1 = fPairCut->GetFirstPartCut()->GetPID();
1208 Int_t id2 = fPairCut->GetSecondPartCut()->GetPID();
1209
1210 if ( (id1==0) || (id2==0) || (id1==id2) ) return kFALSE;
1211
1212 return kTRUE;
1213}