1 #include "AliHBTAnalysis.h"
2 //_________________________________________________________
3 ////////////////////////////////////////////////////////////////////////////
5 // class AliHBTAnalysis
7 // Central Object Of HBTAnalyser:
8 // This class performs main looping within HBT Analysis
9 // User must plug a reader of Type AliReader
10 // User plugs in coorelation and monitor functions
11 // as well as monitor functions
13 // HBT Analysis Tool, which is integral part of AliRoot,
14 // ALICE Off-Line framework:
16 // Piotr.Skowronski@cern.ch
17 // more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
19 ////////////////////////////////////////////////////////////////////////////
20 //_________________________________________________________
27 #include "AliAODParticle.h"
28 #include "AliAODPairCut.h"
29 #include "AliEventCut.h"
31 #include "AliEventBuffer.h"
33 #include "AliReader.h"
34 #include "AliHBTPair.h"
35 #include "AliHBTFunction.h"
36 #include "AliHBTMonitorFunction.h"
39 ClassImp(AliHBTAnalysis)
41 const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100;
42 const UInt_t AliHBTAnalysis::fgkDefaultMixingInfo = 1000;
43 const Int_t AliHBTAnalysis::fgkDefaultBufferSize = 5;
45 AliHBTAnalysis::AliHBTAnalysis():
49 fNParticleFunctions(0),
50 fNParticleAndTrackFunctions(0),
51 fNTrackMonitorFunctions(0),
52 fNParticleMonitorFunctions(0),
53 fNParticleAndTrackMonitorFunctions(0),
54 fTrackFunctions ( new AliHBTOnePairFctn* [fgkFctnArraySize]),
55 fParticleFunctions ( new AliHBTOnePairFctn* [fgkFctnArraySize]),
56 fParticleAndTrackFunctions ( new AliHBTTwoPairFctn* [fgkFctnArraySize]),
57 fParticleMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),
58 fTrackMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),
59 fParticleAndTrackMonitorFunctions ( new AliHBTMonTwoParticleFctn* [fgkFctnArraySize]),
64 fDisplayMixingInfo(fgkDefaultMixingInfo),
66 fProcessOption(kSimulatedAndReconstructed),
76 /*************************************************************************************/
78 AliHBTAnalysis::AliHBTAnalysis(const AliHBTAnalysis& in):
83 fNParticleFunctions(0),
84 fNParticleAndTrackFunctions(0),
85 fNTrackMonitorFunctions(0),
86 fNParticleMonitorFunctions(0),
87 fNParticleAndTrackMonitorFunctions(0),
89 fParticleFunctions(0x0),
90 fParticleAndTrackFunctions(0x0),
91 fParticleMonitorFunctions(0x0),
92 fTrackMonitorFunctions(0x0),
93 fParticleAndTrackMonitorFunctions(0x0),
97 fBufferSize(fgkDefaultBufferSize),
98 fDisplayMixingInfo(fgkDefaultMixingInfo),
100 fProcessOption(kSimulatedAndReconstructed),
101 fNoCorrfctns(kFALSE),
102 fOutputFileName(0x0),
108 Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
110 /*************************************************************************************/
111 AliHBTAnalysis& AliHBTAnalysis::operator=(const AliHBTAnalysis& /*right*/)
114 Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
117 /*************************************************************************************/
118 AliHBTAnalysis::~AliHBTAnalysis()
121 //note that we do not delete functions itself
122 // they should be deleted by whom where created
123 //we only store pointers, and we use "new" only for pointers array
127 if (AliVAODParticle::GetDebug()>5)Info("~AliHBTAnalysis","Is Owner: Attempting to delete functions");
129 if (AliVAODParticle::GetDebug()>5)Info("~AliHBTAnalysis","Delete functions done");
131 delete [] fTrackFunctions;
132 delete [] fParticleFunctions;
133 delete [] fParticleAndTrackFunctions;
135 delete [] fParticleMonitorFunctions;
136 delete [] fTrackMonitorFunctions;
137 delete [] fParticleAndTrackMonitorFunctions;
140 delete fOutputFileName;
143 /*************************************************************************************/
145 Int_t AliHBTAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim)
147 //Processes one event
148 if (fProcEvent == 0x0)
150 Error("ProcessEvent","Analysis <<%s>> option not specified.",GetName());
153 if ( Rejected(aodrec,aodsim) )
155 // Double_t x = 0, y = 0, z = 0;
156 // aodrec->GetPrimaryVertex(x,y,z);
157 // Info("ProcessEvent","Event has vertex at %f %f %f",x,y,z);
158 Info("ProcessEvent","Nch is %d",aodsim->GetNumberOfCharged());
159 Info("ProcessEvent","%s: Event cut rejected this event",GetName());
163 //Move event to the apparent vertex -> must be after the event cut
164 //It is important for any cut that use any spacial coordiantes,
165 //f.g. anti-merging cut in order to preserve the same bahavior of variables (f.g. distance between tracks)
166 Double_t dvx = 0, dvy = 0, dvz = 0;
169 Double_t pvx,pvy,pvz;
170 aodrec->GetPrimaryVertex(pvx,pvy,pvz);
172 dvx = fVertexX - pvx;
173 dvy = fVertexY - pvy;
174 dvz = fVertexZ - pvz;
175 aodrec->Move(dvx,dvy,dvz);
178 Int_t result = (this->*fProcEvent)(aodrec,aodsim);
180 if (aodrec) aodrec->Move(-dvx,-dvy,-dvz);//move event back to the same spacial coordinates
184 /*************************************************************************************/
186 Int_t AliHBTAnalysis::Finish()
192 /*************************************************************************************/
194 void AliHBTAnalysis::DeleteFunctions()
196 //Deletes all functions added to analysis
199 for(ii = 0;ii<fNParticleFunctions;ii++)
201 if (AliVAODParticle::GetDebug()>5)
203 Info("DeleteFunctions","Deleting ParticleFunction %#x",fParticleFunctions[ii]);
204 Info("DeleteFunctions","Deleting ParticleFunction %s",fParticleFunctions[ii]->Name());
206 delete fParticleFunctions[ii];
208 fNParticleFunctions = 0;
210 for(ii = 0;ii<fNTrackFunctions;ii++)
212 if (AliVAODParticle::GetDebug()>5)
214 Info("DeleteFunctions","Deleting TrackFunction %#x",fTrackFunctions[ii]);
215 Info("DeleteFunctions","Deleting TrackFunction %s",fTrackFunctions[ii]->Name());
217 delete fTrackFunctions[ii];
219 fNTrackFunctions = 0;
221 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
223 if (AliVAODParticle::GetDebug()>5)
225 Info("DeleteFunctions","Deleting ParticleAndTrackFunction %#x",fParticleAndTrackFunctions[ii]);
226 Info("DeleteFunctions","Deleting ParticleAndTrackFunction %s",fParticleAndTrackFunctions[ii]->Name());
228 delete fParticleAndTrackFunctions[ii];
230 fNParticleAndTrackFunctions = 0;
232 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
234 if (AliVAODParticle::GetDebug()>5)
236 Info("DeleteFunctions","Deleting ParticleMonitorFunction %#x",fParticleMonitorFunctions[ii]);
237 Info("DeleteFunctions","Deleting ParticleMonitorFunction %s",fParticleMonitorFunctions[ii]->Name());
239 delete fParticleMonitorFunctions[ii];
241 fNParticleMonitorFunctions = 0;
243 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
245 if (AliVAODParticle::GetDebug()>5)
247 Info("DeleteFunctions","Deleting TrackMonitorFunction %#x",fTrackMonitorFunctions[ii]);
248 Info("DeleteFunctions","Deleting TrackMonitorFunction %s",fTrackMonitorFunctions[ii]->Name());
250 delete fTrackMonitorFunctions[ii];
252 fNTrackMonitorFunctions = 0;
254 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
256 if (AliVAODParticle::GetDebug()>5)
258 Info("DeleteFunctions","Deleting ParticleAndTrackMonitorFunction %#x",fParticleAndTrackMonitorFunctions[ii]);
259 Info("DeleteFunctions","Deleting ParticleAndTrackMonitorFunction %s",fParticleAndTrackMonitorFunctions[ii]->Name());
261 delete fParticleAndTrackMonitorFunctions[ii];
263 fNParticleAndTrackMonitorFunctions = 0;
266 /*************************************************************************************/
268 Int_t AliHBTAnalysis::Init()
270 //Initializeation method
271 //calls Init for all functions
273 //Depending on option and pair cut assigns proper analysis method
274 Bool_t nonid = IsNonIdentAnalysis();
275 switch(fProcessOption)
278 if (nonid) fProcEvent = &AliHBTAnalysis::ProcessRecNonId;
279 else fProcEvent = &AliHBTAnalysis::ProcessRec;
284 if (nonid) fProcEvent = &AliHBTAnalysis::ProcessSimNonId;
285 else fProcEvent = &AliHBTAnalysis::ProcessSim;
289 case kSimulatedAndReconstructed:
290 if (nonid) fProcEvent = &AliHBTAnalysis::ProcessRecAndSimNonId;
291 else fProcEvent = &AliHBTAnalysis::ProcessRecAndSim;
295 if (fPartBuffer == 0x0) fPartBuffer = new AliEventBuffer (fBufferSize);
296 else fPartBuffer->Reset();
298 if (fTrackBuffer == 0x0) fTrackBuffer = new AliEventBuffer (fBufferSize);
299 else fTrackBuffer->Reset();
302 fNoCorrfctns = (fNParticleFunctions == 0) && (fNTrackFunctions == 0) && (fNParticleAndTrackFunctions == 0);
305 for(ii = 0;ii<fNParticleFunctions;ii++)
306 fParticleFunctions[ii]->Init();
308 for(ii = 0;ii<fNTrackFunctions;ii++)
309 fTrackFunctions[ii]->Init();
311 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
312 fParticleAndTrackFunctions[ii]->Init();
314 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
315 fParticleMonitorFunctions[ii]->Init();
317 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
318 fTrackMonitorFunctions[ii]->Init();
320 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
321 fParticleAndTrackMonitorFunctions[ii]->Init();
325 /*************************************************************************************/
327 void AliHBTAnalysis::ResetFunctions()
329 //In case fOwner is true, deletes all functions
330 //in other case, just set number of analysis to 0
331 if (fIsOwner) DeleteFunctions();
334 fNParticleFunctions = 0;
335 fNTrackFunctions = 0;
336 fNParticleAndTrackFunctions = 0;
337 fNParticleMonitorFunctions = 0;
338 fNTrackMonitorFunctions = 0;
339 fNParticleAndTrackMonitorFunctions = 0;
342 /*************************************************************************************/
343 Int_t AliHBTAnalysis::ProcessRecAndSim(AliAOD* aodrec, AliAOD* aodsim)
345 //Does analysis for both tracks and particles
346 //mainly for resolution study and analysies with weighting algirithms
348 // cut on particles only -- why?
349 // - PID: when we make resolution analysis we want to take only tracks with correct PID
350 // We need cut on tracks because there are data characteristic
352 AliVAODParticle * part1, * part2;
353 AliVAODParticle * track1, * track2;
355 AliAOD * trackEvent = aodrec, *partEvent = aodsim;
356 AliAOD* trackEvent1 = new AliAOD();
357 AliAOD* partEvent1 = new AliAOD();
359 AliAOD * trackEvent2,*partEvent2;
361 // Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
363 // Int_t nev = fReader->GetNumberOfTrackEvents();
364 static AliHBTPair tpair;
365 static AliHBTPair ppair;
367 AliHBTPair* trackpair = &tpair;
368 AliHBTPair* partpair = &ppair;
370 AliHBTPair * tmptrackpair;//temprary pointers to pairs
371 AliHBTPair * tmppartpair;
377 if ( !partEvent || !trackEvent )
379 Error("ProcessRecAndSim","<<%s>> Can not get event",GetName());
383 if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
385 Error("ProcessRecAndSim",
386 "Number of simulated particles (%d) not equal to number of reconstructed tracks (%d). Skipping Event.",
387 partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
392 for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
394 /***************************************/
395 /****** Looping same events ********/
396 /****** filling numerators ********/
397 /***************************************/
398 if ( (j%fDisplayMixingInfo) == 0)
399 Info("ProcessTracksAndParticles",
400 "Mixing particle %d with particles from the same event",j);
402 part1= partEvent->GetParticle(j);
403 track1= trackEvent->GetParticle(j);
405 Bool_t firstcut = (this->*fkPass1)(part1,track1);
406 if (fBufferSize != 0)
407 if ( (firstcut == kFALSE) || ( (this->*fkPass2)(part1,track1) == kFALSE ) )
409 //accepted by any cut
410 // we have to copy because reader keeps only one event
412 partEvent1->AddParticle(part1);
413 trackEvent1->AddParticle(track1);
416 if (firstcut) continue;
418 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
419 fParticleMonitorFunctions[ii]->Process(part1);
420 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
421 fTrackMonitorFunctions[ii]->Process(track1);
422 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
423 fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
425 if (fNoCorrfctns) continue;
427 for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
429 part2= partEvent->GetParticle(k);
430 if (part1->GetUID() == part2->GetUID()) continue;
431 partpair->SetParticles(part1,part2);
433 track2= trackEvent->GetParticle(k);
434 trackpair->SetParticles(track1,track2);
436 if( (this->*fkPass)(partpair,trackpair) ) //check pair cut
437 { //do not meets crietria of the pair cut, try with swapped pairs
438 if( (this->*fkPass)((AliHBTPair*)partpair->GetSwappedPair(),(AliHBTPair*)trackpair->GetSwappedPair()) )
439 continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
441 { //swaped pair meets all the criteria
442 tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
443 tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
447 {//meets criteria of the pair cut
448 tmptrackpair = trackpair;
449 tmppartpair = partpair;
452 for(ii = 0;ii<fNParticleFunctions;ii++)
453 fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
455 for(ii = 0;ii<fNTrackFunctions;ii++)
456 fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
458 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
459 fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair,tmppartpair);
460 //end of 2nd loop over particles from the same event
461 }//for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
463 /***************************************/
464 /***** Filling denominators *********/
465 /***************************************/
466 if (fBufferSize == 0) continue;
468 fPartBuffer->ResetIter();
469 fTrackBuffer->ResetIter();
471 while (( partEvent2 = fPartBuffer->Next() ))
473 trackEvent2 = fTrackBuffer->Next();
476 if ( (j%fDisplayMixingInfo) == 0)
477 Info("ProcessTracksAndParticles",
478 "Mixing particle %d from current event with particles from event %d",j,-m);
480 for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
482 part2= partEvent2->GetParticle(l);
483 partpair->SetParticles(part1,part2);
485 track2= trackEvent2->GetParticle(l);
486 trackpair->SetParticles(track1,track2);
488 if( (this->*fkPass)(partpair,trackpair) ) //check pair cut
489 { //do not meets crietria of the
490 if( (this->*fkPass)((AliHBTPair*)partpair->GetSwappedPair(),(AliHBTPair*)trackpair->GetSwappedPair()) )
494 tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
495 tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
499 {//meets criteria of the pair cut
500 tmptrackpair = trackpair;
501 tmppartpair = partpair;
504 for(ii = 0;ii<fNParticleFunctions;ii++)
505 fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
507 for(ii = 0;ii<fNTrackFunctions;ii++)
508 fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
510 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
511 fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair);
512 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
515 //end of loop over particles from first event
516 }//for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
517 delete fPartBuffer->Push(partEvent1);
518 delete fTrackBuffer->Push(trackEvent1);
519 //end of loop over events
522 /*************************************************************************************/
524 Int_t AliHBTAnalysis::ProcessSim(AliAOD* /*aodrec*/, AliAOD* aodsim)
526 //Does analysis of simulated data
527 AliVAODParticle * part1, * part2;
529 AliAOD* partEvent = aodsim;
530 AliAOD* partEvent1 = new AliAOD();
536 AliHBTPair* partpair = &ppair;
538 AliHBTPair * tmppartpair;
546 Error("ProcessRecAndSim","Can not get event");
551 for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
553 /***************************************/
554 /****** Looping same events ********/
555 /****** filling numerators ********/
556 /***************************************/
557 if ( (j%fDisplayMixingInfo) == 0)
558 Info("ProcessTracksAndParticles",
559 "Mixing particle %d with particles from the same event",j);
561 part1= partEvent->GetParticle(j);
563 Bool_t firstcut = fPairCut->GetFirstPartCut()->Rejected(part1);
565 if (fBufferSize != 0)
566 if ( (firstcut == kFALSE) || ( fPairCut->GetSecondPartCut()->Rejected(part1) == kFALSE ) )
568 //accepted by any cut
569 // we have to copy because reader keeps only one event
571 partEvent1->AddParticle(part1);
574 if (firstcut) continue;
576 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
577 fParticleMonitorFunctions[ii]->Process(part1);
579 if ( fNParticleFunctions == 0 ) continue;
581 for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
583 part2= partEvent->GetParticle(k);
584 if (part1->GetUID() == part2->GetUID()) continue;
585 partpair->SetParticles(part1,part2);
587 if(fPairCut->Rejected(partpair)) //check pair cut
588 { //do not meets crietria of the
589 if( fPairCut->Rejected((AliHBTPair*)partpair->GetSwappedPair()) ) continue;
590 else tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
594 tmppartpair = partpair;
597 for(ii = 0;ii<fNParticleFunctions;ii++)
598 fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
600 //end of 2nd loop over particles from the same event
601 }//for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
603 /***************************************/
604 /***** Filling denominators *********/
605 /***************************************/
606 if (fBufferSize == 0) continue;
608 fPartBuffer->ResetIter();
610 while (( partEvent2 = fPartBuffer->Next() ))
613 if ( (j%fDisplayMixingInfo) == 0)
614 Info("ProcessParticles",
615 "Mixing particle %d from current event with particles from event %d",j,-m);
616 for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
619 part2= partEvent2->GetParticle(l);
620 partpair->SetParticles(part1,part2);
622 if( fPairCut->Rejected(partpair) ) //check pair cut
623 { //do not meets crietria of the
624 if( fPairCut->Rejected((AliHBTPair*)partpair->GetSwappedPair()) )
628 tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
632 {//meets criteria of the pair cut
633 tmppartpair = partpair;
636 for(ii = 0;ii<fNParticleFunctions;ii++)
637 fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
639 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
642 delete fPartBuffer->Push(partEvent1);
643 //end of loop over events
646 /*************************************************************************************/
647 Int_t AliHBTAnalysis::ProcessRec(AliAOD* aodrec, AliAOD* /*aodsim*/)
649 //Does analysis of reconstructed data
650 AliVAODParticle * track1, * track2;
652 AliAOD* trackEvent = aodrec;
653 AliAOD* trackEvent1 = new AliAOD();
659 AliHBTPair* trackpair = &tpair;
661 AliHBTPair * tmptrackpair;
668 Error("ProcessRecAndSim","Can not get event");
673 for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
675 /***************************************/
676 /****** Looping same events ********/
677 /****** filling numerators ********/
678 /***************************************/
679 if ( (j%fDisplayMixingInfo) == 0)
680 Info("ProcessTracksAndParticles",
681 "Mixing Particle %d with Particles from the same event",j);
683 track1= trackEvent->GetParticle(j);
685 Bool_t firstcut = fPairCut->GetFirstPartCut()->Rejected(track1);
687 if (fBufferSize != 0)
688 if ( (firstcut == kFALSE) || ( fPairCut->GetSecondPartCut()->Rejected(track1) == kFALSE ) )
690 //accepted by any cut
691 // we have to copy because reader keeps only one event
693 trackEvent1->AddParticle(track1);
696 if (firstcut) continue;
698 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
699 fParticleMonitorFunctions[ii]->Process(track1);
701 if ( fNParticleFunctions == 0 ) continue;
703 for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
705 track2= trackEvent->GetParticle(k);
706 if (track1->GetUID() == track2->GetUID()) continue;
707 trackpair->SetParticles(track1,track2);
709 if(fPairCut->Rejected(trackpair)) //check pair cut
710 { //do not meets crietria of the
711 if( fPairCut->Rejected((AliHBTPair*)trackpair->GetSwappedPair()) ) continue;
712 else tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
716 tmptrackpair = trackpair;
719 for(ii = 0;ii<fNTrackFunctions;ii++)
720 fParticleFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
722 //end of 2nd loop over Particles from the same event
723 }//for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
725 /***************************************/
726 /***** Filling denominators *********/
727 /***************************************/
728 if (fBufferSize == 0) continue;
730 fTrackBuffer->ResetIter();
732 while (( trackEvent2 = fTrackBuffer->Next() ))
735 if ( (j%fDisplayMixingInfo) == 0)
736 Info("ProcessParticles",
737 "Mixing Particle %d from current event with Particles from event %d",j,-m);
738 for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++) // ... on all Particles
741 track2= trackEvent2->GetParticle(l);
742 trackpair->SetParticles(track1,track2);
744 if( fPairCut->Rejected(trackpair) ) //check pair cut
745 { //do not meets crietria of the
746 if( fPairCut->Rejected((AliHBTPair*)trackpair->GetSwappedPair()) )
750 tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
754 {//meets criteria of the pair cut
755 tmptrackpair = trackpair;
758 for(ii = 0;ii<fNTrackFunctions;ii++)
759 fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
761 }//for(Int_t l = 0; l<N2;l++) // ... on all Particles
764 delete fTrackBuffer->Push(trackEvent1);
765 //end of loop over events
768 /*************************************************************************************/
770 Int_t AliHBTAnalysis::ProcessRecAndSimNonId(AliAOD* aodrec, AliAOD* aodsim)
772 //Analyzes both reconstructed and simulated data
775 Error("ProcessTracksAndParticlesNonIdentAnal","Reconstructed event is NULL");
781 Error("ProcessTracksAndParticlesNonIdentAnal","Simulated event is NULL");
785 if ( aodrec->GetNumberOfParticles() != aodsim->GetNumberOfParticles() )
787 Error("ProcessTracksAndParticlesNonIdentAnal",
788 "Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
789 aodsim->GetNumberOfParticles() , aodrec->GetNumberOfParticles());
794 AliVAODParticle * part1, * part2;
795 AliVAODParticle * track1, * track2;
797 static AliAOD aodrec1;
798 static AliAOD aodsim1;
800 AliAOD * trackEvent1=&aodrec1,*partEvent1=&aodsim1;//Particle that passes first particle cut, this event
801 trackEvent1->Reset();
803 AliAOD * trackEvent2=0x0,*partEvent2=0x0;//Particle that passes second particle cut, this event
804 AliAOD * trackEvent3=0x0,*partEvent3=0x0;//Particle that passes second particle cut, events from buffer
806 AliAOD* rawtrackEvent = aodrec;//this we get from Reader
807 AliAOD* rawpartEvent = aodsim;//this we get from Reader
809 static AliHBTPair tpair;
810 static AliHBTPair ppair;
812 AliHBTPair* trackpair = &tpair;
813 AliHBTPair* partpair = &ppair;
817 /********************************/
819 /********************************/
820 if ( ( (partEvent2==0x0) || (trackEvent2==0x0)) )//in case fBufferSize == 0 and pointers are created do not eneter
822 partEvent2 = new AliAOD();
823 trackEvent2 = new AliAOD();
826 FilterOut(partEvent1, partEvent2, rawpartEvent, trackEvent1, trackEvent2, rawtrackEvent);
828 for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
830 if ( (j%fDisplayMixingInfo) == 0)
831 Info("ProcessTracksAndParticlesNonIdentAnal",
832 "Mixing particle %d from current event with particles from current event",j);
834 part1= partEvent1->GetParticle(j);
835 track1= trackEvent1->GetParticle(j);
838 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
839 fParticleMonitorFunctions[ii]->Process(part1);
840 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
841 fTrackMonitorFunctions[ii]->Process(track1);
842 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
843 fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
845 if (fNoCorrfctns) continue;
847 /***************************************/
848 /****** filling numerators ********/
849 /****** (particles from event2) ********/
850 /***************************************/
852 for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups
854 part2= partEvent2->GetParticle(k);
855 if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
856 partpair->SetParticles(part1,part2);
858 track2= trackEvent2->GetParticle(k);
859 trackpair->SetParticles(track1,track2);
861 if( (this->*fkPassPairProp)(partpair,trackpair) ) //check pair cut
862 { //do not meets crietria of the pair cut
866 {//meets criteria of the pair cut
867 for(ii = 0;ii<fNParticleFunctions;ii++)
868 fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
870 for(ii = 0;ii<fNTrackFunctions;ii++)
871 fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
873 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
874 fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(trackpair,partpair);
878 if ( fBufferSize == 0) continue;//do not mix diff histograms
879 /***************************************/
880 /***** Filling denominators *********/
881 /***************************************/
882 fPartBuffer->ResetIter();
883 fTrackBuffer->ResetIter();
887 while ( (partEvent3 = fPartBuffer->Next() ) != 0x0)
889 trackEvent3 = fTrackBuffer->Next();
891 if ( (j%fDisplayMixingInfo) == 0)
892 Info("ProcessTracksAndParticlesNonIdentAnal",
893 "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
895 for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
897 part2= partEvent3->GetParticle(k);
898 partpair->SetParticles(part1,part2);
900 track2= trackEvent3->GetParticle(k);
901 trackpair->SetParticles(track1,track2);
903 if( (this->*fkPassPairProp)(partpair,trackpair) ) //check pair cut
904 { //do not meets crietria of the pair cut
908 {//meets criteria of the pair cut
910 for(ii = 0;ii<fNParticleFunctions;ii++)
911 fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
913 for(ii = 0;ii<fNTrackFunctions;ii++)
914 fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
916 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
917 fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(trackpair,partpair);
919 }// for particles event2
921 }//for over particles in event1
923 delete fPartBuffer->Push(partEvent2);
924 delete fTrackBuffer->Push(trackEvent2);
928 /*************************************************************************************/
929 Int_t AliHBTAnalysis::ProcessSimNonId(AliAOD* /*aodrec*/, AliAOD* aodsim)
931 //does analysis of simulated (MC) data in non-identical mode
932 //i.e. when particles selected by first part. cut are a disjunctive set than particles
933 //passed by the second part. cut
940 AliVAODParticle * part1, * part2;
942 static AliAOD aodsim1;
944 AliAOD* partEvent1=&aodsim1;//Particle that passes first particle cut, this event
946 AliAOD* partEvent2=0x0;//Particle that passes second particle cut, this event
947 AliAOD* partEvent3=0x0;//Particle that passes second particle cut, events from buffer
949 AliAOD* rawpartEvent = aodsim;//this we get from Reader
951 static AliHBTPair ppair;
953 AliHBTPair* partpair = &ppair;
957 /********************************/
959 /********************************/
960 if (partEvent2==0x0)//in case fBufferSize == 0 and pointers are created do not eneter
962 partEvent2 = new AliAOD();
965 FilterOut(partEvent1, partEvent2, rawpartEvent);
967 for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
969 if ( (j%fDisplayMixingInfo) == 0)
970 Info("ProcessParticlesNonIdentAnal",
971 "Mixing particle %d from current event with particles from current event",j);
973 part1= partEvent1->GetParticle(j);
976 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
977 fParticleMonitorFunctions[ii]->Process(part1);
979 if (fNParticleFunctions == 0) continue;
981 /***************************************/
982 /****** filling numerators ********/
983 /****** (particles from event2) ********/
984 /***************************************/
986 for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups
988 part2= partEvent2->GetParticle(k);
989 if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
990 partpair->SetParticles(part1,part2);
993 if(fPairCut->PassPairProp(partpair) ) //check pair cut
994 { //do not meets crietria of the pair cut
998 {//meets criteria of the pair cut
999 for(ii = 0;ii<fNParticleFunctions;ii++)
1000 fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
1004 if ( fBufferSize == 0) continue;//do not mix diff histograms
1005 /***************************************/
1006 /***** Filling denominators *********/
1007 /***************************************/
1008 fPartBuffer->ResetIter();
1012 while ( (partEvent3 = fPartBuffer->Next() ) != 0x0)
1015 if ( (j%fDisplayMixingInfo) == 0)
1016 Info("ProcessParticlesNonIdentAnal",
1017 "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
1019 for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
1021 part2= partEvent3->GetParticle(k);
1022 partpair->SetParticles(part1,part2);
1025 if(fPairCut->PassPairProp(partpair) ) //check pair cut
1026 { //do not meets crietria of the pair cut
1030 {//meets criteria of the pair cut
1031 for(ii = 0;ii<fNParticleFunctions;ii++)
1033 fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
1036 }// for particles event2
1038 }//for over particles in event1
1040 delete fPartBuffer->Push(partEvent2);
1044 /*************************************************************************************/
1045 Int_t AliHBTAnalysis::ProcessRecNonId(AliAOD* aodrec, AliAOD* /*aodsim*/)
1047 //Analyzes both reconstructed and simulated data
1053 AliVAODParticle * track1, * track2;
1055 static AliAOD aodrec1;
1056 AliAOD * trackEvent1=&aodrec1;//Particle that passes first particle cut, this event
1057 trackEvent1->Reset();
1058 AliAOD * trackEvent2=0x0;//Particle that passes second particle cut, this event
1059 AliAOD * trackEvent3=0x0;//Particle that passes second particle cut, events from buffer
1060 AliAOD* rawtrackEvent = aodrec;//this we get from Reader
1062 static AliHBTPair tpair;
1064 AliHBTPair* trackpair = &tpair;
1069 /********************************/
1071 /********************************/
1072 if ( trackEvent2==0x0 )//in case fBufferSize == 0 and pointers are created do not eneter
1074 trackEvent2 = new AliAOD();
1077 FilterOut(trackEvent1, trackEvent2, rawtrackEvent);
1079 for (Int_t j = 0; j<trackEvent1->GetNumberOfParticles() ; j++)
1081 if ( (j%fDisplayMixingInfo) == 0)
1082 Info("ProcessTracksNonIdentAnal",
1083 "Mixing particle %d from current event with particles from current event",j);
1085 track1= trackEvent1->GetParticle(j);
1088 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
1089 fTrackMonitorFunctions[ii]->Process(track1);
1091 if (fNTrackFunctions == 0x0) continue;
1093 /***************************************/
1094 /****** filling numerators ********/
1095 /****** (particles from event2) ********/
1096 /***************************************/
1098 for (Int_t k = 0; k < trackEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups
1100 track2= trackEvent2->GetParticle(k);
1101 if (track1->GetUID() == track2->GetUID()) continue;//this is the same particle but with different PID
1102 trackpair->SetParticles(track1,track2);
1105 if( fPairCut->PassPairProp(trackpair)) //check pair cut
1106 { //do not meets crietria of the pair cut
1110 {//meets criteria of the pair cut
1112 for(ii = 0;ii<fNTrackFunctions;ii++)
1113 fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
1117 if ( fBufferSize == 0) continue;//do not mix diff histograms
1118 /***************************************/
1119 /***** Filling denominators *********/
1120 /***************************************/
1121 fTrackBuffer->ResetIter();
1125 while ( (trackEvent3 = fTrackBuffer->Next() ) != 0x0)
1127 if ( (j%fDisplayMixingInfo) == 0)
1128 Info("ProcessTracksNonIdentAnal",
1129 "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
1131 for (Int_t k = 0; k < trackEvent3->GetNumberOfParticles() ; k++)
1133 track2= trackEvent3->GetParticle(k);
1134 trackpair->SetParticles(track1,track2);
1136 if( fPairCut->PassPairProp(trackpair)) //check pair cut
1137 { //do not meets crietria of the pair cut
1141 {//meets criteria of the pair cut
1142 for(ii = 0;ii<fNTrackFunctions;ii++)
1143 fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
1145 }// for particles event2
1147 }//for over particles in event1
1149 delete fTrackBuffer->Push(trackEvent2);
1153 /*************************************************************************************/
1155 void AliHBTAnalysis::Process(Option_t* option)
1157 //default option = "TracksAndParticles"
1158 //Main method of the HBT Analysis Package
1159 //It triggers reading with the global cut (default is an empty cut)
1160 //Than it checks options and data which are read
1161 //if everything is OK, then it calls one of the looping methods
1162 //depending on tfReaderhe option
1163 //These methods differs on what thay are looping on
1166 //--------------------------------------------------------------------
1167 //ProcessTracksAndParticles - "TracksAndParticles"
1169 // it loops over both, tracks(reconstructed) and particles(simulated)
1170 // all function gethered in all 3 lists are called for each (double)pair
1172 //ProcessTracks - "Tracks"
1173 // it loops only on tracks(reconstructed),
1174 // functions ONLY from fTrackFunctions list are called
1176 //ProcessParticles - "Particles"
1177 // it loops only on particles(simulated),
1178 // functions ONLY from fParticleAndTrackFunctions list are called
1183 Error("Process","The reader is not set");
1187 const char *oT = strstr(option,"Tracks");
1188 const char *oP = strstr(option,"Particles");
1190 Bool_t nonid = IsNonIdentAnalysis();
1196 if (nonid) ProcessTracksAndParticlesNonIdentAnal();
1197 else ProcessTracksAndParticles();
1203 if (nonid) ProcessTracksNonIdentAnal();
1204 else ProcessTracks();
1210 if (nonid) ProcessParticlesNonIdentAnal();
1211 else ProcessParticles();
1216 /*************************************************************************************/
1218 void AliHBTAnalysis::ProcessTracksAndParticles()
1220 //Makes analysis for both tracks and particles
1221 //mainly for resolution study and analysies with weighting algirithms
1222 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1223 //the loops are splited
1225 // cut on particles only -- why?
1226 // - PID: when we make resolution analysis we want to take only tracks with correct PID
1227 // We need cut on tracks because there are data characteristic to
1229 AliAOD * trackEvent, *partEvent;
1232 while (fReader->Next() == kFALSE)
1234 partEvent = fReader->GetEventSim();
1235 trackEvent = fReader->GetEventRec();
1236 ProcessRecAndSim(trackEvent,partEvent);
1238 }//while (fReader->Next() == kFALSE)
1241 /*************************************************************************************/
1243 void AliHBTAnalysis::ProcessTracks()
1245 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1246 //the loops are splited
1247 AliAOD * trackEvent;
1249 while (fReader->Next() == kFALSE)
1251 trackEvent = fReader->GetEventRec();
1252 ProcessRec(trackEvent,0x0);
1253 }//while (fReader->Next() == kFALSE)
1256 /*************************************************************************************/
1258 void AliHBTAnalysis::ProcessParticles()
1260 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1261 //the loops are splited
1264 while (fReader->Next() == kFALSE)
1266 partEvent = fReader->GetEventSim();
1267 ProcessSim(0x0,partEvent);
1268 }//while (fReader->Next() == kFALSE)
1270 /*************************************************************************************/
1272 void AliHBTAnalysis::WriteFunctions()
1274 //Calls Write for all defined functions in analysis
1275 //== writes all results
1276 TFile* oututfile = 0x0;
1277 if (fOutputFileName)
1279 oututfile = TFile::Open(*fOutputFileName,"update");
1282 for(ii = 0;ii<fNParticleFunctions;ii++)
1284 if (AliVAODParticle::GetDebug()>5)
1286 Info("WriteFunctions","Writing ParticleFunction %#x",fParticleFunctions[ii]);
1287 Info("WriteFunctions","Writing ParticleFunction %s",fParticleFunctions[ii]->Name());
1289 fParticleFunctions[ii]->Write();
1292 for(ii = 0;ii<fNTrackFunctions;ii++)
1294 if (AliVAODParticle::GetDebug()>5)
1296 Info("WriteFunctions","Writing TrackFunction %#x",fTrackFunctions[ii]);
1297 Info("WriteFunctions","Writing TrackFunction %s",fTrackFunctions[ii]->Name());
1299 fTrackFunctions[ii]->Write();
1302 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
1304 if (AliVAODParticle::GetDebug()>5)
1306 Info("WriteFunctions","Writing ParticleAndTrackFunction %#x",fParticleAndTrackFunctions[ii]);
1307 Info("WriteFunctions","Writing ParticleAndTrackFunction %s",fParticleAndTrackFunctions[ii]->Name());
1309 fParticleAndTrackFunctions[ii]->Write();
1312 for(ii = 0;ii<fNParticleMonitorFunctions;ii++)
1314 if (AliVAODParticle::GetDebug()>5)
1316 Info("WriteFunctions","Writing ParticleMonitorFunction %#x",fParticleMonitorFunctions[ii]);
1317 Info("WriteFunctions","Writing ParticleMonitorFunction %s",fParticleMonitorFunctions[ii]->Name());
1319 fParticleMonitorFunctions[ii]->Write();
1322 for(ii = 0;ii<fNTrackMonitorFunctions;ii++)
1324 if (AliVAODParticle::GetDebug()>5)
1326 Info("WriteFunctions","Writing TrackMonitorFunction %#x",fTrackMonitorFunctions[ii]);
1327 Info("WriteFunctions","Writing TrackMonitorFunction %s",fTrackMonitorFunctions[ii]->Name());
1329 fTrackMonitorFunctions[ii]->Write();
1332 for(ii = 0;ii<fNParticleAndTrackMonitorFunctions;ii++)
1334 if (AliVAODParticle::GetDebug()>5)
1336 Info("WriteFunctions","Writing ParticleAndTrackMonitorFunction %#x",fParticleAndTrackMonitorFunctions[ii]);
1337 Info("WriteFunctions","Writing ParticleAndTrackMonitorFunction %s",fParticleAndTrackMonitorFunctions[ii]->Name());
1339 fParticleAndTrackMonitorFunctions[ii]->Write();
1343 /*************************************************************************************/
1345 void AliHBTAnalysis::SetOutputFileName(const char* fname)
1347 //Sets fiele name where to dump results,
1348 //if not specified reults are written to gDirectory
1351 delete fOutputFileName;
1352 fOutputFileName = 0x0;
1355 if ( strcmp(fname,"") == 0 )
1357 delete fOutputFileName;
1358 fOutputFileName = 0x0;
1361 if (fOutputFileName == 0x0) fOutputFileName = new TString(fname);
1362 else *fOutputFileName = fname;
1364 /*************************************************************************************/
1366 void AliHBTAnalysis::SetGlobalPairCut(AliAODPairCut* cut)
1368 //Sets the global cut
1371 Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
1374 fPairCut = (AliAODPairCut*)cut->Clone();
1377 /*************************************************************************************/
1379 void AliHBTAnalysis::AddTrackFunction(AliHBTOnePairFctn* f)
1381 //Adds track function
1382 if (f == 0x0) return;
1383 if (fNTrackFunctions == fgkFctnArraySize)
1385 Error("AliHBTAnalysis::AddTrackFunction","Can not add this function, not enough place in the array.");
1387 fTrackFunctions[fNTrackFunctions] = f;
1390 /*************************************************************************************/
1392 void AliHBTAnalysis::AddParticleFunction(AliHBTOnePairFctn* f)
1394 //adds particle function
1395 if (f == 0x0) return;
1397 if (fNParticleFunctions == fgkFctnArraySize)
1399 Error("AliHBTAnalysis::AddParticleFunction","Can not add this function, not enough place in the array.");
1401 fParticleFunctions[fNParticleFunctions] = f;
1402 fNParticleFunctions++;
1404 /*************************************************************************************/
1406 void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTTwoPairFctn* f)
1408 //add resolution function
1409 if (f == 0x0) return;
1410 if (fNParticleAndTrackFunctions == fgkFctnArraySize)
1412 Error("AliHBTAnalysis::AddParticleAndTrackFunction","Can not add this function, not enough place in the array.");
1414 fParticleAndTrackFunctions[fNParticleAndTrackFunctions] = f;
1415 fNParticleAndTrackFunctions++;
1417 /*************************************************************************************/
1419 void AliHBTAnalysis::AddParticleMonitorFunction(AliHBTMonOneParticleFctn* f)
1421 //add particle monitoring function
1422 if (f == 0x0) return;
1424 if (fNParticleMonitorFunctions == fgkFctnArraySize)
1426 Error("AliHBTAnalysis::AddParticleMonitorFunction","Can not add this function, not enough place in the array.");
1428 fParticleMonitorFunctions[fNParticleMonitorFunctions] = f;
1429 fNParticleMonitorFunctions++;
1431 /*************************************************************************************/
1433 void AliHBTAnalysis::AddTrackMonitorFunction(AliHBTMonOneParticleFctn* f)
1435 //add track monitoring function
1436 if (f == 0x0) return;
1438 if (fNTrackMonitorFunctions == fgkFctnArraySize)
1440 Error("AliHBTAnalysis::AddTrackMonitorFunction","Can not add this function, not enough place in the array.");
1442 fTrackMonitorFunctions[fNTrackMonitorFunctions] = f;
1443 fNTrackMonitorFunctions++;
1445 /*************************************************************************************/
1447 void AliHBTAnalysis::AddParticleAndTrackMonitorFunction(AliHBTMonTwoParticleFctn* f)
1449 //add resolution monitoring function
1450 if (f == 0x0) return;
1451 if (fNParticleAndTrackMonitorFunctions == fgkFctnArraySize)
1453 Error("AliHBTAnalysis::AddParticleAndTrackMonitorFunction","Can not add this function, not enough place in the array.");
1455 fParticleAndTrackMonitorFunctions[fNParticleAndTrackMonitorFunctions] = f;
1456 fNParticleAndTrackMonitorFunctions++;
1460 /*************************************************************************************/
1461 /*************************************************************************************/
1463 Bool_t AliHBTAnalysis::RunCoherencyCheck()
1465 //Checks if both HBTRuns are similar
1466 //return true if error found
1467 //if they seem to be OK return false
1471 Info("RunCoherencyCheck","Checking HBT Runs Coherency");
1473 //When we use non-buffering reader this is a big waste of time -> We need to read all data to check it
1474 //and reader is implemented safe in this case anyway
1475 // Info("RunCoherencyCheck","Number of events ...");
1476 // if (fReader->GetNumberOfPartEvents() == fReader->GetNumberOfTrackEvents() ) //check whether there is the same number of events
1478 // Info("RunCoherencyCheck","OK. %d found\n",fReader->GetNumberOfTrackEvents());
1481 // { //if not the same - ERROR
1482 // Error("RunCoherencyCheck",
1483 // "Number of simulated events (%d) is not equal to number of reconstructed events(%d)",
1484 // fReader->GetNumberOfPartEvents(),fReader->GetNumberOfTrackEvents());
1488 Info("RunCoherencyCheck","Checking number of Particles AND Particles Types in each event ...");
1492 for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
1494 partEvent= fReader->GetEventSim(i); //gets the "ith" event
1495 trackEvent = fReader->GetEventRec(i);
1497 if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
1498 if ( (partEvent == 0x0) || (partEvent == 0x0) )
1500 Error("RunCoherencyCheck",
1501 "One event is NULL and the other one not. Event Number %d",i);
1505 if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
1507 Error("RunCoherencyCheck",
1508 "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
1509 i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
1513 for (Int_t j = 0; j<partEvent->GetNumberOfParticles(); j++)
1515 if( partEvent->GetParticle(j)->GetPdgCode() != trackEvent->GetParticle(j)->GetPdgCode() )
1517 Error("RunCoherencyCheck",
1518 "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
1519 i,j, partEvent->GetParticle(j)->GetPdgCode(),trackEvent->GetParticle(j)->GetPdgCode() );
1525 Info("RunCoherencyCheck"," Done");
1526 Info("RunCoherencyCheck"," Everything looks OK");
1531 /*************************************************************************************/
1533 void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
1535 //Performs analysis for both, tracks and particles
1536 AliAOD* rawtrackEvent, * rawpartEvent;
1539 Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
1540 Info("ProcessTracksAndParticlesNonIdentAnal","***** NON IDENT MODE ****************");
1541 Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
1543 for (Int_t i = 0;;i++)//infinite loop
1545 if (fReader->Next()) break; //end when no more events available
1547 rawpartEvent = fReader->GetEventSim();
1548 rawtrackEvent = fReader->GetEventRec();
1550 ProcessRecAndSimNonId(rawtrackEvent,rawpartEvent);
1551 }//end of loop over events (1)
1553 /*************************************************************************************/
1555 void AliHBTAnalysis::ProcessTracksNonIdentAnal()
1557 //Process Tracks only with non identical mode
1558 AliAOD * rawtrackEvent;
1561 Info("ProcessTracksNonIdentAnal","**************************************");
1562 Info("ProcessTracksNonIdentAnal","***** NON IDENT MODE ****************");
1563 Info("ProcessTracksNonIdentAnal","**************************************");
1565 for (Int_t i = 0;;i++)//infinite loop
1567 if (fReader->Next()) break; //end when no more events available
1568 rawtrackEvent = fReader->GetEventRec();
1569 ProcessRecNonId(rawtrackEvent,0x0);
1570 }//end of loop over events (1)
1572 /*************************************************************************************/
1574 void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
1576 //process paricles only with non identical mode
1577 AliAOD * rawpartEvent = 0x0;
1580 Info("ProcessParticlesNonIdentAnal","**************************************");
1581 Info("ProcessParticlesNonIdentAnal","***** NON IDENT MODE ****************");
1582 Info("ProcessParticlesNonIdentAnal","**************************************");
1584 for (Int_t i = 0;;i++)//infinite loop
1586 if (fReader->Next()) break; //end when no more events available
1588 rawpartEvent = fReader->GetEventSim();
1589 ProcessSimNonId(0x0,rawpartEvent);
1590 }//end of loop over events (1)
1593 /*************************************************************************************/
1594 void AliHBTAnalysis::FilterOut(AliAOD* outpart1, AliAOD* outpart2, AliAOD* inpart,
1595 AliAOD* outtrack1, AliAOD* outtrack2, AliAOD* intrack) const
1597 //Puts particles accepted as a first particle by global cut in out1
1598 //and as a second particle in out2
1600 AliVAODParticle* part, *track;
1609 for (Int_t i = 0; i < inpart->GetNumberOfParticles(); i++)
1612 part = inpart->GetParticle(i);
1613 track = intrack->GetParticle(i);
1615 if ( ((this->*fkPass1)(part,track)) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
1616 if ( ((this->*fkPass2)(part,track)) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
1618 if (gDebug)//to be removed in real analysis
1619 if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1621 //Particle accpted by both cuts
1622 Error("FilterOut","Particle accepted by both cuts");
1628 outpart1->AddParticle(part);
1629 outtrack1->AddParticle(track);
1635 outpart2->AddParticle(part);
1636 outtrack2->AddParticle(track);
1641 /*************************************************************************************/
1643 void AliHBTAnalysis::FilterOut(AliAOD* out1, AliAOD* out2, AliAOD* in) const
1645 //Puts particles accepted as a first particle by global cut in out1
1646 //and as a second particle in out2
1647 AliVAODParticle* part;
1652 AliAODParticleCut *cut1 = fPairCut->GetFirstPartCut();
1653 AliAODParticleCut *cut2 = fPairCut->GetSecondPartCut();
1657 for (Int_t i = 0; i < in->GetNumberOfParticles(); i++)
1660 part = in->GetParticle(i);
1662 if ( cut1->Rejected(part) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
1663 if ( cut2->Rejected(part) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
1665 if (gDebug)//to be removed in real analysis
1666 if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1668 //Particle accpted by both cuts
1669 Error("FilterOut","Particle accepted by both cuts");
1675 out1->AddParticle(part);
1681 out2->AddParticle(part);
1686 /*************************************************************************************/
1688 Bool_t AliHBTAnalysis::IsNonIdentAnalysis()
1690 //checks if it is possible to use special analysis for non identical particles
1691 //it means - in global pair cut first particle id is different than second one
1692 //and both are different from 0
1693 //in the future is possible to perform more sophisticated check
1694 //if cuts have excluding requirements
1696 if (fPairCut->IsEmpty())
1699 if (fPairCut->GetFirstPartCut()->IsEmpty())
1702 if (fPairCut->GetSecondPartCut()->IsEmpty())
1705 Int_t id1 = fPairCut->GetFirstPartCut()->GetPID();
1706 Int_t id2 = fPairCut->GetSecondPartCut()->GetPID();
1708 if ( (id1==0) || (id2==0) || (id1==id2) )
1713 /*************************************************************************************/
1715 void AliHBTAnalysis::SetApparentVertex(Double_t x, Double_t y, Double_t z)
1717 //Sets apparent vertex
1718 // All events have to be moved to the same vertex position in order to
1719 // to be able to comare any space positions (f.g. anti-merging)
1720 // This method defines this position
1726 /*************************************************************************************/
1728 void AliHBTAnalysis::PressAnyKey()
1730 //small utility function that helps to make comfortable macros
1733 fcntl(0, F_SETFL, O_NONBLOCK);
1734 ::Info("","Press Any Key to continue ...");
1737 nread = read(0, &c, 1);
1738 gSystem->ProcessEvents();
1742 /*************************************************************************************/