1 #include "AliHBTAnalysis.h"
3 #include "AliHBTEvent.h"
4 #include "AliHBTReader.h"
5 #include "AliHBTParticle.h"
6 #include "AliHBTParticleCut.h"
7 #include "AliHBTPair.h"
8 #include "AliHBTPairCut.h"
9 #include "AliHBTFunction.h"
10 #include "AliHBTMonitorFunction.h"
11 #include "AliHBTEventBuffer.h"
13 #include <TBenchmark.h>
16 //_________________________________________________________
17 ///////////////////////////////////////////////////////////
19 //Central Object Of HBTAnalyser:
20 //This class performs main looping within HBT Analysis
21 //User must plug a reader of Type AliHBTReader
22 //User plugs in coorelation and monitor functions
23 //as well as monitor functions
25 //HBT Analysis Tool, which is integral part of AliRoot,
26 //ALICE Off-Line framework:
28 //Piotr.Skowronski@cern.ch
29 //more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
31 //_________________________________________________________
33 ClassImp(AliHBTAnalysis)
35 const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100;
36 const UInt_t AliHBTAnalysis::fgkDefaultMixingInfo = 1000;
37 const Int_t AliHBTAnalysis::fgkDefaultBufferSize = 5;
39 AliHBTAnalysis::AliHBTAnalysis():
42 fNParticleFunctions(0),
43 fNParticleAndTrackFunctions(0),
44 fNTrackMonitorFunctions(0),
45 fNParticleMonitorFunctions(0),
46 fNParticleAndTrackMonitorFunctions(0),
48 fDisplayMixingInfo(fgkDefaultMixingInfo),
52 fTrackFunctions = new AliHBTOnePairFctn* [fgkFctnArraySize];
53 fParticleFunctions = new AliHBTOnePairFctn* [fgkFctnArraySize];
54 fParticleAndTrackFunctions = new AliHBTTwoPairFctn* [fgkFctnArraySize];
56 fParticleMonitorFunctions = new AliHBTMonOneParticleFctn* [fgkFctnArraySize];
57 fTrackMonitorFunctions = new AliHBTMonOneParticleFctn* [fgkFctnArraySize];
58 fParticleAndTrackMonitorFunctions = new AliHBTMonTwoParticleFctn* [fgkFctnArraySize];
60 fPairCut = new AliHBTEmptyPairCut();//empty cut - accepts all particles
62 /*************************************************************************************/
64 AliHBTAnalysis::AliHBTAnalysis(const AliHBTAnalysis& in):
68 fNParticleFunctions(0),
69 fNParticleAndTrackFunctions(0),
70 fNTrackMonitorFunctions(0),
71 fNParticleMonitorFunctions(0),
72 fNParticleAndTrackMonitorFunctions(0),
74 fParticleFunctions(0x0),
75 fParticleAndTrackFunctions(0x0),
76 fParticleMonitorFunctions(0x0),
77 fTrackMonitorFunctions(0x0),
78 fParticleAndTrackMonitorFunctions(0x0),
80 fBufferSize(fgkDefaultBufferSize),
81 fDisplayMixingInfo(fgkDefaultMixingInfo),
85 Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
87 /*************************************************************************************/
88 const AliHBTAnalysis& AliHBTAnalysis::operator=(const AliHBTAnalysis& /*right*/)
91 Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
94 /*************************************************************************************/
95 AliHBTAnalysis::~AliHBTAnalysis()
98 //note that we do not delete functions itself
99 // they should be deleted by whom where created
100 //we only store pointers, and we use "new" only for pointers array
103 delete [] fTrackFunctions;
104 delete [] fParticleFunctions;
105 delete [] fParticleAndTrackFunctions;
107 delete [] fParticleMonitorFunctions;
108 delete [] fTrackMonitorFunctions;
109 delete [] fParticleAndTrackMonitorFunctions;
111 delete fPairCut; // always have an copy of an object - we create we dstroy
114 /*************************************************************************************/
116 void AliHBTAnalysis::DeleteFunctions()
118 //Deletes all functions added to analysis
120 for(ii = 0;ii<fNParticleFunctions;ii++)
121 delete fParticleFunctions[ii];
122 fNParticleFunctions = 0;
124 for(ii = 0;ii<fNTrackFunctions;ii++)
125 delete fTrackFunctions[ii];
126 fNTrackFunctions = 0;
128 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
129 delete fParticleAndTrackFunctions[ii];
130 fNParticleAndTrackFunctions = 0;
132 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
133 delete fParticleMonitorFunctions[ii];
134 fNParticleMonitorFunctions = 0;
136 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
137 delete fTrackMonitorFunctions[ii];
138 fNTrackMonitorFunctions = 0;
140 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
141 delete fParticleAndTrackMonitorFunctions[ii];
142 fNParticleAndTrackMonitorFunctions = 0;
144 /*************************************************************************************/
146 void AliHBTAnalysis::Init()
148 //Initializeation method
149 //calls Init for all functions
151 for(ii = 0;ii<fNParticleFunctions;ii++)
152 fParticleFunctions[ii]->Init();
154 for(ii = 0;ii<fNTrackFunctions;ii++)
155 fTrackFunctions[ii]->Init();
157 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
158 fParticleAndTrackFunctions[ii]->Init();
160 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
161 fParticleMonitorFunctions[ii]->Init();
163 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
164 fTrackMonitorFunctions[ii]->Init();
166 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
167 fParticleAndTrackMonitorFunctions[ii]->Init();
169 /*************************************************************************************/
171 void AliHBTAnalysis::ResetFunctions()
173 //In case fOwner is true, deletes all functions
174 //in other case, just set number of analysis to 0
175 if (fIsOwner) DeleteFunctions();
178 fNParticleFunctions = 0;
179 fNTrackFunctions = 0;
180 fNParticleAndTrackFunctions = 0;
181 fNParticleMonitorFunctions = 0;
182 fNTrackMonitorFunctions = 0;
183 fNParticleAndTrackMonitorFunctions = 0;
186 /*************************************************************************************/
188 void AliHBTAnalysis::Process(Option_t* option)
190 //default option = "TracksAndParticles"
191 //Main method of the HBT Analysis Package
192 //It triggers reading with the global cut (default is an empty cut)
193 //Than it checks options and data which are read
194 //if everything is OK, then it calls one of the looping methods
195 //depending on tfReaderhe option
196 //These methods differs on what thay are looping on
199 //--------------------------------------------------------------------
200 //ProcessTracksAndParticles - "TracksAndParticles"
202 // it loops over both, tracks(reconstructed) and particles(simulated)
203 // all function gethered in all 3 lists are called for each (double)pair
205 //ProcessTracks - "Tracks"
206 // it loops only on tracks(reconstructed),
207 // functions ONLY from fTrackFunctions list are called
209 //ProcessParticles - "Particles"
210 // it loops only on particles(simulated),
211 // functions ONLY from fParticleAndTrackFunctions list are called
216 Error("Process","The reader is not set");
220 const char *oT = strstr(option,"Tracks");
221 const char *oP = strstr(option,"Particles");
223 Bool_t nonid = IsNonIdentAnalysis();
229 if (nonid) ProcessTracksAndParticlesNonIdentAnal();
230 else ProcessTracksAndParticles();
236 if (nonid) ProcessTracksNonIdentAnal();
237 else ProcessTracks();
243 if (nonid) ProcessParticlesNonIdentAnal();
244 else ProcessParticles();
249 /*************************************************************************************/
251 void AliHBTAnalysis::ProcessTracksAndParticles()
253 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
254 //the loops are splited
257 AliHBTParticle * part1, * part2;
258 AliHBTParticle * track1, * track2;
260 AliHBTEvent * trackEvent, *partEvent;
261 AliHBTEvent * trackEvent1 = 0x0,*partEvent1 = 0x0;
262 AliHBTEvent * trackEvent2,*partEvent2;
264 // Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
266 // Int_t nev = fReader->GetNumberOfTrackEvents();
267 AliHBTPair * trackpair = new AliHBTPair();
268 AliHBTPair * partpair = new AliHBTPair();
270 AliHBTPair * tmptrackpair;//temprary pointers to pairs
271 AliHBTPair * tmppartpair;
273 AliHBTEventBuffer partbuffer(fBufferSize);
274 AliHBTEventBuffer trackbuffer(fBufferSize);
279 Bool_t nocorrfctns = (fNParticleFunctions == 0) && (fNTrackFunctions == 0) && (fNParticleAndTrackFunctions == 0);
283 while (fReader->Next() == kFALSE)
286 partEvent= fReader->GetParticleEvent();
287 trackEvent = fReader->GetTrackEvent();
289 if ( !partEvent || !trackEvent )
291 Error("ProcessTracksAndParticles","Can not get event");
295 if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
297 Fatal("ProcessTracksAndParticles",
298 "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
299 i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
303 if(partEvent1 == 0x0)
305 partEvent1 = new AliHBTEvent();
306 partEvent1->SetOwner(kTRUE);
308 trackEvent1 = new AliHBTEvent();
309 trackEvent1->SetOwner(kTRUE);
314 trackEvent1->Reset();
317 for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
319 /***************************************/
320 /****** Looping same events ********/
321 /****** filling numerators ********/
322 /***************************************/
323 if ( (j%fDisplayMixingInfo) == 0)
324 Info("ProcessTracksAndParticles",
325 "Mixing particle %d from event %d with particles from event %d",j,i,i);
327 part1= partEvent->GetParticle(j);
328 track1= trackEvent->GetParticle(j);
330 //PID imperfections ???
331 // if( part1->GetPdgCode() != track1->GetPdgCode() )
333 // Fatal("ProcessTracksAndParticles",
334 // "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
335 // i,j, part1->GetPdgCode(),track1->GetPdgCode() );
339 Bool_t firstcut = fPairCut->GetFirstPartCut()->Pass(part1);
340 if (fBufferSize != 0)
341 if ( (firstcut == kFALSE) || (fPairCut->GetSecondPartCut()->Pass(part1) == kFALSE) )
343 //accepted by any cut
344 // we have to copy because reader keeps only one event
346 partEvent1->AddParticle(new AliHBTParticle(*part1));
347 trackEvent1->AddParticle(new AliHBTParticle(*track1));
350 if (firstcut) continue;
352 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
353 fParticleMonitorFunctions[ii]->Process(part1);
354 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
355 fTrackMonitorFunctions[ii]->Process(track1);
356 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
357 fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
359 if (nocorrfctns) continue;
361 for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
363 part2= partEvent->GetParticle(k);
364 if (part1->GetUID() == part2->GetUID()) continue;
365 partpair->SetParticles(part1,part2);
367 track2= trackEvent->GetParticle(k);
368 trackpair->SetParticles(track1,track2);
370 if(fPairCut->Pass(trackpair) ) //check pair cut
371 { //do not meets crietria of the pair cut, try with swapped pairs
372 if( fPairCut->Pass(trackpair->GetSwapedPair()) )
373 continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
375 { //swaped pair meets all the criteria
376 tmppartpair = partpair->GetSwapedPair();
377 tmptrackpair = trackpair->GetSwapedPair();
381 {//meets criteria of the pair cut
382 tmptrackpair = trackpair;
383 tmppartpair = partpair;
385 for(ii = 0;ii<fNParticleFunctions;ii++)
386 fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
388 for(ii = 0;ii<fNTrackFunctions;ii++)
389 fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
391 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
392 fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair,tmppartpair);
393 //end of 2nd loop over particles from the same event
394 }//for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
396 /***************************************/
397 /***** Filling denominators *********/
398 /***************************************/
399 if (fBufferSize == 0) continue;
401 partbuffer.ResetIter();
402 trackbuffer.ResetIter();
404 while (( partEvent2 = partbuffer.Next() ))
406 trackEvent2 = trackbuffer.Next();
409 if ( (j%fDisplayMixingInfo) == 0)
410 Info("ProcessTracksAndParticles",
411 "Mixing particle %d from event %d with particles from event %d",j,i,i-m);
413 for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
415 part2= partEvent2->GetParticle(l);
416 partpair->SetParticles(part1,part2);
418 track2= trackEvent2->GetParticle(l);
419 trackpair->SetParticles(track1,track2);
421 if( fPairCut->Pass(partpair) ) //check pair cut
422 { //do not meets crietria of the
423 if( fPairCut->Pass(partpair->GetSwapedPair()) )
427 tmppartpair = partpair->GetSwapedPair();
428 tmptrackpair = trackpair->GetSwapedPair();
432 {//meets criteria of the pair cut
433 tmptrackpair = trackpair;
434 tmppartpair = partpair;
436 for(ii = 0;ii<fNParticleFunctions;ii++)
437 fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
439 for(ii = 0;ii<fNTrackFunctions;ii++)
440 fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
442 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
443 fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair);
444 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
447 //end of loop over particles from first event
448 }//for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
449 partEvent1 = partbuffer.Push(partEvent1);
450 trackEvent1 = trackbuffer.Push(trackEvent1);
451 //end of loop over events
452 }//while (fReader->Next() == kFALSE)
454 /*************************************************************************************/
456 void AliHBTAnalysis::ProcessTracks()
458 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
459 //the loops are splited
460 AliHBTParticle * track1, * track2;
461 AliHBTEvent * trackEvent;
462 AliHBTEvent * trackEvent1 = 0x0;
463 AliHBTEvent * trackEvent2;
467 AliHBTPair * trackpair = new AliHBTPair();
468 AliHBTPair * tmptrackpair; //temporary pointer
470 AliHBTEventBuffer trackbuffer(fBufferSize);
474 while (fReader->Next() == kFALSE)
477 trackEvent = fReader->GetTrackEvent();
478 if (!trackEvent) continue;
480 if(trackEvent1 == 0x0)
482 trackEvent1 = new AliHBTEvent();
483 trackEvent1->SetOwner(kTRUE);
487 trackEvent1->Reset();
490 for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
492 /***************************************/
493 /****** Looping same events ********/
494 /****** filling numerators ********/
495 /***************************************/
496 if ( (j%fDisplayMixingInfo) == 0)
497 Info("ProcessTracks",
498 "Mixing particle %d from event %d with particles from event %d",j,i,i);
500 track1= trackEvent->GetParticle(j);
501 Bool_t firstcut = fPairCut->GetFirstPartCut()->Pass(track1);
503 if (fBufferSize != 0)
504 if ( (firstcut == kFALSE) || (fPairCut->GetSecondPartCut()->Pass(track1) == kFALSE) )
506 //accepted by any cut
507 // we have to copy because reader keeps only one event
508 trackEvent1->AddParticle(new AliHBTParticle(*track1));
511 if (firstcut) continue;
513 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
514 fTrackMonitorFunctions[ii]->Process(track1);
516 if ( fNTrackFunctions ==0 ) continue;
518 for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
520 track2= trackEvent->GetParticle(k);
521 if (track1->GetUID() == track2->GetUID()) continue;
523 trackpair->SetParticles(track1,track2);
524 if(fPairCut->Pass(trackpair)) //check pair cut
525 { //do not meets crietria of the
526 if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
527 else tmptrackpair = trackpair->GetSwapedPair();
531 tmptrackpair = trackpair;
533 for(ii = 0;ii<fNTrackFunctions;ii++)
534 fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
536 /***************************************/
537 /***** Filling denominators *********/
538 /***************************************/
540 if (fBufferSize == 0) continue;
542 trackbuffer.ResetIter();
545 while (( trackEvent2 = trackbuffer.Next() ))
548 if ( (j%fDisplayMixingInfo) == 0)
549 Info("ProcessTracks",
550 "Mixing track %d from event %d with tracks from event %d",j,i,i-m);
551 for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++) // ... on all particles
554 track2= trackEvent2->GetParticle(l);
555 trackpair->SetParticles(track1,track2);
557 if( fPairCut->Pass(trackpair) ) //check pair cut
558 { //do not meets crietria of the
559 if( fPairCut->Pass(trackpair->GetSwapedPair()) )
563 tmptrackpair = trackpair->GetSwapedPair();
567 {//meets criteria of the pair cut
568 tmptrackpair = trackpair;
571 for(ii = 0;ii<fNTrackFunctions;ii++)
572 fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
574 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
578 trackEvent1 = trackbuffer.Push(trackEvent1);
579 }//while (fReader->Next() == kFALSE)
582 /*************************************************************************************/
584 void AliHBTAnalysis::ProcessParticles()
586 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
587 //the loops are splited
588 AliHBTParticle * part1, * part2;
589 AliHBTEvent * partEvent;
590 AliHBTEvent * partEvent1 = 0x0;
591 AliHBTEvent * partEvent2;
595 AliHBTPair * partpair = new AliHBTPair();
596 AliHBTPair * tmppartpair; //temporary pointer
598 AliHBTEventBuffer partbuffer(fBufferSize);
602 while (fReader->Next() == kFALSE)
605 partEvent = fReader->GetParticleEvent();
606 if (!partEvent) continue;
608 if(partEvent1 == 0x0)
610 partEvent1 = new AliHBTEvent();
611 partEvent1->SetOwner(kTRUE);
618 for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
620 /***************************************/
621 /****** Looping same events ********/
622 /****** filling numerators ********/
623 /***************************************/
624 if ( (j%fDisplayMixingInfo) == 0)
625 Info("ProcessParticles",
626 "Mixing particle %d from event %d with particles from event %d",j,i,i);
628 part1 = partEvent->GetParticle(j);
629 Bool_t firstcut = fPairCut->GetFirstPartCut()->Pass(part1);
631 if (fBufferSize != 0) //useless in case
632 if ( (firstcut == kFALSE) || (fPairCut->GetSecondPartCut()->Pass(part1) == kFALSE) )
634 //accepted by any cut
635 // we have to copy because reader keeps only one event
636 partEvent1->AddParticle(new AliHBTParticle(*part1));
639 if (firstcut) continue;
641 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
642 fParticleMonitorFunctions[ii]->Process(part1);
644 if ( fNParticleFunctions == 0 ) continue;
646 for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
648 part2= partEvent->GetParticle(k);
649 if (part1->GetUID() == part2->GetUID()) continue;
651 partpair->SetParticles(part1,part2);
652 if(fPairCut->Pass(partpair)) //check pair cut
653 { //do not meets crietria of the
654 if( fPairCut->Pass(partpair->GetSwapedPair()) ) continue;
655 else tmppartpair = partpair->GetSwapedPair();
659 tmppartpair = partpair;
661 for(ii = 0;ii<fNParticleFunctions;ii++)
662 fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
664 /***************************************/
665 /***** Filling denominators *********/
666 /***************************************/
668 if (fBufferSize == 0) continue;
670 partbuffer.ResetIter();
673 while (( partEvent2 = partbuffer.Next() ))
676 if ( (j%fDisplayMixingInfo) == 0)
677 Info("ProcessParticles",
678 "Mixing particle %d from event %d with particles from event %d",j,i,i-m);
679 for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
682 part2= partEvent2->GetParticle(l);
683 partpair->SetParticles(part1,part2);
685 if( fPairCut->Pass(partpair) ) //check pair cut
686 { //do not meets crietria of the
687 if( fPairCut->Pass(partpair->GetSwapedPair()) )
691 tmppartpair = partpair->GetSwapedPair();
695 {//meets criteria of the pair cut
696 tmppartpair = partpair;
699 for(ii = 0;ii<fNParticleFunctions;ii++)
700 fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
702 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
705 partEvent1 = partbuffer.Push(partEvent1);
706 }//while (fReader->Next() == kFALSE)
708 /*************************************************************************************/
710 void AliHBTAnalysis::WriteFunctions()
712 //Calls Write for all defined functions in analysis
713 //== writes all results
715 for(ii = 0;ii<fNParticleFunctions;ii++)
716 fParticleFunctions[ii]->Write();
718 for(ii = 0;ii<fNTrackFunctions;ii++)
719 fTrackFunctions[ii]->Write();
721 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
722 fParticleAndTrackFunctions[ii]->Write();
724 for(ii = 0;ii<fNParticleMonitorFunctions;ii++)
725 fParticleMonitorFunctions[ii]->Write();
727 for(ii = 0;ii<fNTrackMonitorFunctions;ii++)
728 fTrackMonitorFunctions[ii]->Write();
730 for(ii = 0;ii<fNParticleAndTrackMonitorFunctions;ii++)
731 fParticleAndTrackMonitorFunctions[ii]->Write();
733 /*************************************************************************************/
735 void AliHBTAnalysis::SetGlobalPairCut(AliHBTPairCut* cut)
737 //Sets the global cut
740 Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
743 fPairCut = (AliHBTPairCut*)cut->Clone();
746 /*************************************************************************************/
748 void AliHBTAnalysis::AddTrackFunction(AliHBTOnePairFctn* f)
750 //Adds track function
751 if (f == 0x0) return;
752 if (fNTrackFunctions == fgkFctnArraySize)
754 Error("AliHBTAnalysis::AddTrackFunction","Can not add this function, not enough place in the array.");
756 fTrackFunctions[fNTrackFunctions] = f;
759 /*************************************************************************************/
761 void AliHBTAnalysis::AddParticleFunction(AliHBTOnePairFctn* f)
763 //adds particle function
764 if (f == 0x0) return;
766 if (fNParticleFunctions == fgkFctnArraySize)
768 Error("AliHBTAnalysis::AddParticleFunction","Can not add this function, not enough place in the array.");
770 fParticleFunctions[fNParticleFunctions] = f;
771 fNParticleFunctions++;
773 /*************************************************************************************/
775 void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTTwoPairFctn* f)
777 //add resolution function
778 if (f == 0x0) return;
779 if (fNParticleAndTrackFunctions == fgkFctnArraySize)
781 Error("AliHBTAnalysis::AddParticleAndTrackFunction","Can not add this function, not enough place in the array.");
783 fParticleAndTrackFunctions[fNParticleAndTrackFunctions] = f;
784 fNParticleAndTrackFunctions++;
786 /*************************************************************************************/
788 void AliHBTAnalysis::AddParticleMonitorFunction(AliHBTMonOneParticleFctn* f)
790 //add particle monitoring function
791 if (f == 0x0) return;
793 if (fNParticleMonitorFunctions == fgkFctnArraySize)
795 Error("AliHBTAnalysis::AddParticleMonitorFunction","Can not add this function, not enough place in the array.");
797 fParticleMonitorFunctions[fNParticleMonitorFunctions] = f;
798 fNParticleMonitorFunctions++;
800 /*************************************************************************************/
802 void AliHBTAnalysis::AddTrackMonitorFunction(AliHBTMonOneParticleFctn* f)
804 //add track monitoring function
805 if (f == 0x0) return;
807 if (fNTrackMonitorFunctions == fgkFctnArraySize)
809 Error("AliHBTAnalysis::AddTrackMonitorFunction","Can not add this function, not enough place in the array.");
811 fTrackMonitorFunctions[fNTrackMonitorFunctions] = f;
812 fNTrackMonitorFunctions++;
814 /*************************************************************************************/
816 void AliHBTAnalysis::AddParticleAndTrackMonitorFunction(AliHBTMonTwoParticleFctn* f)
818 //add resolution monitoring function
819 if (f == 0x0) return;
820 if (fNParticleAndTrackMonitorFunctions == fgkFctnArraySize)
822 Error("AliHBTAnalysis::AddParticleAndTrackMonitorFunction","Can not add this function, not enough place in the array.");
824 fParticleAndTrackMonitorFunctions[fNParticleAndTrackMonitorFunctions] = f;
825 fNParticleAndTrackMonitorFunctions++;
829 /*************************************************************************************/
830 /*************************************************************************************/
832 Bool_t AliHBTAnalysis::RunCoherencyCheck()
834 //Checks if both HBTRuns are similar
835 //return true if error found
836 //if they seem to be OK return false
838 Info("RunCoherencyCheck","Checking HBT Runs Coherency");
840 Info("RunCoherencyCheck","Number of events ...");
841 if (fReader->GetNumberOfPartEvents() == fReader->GetNumberOfTrackEvents() ) //check whether there is the same number of events
843 Info("RunCoherencyCheck","OK. %d found\n",fReader->GetNumberOfTrackEvents());
846 { //if not the same - ERROR
847 Error("AliHBTAnalysis::RunCoherencyCheck()",
848 "Number of simulated events (%d) is not equal to number of reconstructed events(%d)",
849 fReader->GetNumberOfPartEvents(),fReader->GetNumberOfTrackEvents());
853 Info("RunCoherencyCheck","Checking number of Particles AND Particles Types in each event ...");
855 AliHBTEvent *partEvent;
856 AliHBTEvent *trackEvent;
857 for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
859 partEvent= fReader->GetParticleEvent(i); //gets the "ith" event
860 trackEvent = fReader->GetTrackEvent(i);
862 if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
863 if ( (partEvent == 0x0) || (partEvent == 0x0) )
865 Error("AliHBTAnalysis::RunCoherencyCheck()",
866 "One event is NULL and the other one not. Event Number %d",i);
870 if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
872 Error("AliHBTAnalysis::RunCoherencyCheck()",
873 "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
874 i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
878 for (Int_t j = 0; j<partEvent->GetNumberOfParticles(); j++)
880 if( partEvent->GetParticle(j)->GetPdgCode() != trackEvent->GetParticle(j)->GetPdgCode() )
882 Error("AliHBTAnalysis::RunCoherencyCheck()",
883 "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
884 i,j, partEvent->GetParticle(j)->GetPdgCode(),trackEvent->GetParticle(j)->GetPdgCode() );
890 Info("RunCoherencyCheck"," Done");
891 Info("RunCoherencyCheck"," Everything looks OK");
895 /*************************************************************************************/
897 void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
899 //Performs analysis for both, tracks and particles
901 AliHBTParticle * part1, * part2;
902 AliHBTParticle * track1, * track2;
904 AliHBTEvent * trackEvent1=0x0,*partEvent1=0x0;
905 AliHBTEvent * trackEvent2=0x0,*partEvent2=0x0;
906 AliHBTEvent * trackEvent3=0x0,*partEvent3=0x0;
908 AliHBTEvent * rawtrackEvent, *rawpartEvent;//this we get from Reader
910 AliHBTPair * trackpair = new AliHBTPair();
911 AliHBTPair * partpair = new AliHBTPair();
913 AliHBTEventBuffer partbuffer(fBufferSize);
914 AliHBTEventBuffer trackbuffer(fBufferSize);
918 trackEvent1 = new AliHBTEvent();
919 partEvent1 = new AliHBTEvent();
920 trackEvent1->SetOwner(kFALSE);
921 partEvent1->SetOwner(kFALSE);;
923 Bool_t nocorrfctns = (fNParticleFunctions == 0) && (fNTrackFunctions == 0) && (fNParticleAndTrackFunctions == 0);
927 Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
928 Info("ProcessTracksAndParticlesNonIdentAnal","***** NON IDENT MODE ****************");
929 Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
931 for (Int_t i = 0;;i++)//infinite loop
933 if (fReader->Next()) break; //end when no more events available
935 rawpartEvent = fReader->GetParticleEvent();
936 rawtrackEvent = fReader->GetTrackEvent();
937 if ( (rawpartEvent == 0x0) || (rawtrackEvent == 0x0) ) continue;//in case of any error
939 if ( rawpartEvent->GetNumberOfParticles() != rawtrackEvent->GetNumberOfParticles() )
941 Fatal("ProcessTracksAndParticlesNonIdentAnal",
942 "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
943 i,rawpartEvent->GetNumberOfParticles() , rawtrackEvent->GetNumberOfParticles());
947 /********************************/
949 /********************************/
950 if ( ( (partEvent2==0x0) || (trackEvent2==0x0)) )//in case fBufferSize == 0 and pointers are created do not eneter
952 partEvent2 = new AliHBTEvent();
953 trackEvent2 = new AliHBTEvent();
954 partEvent2->SetOwner(kTRUE);
955 trackEvent2->SetOwner(kTRUE);
958 FilterOut(partEvent1, partEvent2, rawpartEvent, trackEvent1, trackEvent2, rawtrackEvent);
960 for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
962 if ( (j%fDisplayMixingInfo) == 0)
963 Info("ProcessTracksAndParticlesNonIdentAnal",
964 "Mixing particle %d from event %d with particles from event %d",j,i,i);
966 part1= partEvent1->GetParticle(j);
967 track1= trackEvent1->GetParticle(j);
969 //PID reconstruction imperfections
970 // if( part1->GetPdgCode() != track1->GetPdgCode() )
972 // Fatal("ProcessTracksAndParticlesNonIdentAnal",
973 // "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
974 // i,j, part1->GetPdgCode(),track1->GetPdgCode() );
978 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
979 fParticleMonitorFunctions[ii]->Process(part1);
980 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
981 fTrackMonitorFunctions[ii]->Process(track1);
982 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
983 fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
985 if (nocorrfctns) continue;
987 /***************************************/
988 /****** filling numerators ********/
989 /****** (particles from event2) ********/
990 /***************************************/
992 for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups
994 part2= partEvent2->GetParticle(k);
995 if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
996 partpair->SetParticles(part1,part2);
998 track2= trackEvent2->GetParticle(k);
999 trackpair->SetParticles(track1,track2);
1001 if( (fPairCut->PassPairProp(trackpair)) ) //check pair cut
1002 { //do not meets crietria of the pair cut
1006 {//meets criteria of the pair cut
1007 for(ii = 0;ii<fNParticleFunctions;ii++)
1008 fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
1010 for(ii = 0;ii<fNTrackFunctions;ii++)
1011 fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
1013 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
1014 fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(trackpair,partpair);
1018 if ( fBufferSize == 0) continue;//do not mix diff histograms
1019 /***************************************/
1020 /***** Filling denominators *********/
1021 /***************************************/
1022 partbuffer.ResetIter();
1023 trackbuffer.ResetIter();
1027 while ( (partEvent3 = partbuffer.Next() ) != 0x0)
1029 trackEvent3 = trackbuffer.Next();
1031 if ( (j%fDisplayMixingInfo) == 0)
1032 Info("ProcessTracksAndParticlesNonIdentAnal",
1033 "Mixing particle %d from event %d with particles from event %d",j,i,i-(++nmonitor));
1035 for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
1037 part2= partEvent3->GetParticle(k);
1038 partpair->SetParticles(part1,part2);
1040 track2= trackEvent3->GetParticle(k);
1041 trackpair->SetParticles(track1,track2);
1043 if( (fPairCut->PassPairProp(trackpair)) ) //check pair cut
1044 { //do not meets crietria of the pair cut
1048 {//meets criteria of the pair cut
1050 for(ii = 0;ii<fNParticleFunctions;ii++)
1051 fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
1053 for(ii = 0;ii<fNTrackFunctions;ii++)
1054 fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
1056 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
1057 fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(trackpair,partpair);
1059 }// for particles event2
1061 }//for over particles in event1
1063 partEvent2 = partbuffer.Push(partEvent2);
1064 trackEvent2 = trackbuffer.Push(trackEvent2);
1066 }//end of loop over events (1)
1068 partbuffer.SetOwner(kTRUE);
1069 trackbuffer.SetOwner(kTRUE);
1079 /*************************************************************************************/
1081 void AliHBTAnalysis::ProcessTracksNonIdentAnal()
1083 //Process Tracks only with non identical mode
1084 AliHBTParticle * track1, * track2;
1086 AliHBTEvent * trackEvent1=0x0;
1087 AliHBTEvent * trackEvent2=0x0;
1088 AliHBTEvent * trackEvent3=0x0;
1090 AliHBTEvent * rawtrackEvent;
1092 AliHBTPair * trackpair = new AliHBTPair();
1093 AliHBTEventBuffer trackbuffer(fBufferSize);
1097 trackEvent1 = new AliHBTEvent();
1098 trackEvent1->SetOwner(kFALSE);
1102 Info("ProcessTracksNonIdentAnal","**************************************");
1103 Info("ProcessTracksNonIdentAnal","***** NON IDENT MODE ****************");
1104 Info("ProcessTracksNonIdentAnal","**************************************");
1107 for (Int_t i = 0;;i++)//infinite loop
1109 if (fReader->Next()) break; //end when no more events available
1110 rawtrackEvent = fReader->GetTrackEvent();
1112 if (rawtrackEvent == 0x0) continue;//in case of any error
1114 /********************************/
1116 /********************************/
1117 if ( trackEvent2==0x0 )//in case fBufferSize == 0 and pointers are created do not eneter
1119 trackEvent2 = new AliHBTEvent();
1120 trackEvent2->SetOwner(kTRUE);
1123 FilterOut(trackEvent1, trackEvent2, rawtrackEvent);
1125 for (Int_t j = 0; j<trackEvent1->GetNumberOfParticles() ; j++)
1127 if ( (j%fDisplayMixingInfo) == 0)
1128 Info("ProcessTracksNonIdentAnal",
1129 "Mixing particle %d from event %d with particles from event %d",j,i,i);
1131 track1= trackEvent1->GetParticle(j);
1133 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
1134 fTrackMonitorFunctions[ii]->Process(track1);
1136 if (fNTrackFunctions == 0x0) continue;
1138 /***************************************/
1139 /****** filling numerators ********/
1140 /****** (particles from event2) ********/
1141 /***************************************/
1142 for (Int_t k = 0; k < trackEvent2->GetNumberOfParticles() ; k++)
1144 track2= trackEvent2->GetParticle(k);
1145 if (track1->GetUID() == track2->GetUID()) continue;//this is the same particle but with different PID
1146 trackpair->SetParticles(track1,track2);
1149 if( fPairCut->PassPairProp(trackpair)) //check pair cut
1150 { //do not meets crietria of the pair cut
1154 {//meets criteria of the pair cut
1156 for(ii = 0;ii<fNTrackFunctions;ii++)
1157 fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
1160 if ( fBufferSize == 0) continue;//do not mix diff histograms
1161 /***************************************/
1162 /***** Filling denominators *********/
1163 /***************************************/
1164 trackbuffer.ResetIter();
1167 while ( (trackEvent3 = trackbuffer.Next() ) != 0x0)
1170 if ( (j%fDisplayMixingInfo) == 0)
1171 Info("ProcessTracksNonIdentAnal",
1172 "Mixing particle %d from event %d with particles from event %d",j,i,i-(++nmonitor));
1174 for (Int_t k = 0; k < trackEvent3->GetNumberOfParticles() ; k++)
1177 track2= trackEvent3->GetParticle(k);
1178 trackpair->SetParticles(track1,track2);
1181 if( fPairCut->PassPairProp(trackpair)) //check pair cut
1182 { //do not meets crietria of the pair cut
1186 {//meets criteria of the pair cut
1187 for(ii = 0;ii<fNTrackFunctions;ii++)
1188 fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
1190 }// for particles event2
1192 }//for over particles in event1
1194 trackEvent2 = trackbuffer.Push(trackEvent2);
1196 }//end of loop over events (1)
1198 trackbuffer.SetOwner(kTRUE);
1203 /*************************************************************************************/
1205 void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
1207 //process paricles only with non identical mode
1208 AliHBTParticle * part1 = 0x0, * part2 = 0x0;
1210 AliHBTEvent * partEvent1 = 0x0;
1211 AliHBTEvent * partEvent2 = 0x0;
1212 AliHBTEvent * partEvent3 = 0x0;
1214 AliHBTEvent * rawpartEvent = 0x0;
1216 AliHBTPair * partpair = new AliHBTPair();
1217 AliHBTEventBuffer partbuffer(fBufferSize);
1221 partEvent1 = new AliHBTEvent();
1222 partEvent1->SetOwner(kFALSE);
1226 Info("ProcessParticlesNonIdentAnal","**************************************");
1227 Info("ProcessParticlesNonIdentAnal","***** NON IDENT MODE ****************");
1228 Info("ProcessParticlesNonIdentAnal","**************************************");
1230 for (Int_t i = 0;;i++)//infinite loop
1232 if (fReader->Next()) break; //end when no more events available
1234 rawpartEvent = fReader->GetParticleEvent();
1235 if ( rawpartEvent == 0x0 ) continue;//in case of any error
1237 /********************************/
1239 /********************************/
1240 if (partEvent2==0x0)//in case fBufferSize == 0 and pointers are created do not eneter
1242 partEvent2 = new AliHBTEvent();
1243 partEvent2->SetOwner(kTRUE);
1246 FilterOut(partEvent1, partEvent2, rawpartEvent);
1248 for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
1250 if ( (j%fDisplayMixingInfo) == 0)
1251 Info("ProcessParticlesNonIdentAnal",
1252 "Mixing particle %d from event %d with particles from event %d",j,i,i);
1254 part1= partEvent1->GetParticle(j);
1256 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
1257 fParticleMonitorFunctions[ii]->Process(part1);
1259 if (fNParticleFunctions == 0) continue;
1261 /***************************************/
1262 /****** filling numerators ********/
1263 /****** (particles from event2) ********/
1264 /***************************************/
1265 for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++)
1267 part2= partEvent2->GetParticle(k);
1268 if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
1269 partpair->SetParticles(part1,part2);
1271 if(fPairCut->PassPairProp(partpair) ) //check pair cut
1272 { //do not meets crietria of the pair cut
1276 {//meets criteria of the pair cut
1277 for(ii = 0;ii<fNParticleFunctions;ii++)
1278 fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
1281 if ( fBufferSize == 0) continue;//do not mix diff histograms
1282 /***************************************/
1283 /***** Filling denominators *********/
1284 /***************************************/
1285 partbuffer.ResetIter();
1288 while ( (partEvent3 = partbuffer.Next() ) != 0x0)
1290 if ( (j%fDisplayMixingInfo) == 0)
1291 Info("ProcessParticlesNonIdentAnal",
1292 "Mixing particle %d from event %d with particles from event %d",j,i,i-(++nmonitor));
1294 for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
1296 part2= partEvent3->GetParticle(k);
1297 partpair->SetParticles(part1,part2);
1299 if(fPairCut->PassPairProp(partpair) ) //check pair cut
1300 { //do not meets crietria of the pair cut
1304 {//meets criteria of the pair cut
1305 for(ii = 0;ii<fNParticleFunctions;ii++)
1307 fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
1310 }// for particles event2
1312 }//for over particles in event1
1313 partEvent2 = partbuffer.Push(partEvent2);
1314 }//end of loop over events (1)
1316 partbuffer.SetOwner(kTRUE);
1322 /*************************************************************************************/
1323 void AliHBTAnalysis::FilterOut(AliHBTEvent* outpart1, AliHBTEvent* outpart2, AliHBTEvent* inpart,
1324 AliHBTEvent* outtrack1, AliHBTEvent* outtrack2, AliHBTEvent* intrack)
1326 //Puts particles accepted as a first particle by global cut in out1
1327 //and as a second particle in out2
1329 AliHBTParticle* part, *track;
1336 AliHBTParticleCut *cut1 = fPairCut->GetFirstPartCut();
1337 AliHBTParticleCut *cut2 = fPairCut->GetSecondPartCut();
1341 for (Int_t i = 0; i < inpart->GetNumberOfParticles(); i++)
1344 part = inpart->GetParticle(i);
1345 track = intrack->GetParticle(i);
1347 if ( (cut1->Pass(track)) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
1348 if ( (cut2->Pass(track)) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
1350 if (gDebug)//to be removed in real analysis
1351 if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1353 //Particle accpted by both cuts
1354 Error("FilterOut","Particle accepted by both cuts");
1360 outpart1->AddParticle(part);
1361 outtrack1->AddParticle(track);
1367 outpart2->AddParticle(new AliHBTParticle(*part));
1368 outtrack2->AddParticle(new AliHBTParticle(*track));
1373 /*************************************************************************************/
1375 void AliHBTAnalysis::FilterOut(AliHBTEvent* out1, AliHBTEvent* out2, AliHBTEvent* in)
1377 //Puts particles accepted as a first particle by global cut in out1
1378 //and as a second particle in out2
1379 AliHBTParticle* part;
1384 AliHBTParticleCut *cut1 = fPairCut->GetFirstPartCut();
1385 AliHBTParticleCut *cut2 = fPairCut->GetSecondPartCut();
1389 for (Int_t i = 0; i < in->GetNumberOfParticles(); i++)
1392 part = in->GetParticle(i);
1394 if ( cut1->Pass(part) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
1395 if ( cut2->Pass(part) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
1397 if (gDebug)//to be removed in real analysis
1398 if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1400 //Particle accpted by both cuts
1401 Error("FilterOut","Particle accepted by both cuts");
1407 out1->AddParticle(part);
1413 out2->AddParticle(part);
1418 /*************************************************************************************/
1420 Bool_t AliHBTAnalysis::IsNonIdentAnalysis()
1422 //checks if it is possible to use special analysis for non identical particles
1423 //it means - in global pair cut first particle id is different than second one
1424 //and both are different from 0
1425 //in the future is possible to perform more sophisticated check
1426 //if cuts have excluding requirements
1428 if (fPairCut->IsEmpty())
1431 if (fPairCut->GetFirstPartCut()->IsEmpty())
1434 if (fPairCut->GetSecondPartCut()->IsEmpty())
1437 Int_t id1 = fPairCut->GetFirstPartCut()->GetPID();
1438 Int_t id2 = fPairCut->GetSecondPartCut()->GetPID();
1440 if ( (id1==0) || (id2==0) || (id1==id2) )