2 #include "AliHBTAnalysis.h"
7 #include "AliHBTEvent.h"
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 #include "AliHBTMonitorFunction.h"
16 #include <TBenchmark.h>
21 ClassImp(AliHBTAnalysis)
23 const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100;
24 const Int_t AliHBTAnalysis::fgkHbtAnalyzeAll = 0;
26 AliHBTAnalysis::AliHBTAnalysis()
30 fTrackFunctions = new AliHBTOnePairFctn* [fgkFctnArraySize];
31 fParticleFunctions = new AliHBTOnePairFctn* [fgkFctnArraySize];
32 fParticleAndTrackFunctions = new AliHBTTwoPairFctn* [fgkFctnArraySize];
34 fParticleMonitorFunctions = new AliHBTMonOneParticleFctn* [fgkFctnArraySize];
35 fTrackMonitorFunctions = new AliHBTMonOneParticleFctn* [fgkFctnArraySize];
36 fParticleAndTrackMonitorFunctions = new AliHBTMonTwoParticleFctn* [fgkFctnArraySize];
39 fNParticleFunctions = 0;
40 fNParticleAndTrackFunctions = 0;
42 fNParticleMonitorFunctions = 0;
43 fNTrackMonitorFunctions = 0;
44 fNParticleAndTrackMonitorFunctions = 0;
46 fPairCut = new AliHBTEmptyPairCut();//empty cut - accepts all particles
50 /*************************************************************************************/
52 AliHBTAnalysis::~AliHBTAnalysis()
55 //note that we do not delete functions itself
56 // they should be deleted by whom where created
57 //we only store pointers, and we use "new" only for pointers array
58 delete [] fTrackFunctions;
59 delete [] fParticleFunctions;
60 delete [] fParticleAndTrackFunctions;
62 delete [] fParticleMonitorFunctions;
63 delete [] fTrackMonitorFunctions;
64 delete [] fParticleAndTrackMonitorFunctions;
66 delete fPairCut; // always have an copy of an object - we create we dstroy
69 /*************************************************************************************/
71 void AliHBTAnalysis::Process(Option_t* option)
73 //default option = "TracksAndParticles"
74 //Main method of the HBT Analysis Package
75 //It triggers reading with the global cut (default is an empty cut)
76 //Than it checks options and data which are read
77 //if everything is OK, then it calls one of the looping methods
78 //depending on tfReaderhe option
79 //These methods differs on what thay are looping on
82 //--------------------------------------------------------------------
83 //ProcessTracksAndParticles - "TracksAndParticles"
85 // it loops over both, tracks(reconstructed) and particles(simulated)
86 // all function gethered in all 3 lists are called for each (double)pair
88 //ProcessTracks - "Tracks"
89 // it loops only on tracks(reconstructed),
90 // functions ONLY from fTrackFunctions list are called
92 //ProcessParticles - "Particles"
93 // it loops only on particles(simulated),
94 // functions ONLY from fParticleAndTrackFunctions list are called
99 Error("Process","The reader is not set");
104 const char *oT = strstr(option,"Tracks");
105 const char *oP = strstr(option,"Particles");
107 Bool_t nonid = IsNonIdentAnalysis();
109 cout<<"nonid = "<<nonid<<endl;
113 if (fReader->GetNumberOfPartEvents() <1)
115 Error("Process","There is no Particles. Maybe change the option?");
118 if (fReader->GetNumberOfTrackEvents() <1)
120 Error("Process","There is no Tracks. Maybe change the option?");
124 if ( RunCoherencyCheck() )
127 "Coherency check not passed. Maybe change the option?\n");
130 if (nonid) ProcessTracksAndParticlesNonIdentAnal();
131 else ProcessTracksAndParticles();
137 if (fReader->GetNumberOfTrackEvents() <1)
139 Error("Process","There is no data to analyze.");
142 if (nonid) ProcessTracksNonIdentAnal();
143 else ProcessTracks();
149 if (fReader->GetNumberOfPartEvents() <1)
151 Error("Process","There is no data to analyze.");
154 if (nonid) ProcessParticlesNonIdentAnal();
155 else ProcessParticles();
157 // cout<<"NON ID"<<endl;
158 // ProcessParticlesNonIdentAnal();
159 // cout<<"\n\n\n NORMAL"<<endl;
160 // ProcessParticles();
167 /*************************************************************************************/
169 void AliHBTAnalysis::ProcessTracksAndParticles()
172 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
173 //the loops are splited
176 AliHBTParticle * part1, * part2;
177 AliHBTParticle * track1, * track2;
179 AliHBTEvent * trackEvent, *partEvent;
180 AliHBTEvent * trackEvent2,*partEvent2;
182 // Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
184 Int_t Nev = fReader->GetNumberOfTrackEvents();
186 /***************************************/
187 /****** Looping same events ********/
188 /****** filling numerators ********/
189 /***************************************/
190 AliHBTPair * trackpair = new AliHBTPair();
191 AliHBTPair * partpair = new AliHBTPair();
193 AliHBTPair * tmptrackpair;//temprary pointers to pairs
194 AliHBTPair * tmppartpair;
198 for (Int_t i = 0;i<Nev;i++)
200 partEvent= fReader->GetParticleEvent(i);
201 trackEvent = fReader->GetTrackEvent(i);
203 if (!partEvent) continue;
207 for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
209 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
211 part1= partEvent->GetParticle(j);
212 track1= trackEvent->GetParticle(j);
215 for(zz = 0; zz<fNParticleMonitorFunctions; zz++)
216 fParticleMonitorFunctions[zz]->ProcessSameEventParticles(part1);
217 for(zz = 0; zz<fNTrackMonitorFunctions; zz++)
218 fTrackMonitorFunctions[zz]->ProcessSameEventParticles(track1);
219 for(zz = 0; zz<fNParticleAndTrackMonitorFunctions; zz++)
220 fParticleAndTrackMonitorFunctions[zz]->ProcessSameEventParticles(track1,part1);
222 if ( (fNParticleFunctions == 0) && (fNTrackFunctions ==0) && (fNParticleAndTrackFunctions == 0))
225 for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
227 part2= partEvent->GetParticle(k);
228 partpair->SetParticles(part1,part2);
230 track2= trackEvent->GetParticle(k);
231 trackpair->SetParticles(track1,track2);
233 if(fPairCut->Pass(partpair) || (fPairCut->Pass(trackpair))) //check pair cut
234 { //do not meets crietria of the pair cut, try with swaped pairs
235 if( ( fPairCut->Pass(partpair->GetSwapedPair()) ) || ( fPairCut->Pass(trackpair->GetSwapedPair()) ) )
236 continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
238 { //swaped pair meets all the criteria
239 tmppartpair = partpair->GetSwapedPair();
240 tmptrackpair = trackpair->GetSwapedPair();
244 {//meets criteria of the pair cut
245 tmptrackpair = trackpair;
246 tmppartpair = partpair;
249 for(ii = 0;ii<fNParticleFunctions;ii++)
250 fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
252 for(ii = 0;ii<fNTrackFunctions;ii++)
253 fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
255 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
256 fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair,tmppartpair);
262 /***** Filling denominators *********/
263 /***************************************/
264 for (Int_t i = 0;i<Nev-1;i++) //In each event (but last) ....
267 if ((fNParticleFunctions == 0) && (fNTrackFunctions ==0) && (fNParticleAndTrackFunctions == 0))
270 partEvent= fReader->GetParticleEvent(i);
271 if (!partEvent) continue;
273 trackEvent = fReader->GetTrackEvent(i);
277 for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
280 part1= partEvent->GetParticle(j);
282 track1= trackEvent->GetParticle(j);
286 if ( ((i+fBufferSize) >= Nev) ||( fBufferSize < 0) ) //if buffer size is negative
287 //or current event+buffersize is greater
288 //than max nuber of events
290 NNN = Nev; //set the max event number
294 NNN = i+fBufferSize; //set the current event number + fBufferSize
297 for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
300 partEvent2= fReader->GetParticleEvent(k);
301 if (!partEvent2) continue;
303 trackEvent2 = fReader->GetTrackEvent(k);
305 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
307 for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
310 // if (N>MAXCOMB) break;
312 part2= partEvent2->GetParticle(l);
313 partpair->SetParticles(part1,part2);
315 track2= trackEvent2->GetParticle(l);
316 trackpair->SetParticles(track1,track2);
318 if( (fPairCut->Pass(partpair)) || (fPairCut->Pass(trackpair)) ) //check pair cut
319 { //do not meets crietria of the
320 if( ( fPairCut->Pass(partpair->GetSwapedPair()) ) || ( fPairCut->Pass(trackpair->GetSwapedPair()) ) )
324 tmppartpair = partpair->GetSwapedPair();
325 tmptrackpair = trackpair->GetSwapedPair();
329 {//meets criteria of the pair cut
330 tmptrackpair = trackpair;
331 tmppartpair = partpair;
334 for(ii = 0;ii<fNParticleFunctions;ii++)
335 fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
337 for(ii = 0;ii<fNTrackFunctions;ii++)
338 fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
340 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
341 fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair);
344 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
345 }//for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
350 /***************************************/
354 /*************************************************************************************/
356 void AliHBTAnalysis::ProcessTracks()
358 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
359 //the loops are splited
360 AliHBTParticle * track1, * track2;
362 AliHBTEvent * trackEvent;
363 AliHBTEvent * trackEvent2;
365 // Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
367 Int_t Nev = fReader->GetNumberOfTrackEvents();
369 /***************************************/
370 /****** Looping same events ********/
371 /****** filling numerators ********/
372 /***************************************/
373 AliHBTPair * trackpair = new AliHBTPair();
374 AliHBTPair * tmptrackpair; //temporary pointer
376 for (Int_t i = 0;i<Nev;i++)
378 trackEvent = fReader->GetTrackEvent(i);
379 if (!trackEvent) continue;
382 for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
384 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
386 track1= trackEvent->GetParticle(j);
389 for(zz = 0; zz<fNTrackMonitorFunctions; zz++)
390 fTrackMonitorFunctions[zz]->ProcessSameEventParticles(track1);
392 if ( fNTrackFunctions ==0 )
395 for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
397 track2= trackEvent->GetParticle(k);
398 trackpair->SetParticles(track1,track2);
399 if(fPairCut->Pass(trackpair)) //check pair cut
400 { //do not meets crietria of the
401 if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
402 else tmptrackpair = trackpair->GetSwapedPair();
406 tmptrackpair = trackpair;
410 for(ii = 0;ii<fNTrackFunctions;ii++)
411 fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
418 /***************************************/
419 /***** Filling diff histogram *********/
420 /***************************************/
421 for (Int_t i = 0;i<Nev-1;i++) //In each event (but last) ....
423 if ( fNTrackFunctions ==0 )
426 trackEvent = fReader->GetTrackEvent(i);
427 if (!trackEvent) continue;
430 for (Int_t j = 0; j< trackEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
432 // if (N>MAXCOMB) break;
434 track1= trackEvent->GetParticle(j);
438 if ( ((i+fBufferSize) >= Nev) ||( fBufferSize < 0) ) //if buffer size is negative
439 //or current event+buffersize is greater
440 //than max nuber of events
442 NNN = Nev; //set the max event number
446 NNN = i+fBufferSize; //set the current event number + fBufferSize
449 for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
452 trackEvent2 = fReader->GetTrackEvent(k);
453 if (!trackEvent2) continue;
455 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
457 for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++) // ... on all particles
460 // if (N>MAXCOMB) break;
461 track2= trackEvent2->GetParticle(l);
462 trackpair->SetParticles(track1,track2);
464 if(fPairCut->Pass(trackpair)) //check pair cut
465 { //do not meets crietria of the
466 if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
467 else tmptrackpair = trackpair->GetSwapedPair();
471 tmptrackpair = trackpair;
474 for(ii = 0;ii<fNTrackFunctions;ii++)
475 fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
477 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
478 }//for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
483 /***************************************/
488 /*************************************************************************************/
490 void AliHBTAnalysis::ProcessParticles()
492 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
493 //the loops are splited
494 AliHBTParticle * part1, * part2;
496 AliHBTEvent * partEvent;
497 AliHBTEvent * partEvent2;
499 AliHBTPair * partpair = new AliHBTPair();
500 AliHBTPair * tmppartpair; //temporary pointer to the pair
502 // Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
504 Int_t Nev = fReader->GetNumberOfPartEvents();
507 /***************************************/
508 /****** Looping same events ********/
509 /****** filling numerators ********/
510 /***************************************/
511 // gBenchmark->Start("id");
513 for (Int_t i = 0;i<Nev;i++)
515 partEvent= fReader->GetParticleEvent(i);
516 if (!partEvent) continue;
519 for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
521 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
523 part1= partEvent->GetParticle(j);
526 for(zz = 0; zz<fNParticleMonitorFunctions; zz++)
527 fParticleMonitorFunctions[zz]->ProcessSameEventParticles(part1);
529 if ( fNParticleFunctions ==0 )
532 for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
534 part2= partEvent->GetParticle(k);
535 partpair->SetParticles(part1,part2);
537 if( fPairCut->Pass(partpair) ) //check pair cut
538 { //do not meets crietria of the pair cut, try with swaped pairs
539 if( fPairCut->Pass(partpair->GetSwapedPair() ) )
540 continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
542 { //swaped pair meets all the criteria
543 tmppartpair = partpair->GetSwapedPair();
548 tmppartpair = partpair;
552 for(ii = 0;ii<fNParticleFunctions;ii++)
553 fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
558 /***************************************/
559 /***** Filling diff histogram *********/
560 /***************************************/
561 for (Int_t i = 0;i<Nev-1;i++) //In each event (but last)....
563 if ( fNParticleFunctions ==0 )
566 partEvent= fReader->GetParticleEvent(i);
567 if (!partEvent) continue;
571 for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
573 // if (N>MAXCOMB) break;
575 part1= partEvent->GetParticle(j);
579 if ( ((i+fBufferSize) >= Nev) ||( fBufferSize < 0) ) //if buffer size is negative
580 //or current event+buffersize is greater
581 //than max nuber of events
583 NNN = Nev; //set the max event number
587 NNN = i+fBufferSize; //set the current event number + fBufferSize
590 // cout<<"NNN = "<<NNN<<endl;
591 for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
594 partEvent2= fReader->GetParticleEvent(k);
595 if (!partEvent2) continue;
597 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
599 for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
602 // if (N>MAXCOMB) break;
603 part2= partEvent2->GetParticle(l);
604 partpair->SetParticles(part1,part2);
606 if(fPairCut->Pass(partpair)) //check pair cut
607 { //do not meets crietria of the
608 if( fPairCut->Pass(partpair->GetSwapedPair()) ) continue;
609 else tmppartpair = partpair->GetSwapedPair();
613 tmppartpair = partpair;
616 for(ii = 0;ii<fNParticleFunctions;ii++)
617 fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
619 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
620 }//for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
625 /***************************************/
627 // gBenchmark->Stop("id");
628 // gBenchmark->Show("id");
632 /*************************************************************************************/
635 void AliHBTAnalysis::WriteFunctions()
638 for(ii = 0;ii<fNParticleFunctions;ii++)
639 fParticleFunctions[ii]->Write();
641 for(ii = 0;ii<fNTrackFunctions;ii++)
642 fTrackFunctions[ii]->Write();
644 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
645 fParticleAndTrackFunctions[ii]->Write();
647 for(ii = 0;ii<fNParticleMonitorFunctions;ii++)
648 fParticleMonitorFunctions[ii]->Write();
650 for(ii = 0;ii<fNTrackMonitorFunctions;ii++)
651 fTrackMonitorFunctions[ii]->Write();
653 for(ii = 0;ii<fNParticleAndTrackMonitorFunctions;ii++)
654 fParticleAndTrackMonitorFunctions[ii]->Write();
656 /*************************************************************************************/
658 void AliHBTAnalysis::SetGlobalPairCut(AliHBTPairCut* cut)
662 Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
665 fPairCut = (AliHBTPairCut*)cut->Clone();
668 /*************************************************************************************/
670 void AliHBTAnalysis::AddTrackFunction(AliHBTOnePairFctn* f)
672 if (f == 0x0) return;
673 if (fNTrackFunctions == fgkFctnArraySize)
675 Error("AliHBTAnalysis::AddTrackFunction","Can not add this function, not enough place in the array.");
677 fTrackFunctions[fNTrackFunctions] = f;
680 /*************************************************************************************/
682 void AliHBTAnalysis::AddParticleFunction(AliHBTOnePairFctn* f)
684 if (f == 0x0) return;
686 if (fNParticleFunctions == fgkFctnArraySize)
688 Error("AliHBTAnalysis::AddParticleFunction","Can not add this function, not enough place in the array.");
690 fParticleFunctions[fNParticleFunctions] = f;
691 fNParticleFunctions++;
693 /*************************************************************************************/
695 void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTTwoPairFctn* f)
697 if (f == 0x0) return;
698 if (fNParticleAndTrackFunctions == fgkFctnArraySize)
700 Error("AliHBTAnalysis::AddParticleAndTrackFunction","Can not add this function, not enough place in the array.");
702 fParticleAndTrackFunctions[fNParticleAndTrackFunctions] = f;
703 fNParticleAndTrackFunctions++;
705 /*************************************************************************************/
707 void AliHBTAnalysis::AddParticleMonitorFunction(AliHBTMonOneParticleFctn* f)
709 if (f == 0x0) return;
711 if (fNParticleMonitorFunctions == fgkFctnArraySize)
713 Error("AliHBTAnalysis::AddParticleMonitorFunction","Can not add this function, not enough place in the array.");
715 fParticleMonitorFunctions[fNParticleMonitorFunctions] = f;
716 fNParticleMonitorFunctions++;
718 /*************************************************************************************/
720 void AliHBTAnalysis::AddTrackMonitorFunction(AliHBTMonOneParticleFctn* f)
722 if (f == 0x0) return;
724 if (fNTrackMonitorFunctions == fgkFctnArraySize)
726 Error("AliHBTAnalysis::AddTrackMonitorFunction","Can not add this function, not enough place in the array.");
728 fTrackMonitorFunctions[fNTrackMonitorFunctions] = f;
729 fNTrackMonitorFunctions++;
731 /*************************************************************************************/
733 void AliHBTAnalysis::AddParticleAndTrackMonitorFunction(AliHBTMonTwoParticleFctn* f)
735 if (f == 0x0) return;
736 if (fNParticleAndTrackMonitorFunctions == fgkFctnArraySize)
738 Error("AliHBTAnalysis::AddParticleAndTrackMonitorFunction","Can not add this function, not enough place in the array.");
740 fParticleAndTrackMonitorFunctions[fNParticleAndTrackMonitorFunctions] = f;
741 fNParticleAndTrackMonitorFunctions++;
745 /*************************************************************************************/
746 /*************************************************************************************/
748 Bool_t AliHBTAnalysis::RunCoherencyCheck()
750 //Checks if both HBTRuns are similar
751 //return true if error found
752 //if they seem to be OK return false
754 cout<<"Checking HBT Runs Coherency"<<endl;
756 cout<<"Number of events ..."; fflush(stdout);
758 if (fReader->GetNumberOfPartEvents() == fReader->GetNumberOfTrackEvents() ) //check whether there is the same number of events
760 cout<<"OK. "<<fReader->GetNumberOfTrackEvents()<<" found"<<endl;
763 { //if not the same - ERROR
764 Error("AliHBTAnalysis::RunCoherencyCheck()",
765 "Number of simulated events (%d) is not equal to number of reconstructed events(%d)",
766 fReader->GetNumberOfPartEvents(),fReader->GetNumberOfTrackEvents());
770 cout<<"Checking number of Particles AND Particles Types in each event ...";fflush(stdout);
772 AliHBTEvent *partEvent;
773 AliHBTEvent *trackEvent;
774 for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
776 partEvent= fReader->GetParticleEvent(i); //gets the "ith" event
777 trackEvent = fReader->GetTrackEvent(i);
779 if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
780 if ( (partEvent == 0x0) || (partEvent == 0x0) )
782 Error("AliHBTAnalysis::RunCoherencyCheck()",
783 "One event is NULL and the other one not. Event Number %d",i);
787 if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
789 Error("AliHBTAnalysis::RunCoherencyCheck()",
790 "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
791 i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
795 for (Int_t j = 0; j<partEvent->GetNumberOfParticles(); j++)
797 if( partEvent->GetParticle(j)->GetPdgCode() != trackEvent->GetParticle(j)->GetPdgCode() )
799 Error("AliHBTAnalysis::RunCoherencyCheck()",
800 "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
801 i,j, partEvent->GetParticle(j)->GetPdgCode(),trackEvent->GetParticle(j)->GetPdgCode() );
808 cout<<" Everything looks OK"<<endl;
812 /*************************************************************************************/
814 void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
816 AliHBTParticle * part1, * part2;
817 AliHBTParticle * track1, * track2;
819 AliHBTEvent * trackEvent1,*partEvent1;
820 AliHBTEvent * trackEvent2,*partEvent2;
821 AliHBTEvent * trackEvent3,*partEvent3;
823 AliHBTEvent * rawtrackEvent, *rawpartEvent;
824 // AliHBTEvent * rawtrackEvent2,*rawpartEvent2;
827 Int_t Nev = fReader->GetNumberOfTrackEvents();
829 AliHBTPair * trackpair = new AliHBTPair();
830 AliHBTPair * partpair = new AliHBTPair();
836 trackEvent1 = new AliHBTEvent();
837 partEvent1 = new AliHBTEvent();
838 trackEvent1->SetOwner(kFALSE);
839 partEvent1->SetOwner(kFALSE);;
841 cout<<"**************************************"<<endl;
842 cout<<"***** NON IDENT MODE ****************"<<endl;
843 cout<<"**************************************"<<endl;
844 for (Int_t i = 0;i<Nev;i++)
846 rawpartEvent = fReader->GetParticleEvent(i);
847 rawtrackEvent = fReader->GetTrackEvent(i);
848 if ( (rawpartEvent == 0x0) || (rawtrackEvent == 0x0) ) continue;//in case of any error
850 /********************************/
852 /********************************/
853 if ((ninbuffer > fBufferSize) && (fBufferSize > 0))
854 {//if we have in buffer desired number of events, use the last. If fBufferSize<0 mix all events for background
855 partEvent2 = (AliHBTEvent*)pbuffer.Remove(pbuffer.Last()); //remove from the end to be reset, filled and put on the beginning
856 trackEvent2 = (AliHBTEvent*)tbuffer.Remove(tbuffer.Last());
861 partEvent2 = new AliHBTEvent();
862 trackEvent2 = new AliHBTEvent();
863 partEvent2->SetOwner(kFALSE);
864 trackEvent2->SetOwner(kFALSE);
866 FilterOut(partEvent1, partEvent2, rawpartEvent, trackEvent1, trackEvent2, rawtrackEvent);
868 for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
870 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
872 part1= partEvent1->GetParticle(j);
873 track1= trackEvent1->GetParticle(j);
876 for(zz = 0; zz<fNParticleMonitorFunctions; zz++)
877 fParticleMonitorFunctions[zz]->ProcessSameEventParticles(part1);
878 for(zz = 0; zz<fNTrackMonitorFunctions; zz++)
879 fTrackMonitorFunctions[zz]->ProcessSameEventParticles(track1);
880 for(zz = 0; zz<fNParticleAndTrackMonitorFunctions; zz++)
881 fParticleAndTrackMonitorFunctions[zz]->ProcessSameEventParticles(track1,part1);
883 /***************************************/
884 /****** filling numerators ********/
885 /****** (particles from event2) ********/
886 /***************************************/
887 for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++)
889 part2= partEvent2->GetParticle(k);
890 partpair->SetParticles(part1,part2);
892 track2= trackEvent2->GetParticle(k);
893 trackpair->SetParticles(track1,track2);
896 if( (fPairCut->PassPairProp(partpair)) || (fPairCut->PassPairProp(trackpair))) //check pair cut
897 { //do not meets crietria of the pair cut
901 {//meets criteria of the pair cut
903 for(ii = 0;ii<fNParticleFunctions;ii++)
904 fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
906 for(ii = 0;ii<fNTrackFunctions;ii++)
907 fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
909 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
910 fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(trackpair,partpair);
915 /***************************************/
916 /***** Filling denominators *********/
917 /***************************************/
918 TIter piter(&pbuffer);
919 TIter titer(&tbuffer);
922 while ( (partEvent3 = (AliHBTEvent*)piter.Next()) != 0x0)
924 trackEvent3 = (AliHBTEvent*)titer.Next();
925 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i - (++nmonitor)<<endl;
927 for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
929 part2= partEvent3->GetParticle(k);
930 partpair->SetParticles(part1,part2);
932 track2= trackEvent3->GetParticle(k);
933 trackpair->SetParticles(track1,track2);
936 if( (fPairCut->PassPairProp(partpair)) || (fPairCut->PassPairProp(trackpair))) //check pair cut
937 { //do not meets crietria of the pair cut
941 {//meets criteria of the pair cut
943 for(ii = 0;ii<fNParticleFunctions;ii++)
944 fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
946 for(ii = 0;ii<fNTrackFunctions;ii++)
947 fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
949 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
950 fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(trackpair,partpair);
952 }// for particles event2
954 }//for over particles in event1
956 pbuffer.AddFirst(partEvent2);
957 tbuffer.AddFirst(trackEvent2);
960 }//end of loop over events (1)
962 pbuffer.SetOwner(); //to clean stored events by the way of buffer destruction
965 /*************************************************************************************/
967 void AliHBTAnalysis::ProcessTracksNonIdentAnal()
969 AliHBTParticle * track1, * track2;
971 AliHBTEvent * trackEvent1;
972 AliHBTEvent * trackEvent2;
973 AliHBTEvent * trackEvent3;
975 AliHBTEvent * rawtrackEvent;
976 // AliHBTEvent * rawtrackEvent2,*rawpartEvent2;
979 Int_t Nev = fReader->GetNumberOfTrackEvents();
981 AliHBTPair * trackpair = new AliHBTPair();
986 trackEvent1 = new AliHBTEvent();
987 trackEvent1->SetOwner(kFALSE);
989 cout<<"**************************************"<<endl;
990 cout<<"***** NON IDENT MODE ****************"<<endl;
991 cout<<"**************************************"<<endl;
992 for (Int_t i = 0;i<Nev;i++)
994 rawtrackEvent = fReader->GetTrackEvent(i);
995 if (rawtrackEvent == 0x0) continue;//in case of any error
997 /********************************/
999 /********************************/
1000 if ((ninbuffer > fBufferSize) && (fBufferSize > 0))
1001 {//if we have in buffer desired number of events, use the last. If fBufferSize<0 mix all events for background
1002 trackEvent2 = (AliHBTEvent*)tbuffer.Remove(tbuffer.Last());
1007 trackEvent2 = new AliHBTEvent();
1008 trackEvent2->SetOwner(kFALSE);
1010 FilterOut(trackEvent1, trackEvent2, rawtrackEvent);
1012 for (Int_t j = 0; j<trackEvent1->GetNumberOfParticles() ; j++)
1014 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
1016 track1= trackEvent1->GetParticle(j);
1019 for(zz = 0; zz<fNTrackMonitorFunctions; zz++)
1020 fTrackMonitorFunctions[zz]->ProcessSameEventParticles(track1);
1022 /***************************************/
1023 /****** filling numerators ********/
1024 /****** (particles from event2) ********/
1025 /***************************************/
1026 for (Int_t k = 0; k < trackEvent2->GetNumberOfParticles() ; k++)
1028 track2= trackEvent2->GetParticle(k);
1029 trackpair->SetParticles(track1,track2);
1032 if( fPairCut->PassPairProp(trackpair)) //check pair cut
1033 { //do not meets crietria of the pair cut
1037 {//meets criteria of the pair cut
1039 for(ii = 0;ii<fNTrackFunctions;ii++)
1040 fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
1043 /***************************************/
1044 /***** Filling denominators *********/
1045 /***************************************/
1046 TIter titer(&tbuffer);
1049 while ( (trackEvent3 = (AliHBTEvent*)titer.Next()) != 0x0)
1052 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i - (++nmonitor)<<endl;
1054 for (Int_t k = 0; k < trackEvent3->GetNumberOfParticles() ; k++)
1057 track2= trackEvent3->GetParticle(k);
1058 trackpair->SetParticles(track1,track2);
1061 if( fPairCut->PassPairProp(trackpair)) //check pair cut
1062 { //do not meets crietria of the pair cut
1066 {//meets criteria of the pair cut
1068 for(ii = 0;ii<fNTrackFunctions;ii++)
1069 fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
1072 }// for particles event2
1074 }//for over particles in event1
1076 tbuffer.AddFirst(trackEvent2);
1079 }//end of loop over events (1)
1083 /*************************************************************************************/
1085 void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
1087 AliHBTParticle * part1 = 0x0, * part2 = 0x0;
1089 AliHBTEvent * partEvent1 = 0x0;
1090 AliHBTEvent * partEvent2 = 0x0;
1091 AliHBTEvent * partEvent3 = 0x0;
1093 AliHBTEvent * rawpartEvent = 0x0;
1095 Int_t Nev = fReader->GetNumberOfPartEvents();
1097 AliHBTPair * partpair = new AliHBTPair();
1100 Int_t ninbuffer = 0;
1102 partEvent1 = new AliHBTEvent();
1103 partEvent1->SetOwner(kFALSE);;
1105 cout<<"**************************************"<<endl;
1106 cout<<"***** PART NON IDENT MODE ******"<<endl;
1107 cout<<"**************************************"<<endl;
1109 // gBenchmark->Start("non_id");
1110 for (Int_t i = 0;i<Nev;i++)
1112 rawpartEvent = fReader->GetParticleEvent(i);
1113 if ( rawpartEvent == 0x0 ) continue;//in case of any error
1115 /********************************/
1117 /********************************/
1118 if ((ninbuffer > fBufferSize) && (fBufferSize > 0))
1119 {//if we have in buffer desired number of events, use the last. If fBufferSize<0 mix all events for background
1120 partEvent2 = (AliHBTEvent*)pbuffer.Remove(pbuffer.Last()); //remove from the end to be reset, filled and put on the beginning
1125 partEvent2 = new AliHBTEvent();
1126 partEvent2->SetOwner(kFALSE);
1128 FilterOut(partEvent1, partEvent2, rawpartEvent);
1130 for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
1132 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
1134 part1= partEvent1->GetParticle(j);
1137 for(zz = 0; zz<fNParticleMonitorFunctions; zz++)
1138 fParticleMonitorFunctions[zz]->ProcessSameEventParticles(part1);
1140 /***************************************/
1141 /****** filling numerators ********/
1142 /****** (particles from event2) ********/
1143 /***************************************/
1144 for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++)
1146 part2= partEvent2->GetParticle(k);
1147 partpair->SetParticles(part1,part2);
1149 if(fPairCut->PassPairProp(partpair) ) //check pair cut
1150 { //do not meets crietria of the pair cut
1154 {//meets criteria of the pair cut
1156 for(ii = 0;ii<fNParticleFunctions;ii++)
1158 // if ((k%100) == 0) cout<<".";
1159 fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
1163 /***************************************/
1164 /***** Filling denominators *********/
1165 /***************************************/
1166 TIter piter(&pbuffer);
1169 while ( (partEvent3 = (AliHBTEvent*)piter.Next()) != 0x0)
1171 if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i - (++nmonitor)<<endl;
1173 for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
1175 part2= partEvent3->GetParticle(k);
1176 partpair->SetParticles(part1,part2);
1179 if(fPairCut->PassPairProp(partpair) ) //check pair cut
1180 { //do not meets crietria of the pair cut
1184 {//meets criteria of the pair cut
1186 for(ii = 0;ii<fNParticleFunctions;ii++)
1188 // if ((k%100) == 0) cout<<"*";
1189 fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
1192 }// for particles event2
1194 }//for over particles in event1
1196 pbuffer.AddFirst(partEvent2);
1199 }//end of loop over events (1)
1201 // gBenchmark->Stop("non_id");
1202 // gBenchmark->Show("non_id");
1203 pbuffer.SetOwner();//to delete stered events.
1206 /*************************************************************************************/
1207 void AliHBTAnalysis::FilterOut(AliHBTEvent* outpart1, AliHBTEvent* outpart2, AliHBTEvent* inpart,
1208 AliHBTEvent* outtrack1, AliHBTEvent* outtrack2, AliHBTEvent* intrack)
1210 //Puts particles accepted as a first particle by global cut in out1
1211 //and as a second particle in out2
1213 AliHBTParticle* part, *track;
1220 AliHBTParticleCut *cut1 = fPairCut->GetFirstPartCut();
1221 AliHBTParticleCut *cut2 = fPairCut->GetSecondPartCut();
1225 for (Int_t i = 0; i < inpart->GetNumberOfParticles(); i++)
1228 part = inpart->GetParticle(i);
1229 track = intrack->GetParticle(i);
1231 if ( (cut1->Pass(part)) || ( cut1->Pass(track) ) ) in1 = kFALSE; //if any, part or track, is rejected by cut1, in1 is false
1232 if ( (cut2->Pass(part)) || ( cut2->Pass(track) ) ) in2 = kFALSE; //if any, part or track, is rejected by cut2, in2 is false
1234 if (gDebug)//to be removed in real analysis
1235 if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1237 //Particle accpted by both cuts
1238 Error("FilterOut","Particle accepted by both cuts");
1244 outpart1->AddParticle(part);
1245 outtrack1->AddParticle(track);
1251 outpart2->AddParticle(part);
1252 outtrack2->AddParticle(track);
1260 /*************************************************************************************/
1261 void AliHBTAnalysis::FilterOut(AliHBTEvent* out1, AliHBTEvent* out2, AliHBTEvent* in)
1263 //Puts particles accepted as a first particle by global cut in out1
1264 //and as a second particle in out2
1265 AliHBTParticle* part;
1270 AliHBTParticleCut *cut1 = fPairCut->GetFirstPartCut();
1271 AliHBTParticleCut *cut2 = fPairCut->GetSecondPartCut();
1275 for (Int_t i = 0; i < in->GetNumberOfParticles(); i++)
1278 part = in->GetParticle(i);
1280 if ( cut1->Pass(part) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
1281 if ( cut2->Pass(part) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
1283 if (gDebug)//to be removed in real analysis
1284 if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1286 //Particle accpted by both cuts
1287 Error("FilterOut","Particle accepted by both cuts");
1293 out1->AddParticle(part);
1299 out2->AddParticle(part);
1304 /*************************************************************************************/
1306 Bool_t AliHBTAnalysis::IsNonIdentAnalysis()
1308 //checks if it is possible to use special analysis for non identical particles
1309 //it means - in global pair cut first particle id is different than second one
1310 //and both are different from 0
1311 //in the future is possible to perform more sophisticated check
1312 //if cuts have excluding requirements
1314 if (fPairCut->IsEmpty())
1317 if (fPairCut->GetFirstPartCut()->IsEmpty())
1320 if (fPairCut->GetSecondPartCut()->IsEmpty())
1323 Int_t id1 = fPairCut->GetFirstPartCut()->GetPID();
1324 Int_t id2 = fPairCut->GetSecondPartCut()->GetPID();
1326 if ( (id1==0) || (id2==0) || (id1==id2) )