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 //_________________________________________________________
22 #include <Riostream.h>
29 #include "AliAODParticle.h"
30 #include "AliAODPairCut.h"
31 #include "AliEventCut.h"
33 #include "AliEventBuffer.h"
35 #include "AliReader.h"
36 #include "AliHBTPair.h"
37 #include "AliHBTFunction.h"
38 #include "AliHBTMonitorFunction.h"
41 ClassImp(AliHBTAnalysis)
43 const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100;
44 const UInt_t AliHBTAnalysis::fgkDefaultMixingInfo = 1000;
45 const Int_t AliHBTAnalysis::fgkDefaultBufferSize = 5;
47 AliHBTAnalysis::AliHBTAnalysis():
51 fNParticleFunctions(0),
52 fNParticleAndTrackFunctions(0),
53 fNTrackMonitorFunctions(0),
54 fNParticleMonitorFunctions(0),
55 fNParticleAndTrackMonitorFunctions(0),
56 fTrackFunctions ( new AliHBTOnePairFctn* [fgkFctnArraySize]),
57 fParticleFunctions ( new AliHBTOnePairFctn* [fgkFctnArraySize]),
58 fParticleAndTrackFunctions ( new AliHBTTwoPairFctn* [fgkFctnArraySize]),
59 fParticleMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),
60 fTrackMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),
61 fParticleAndTrackMonitorFunctions ( new AliHBTMonTwoParticleFctn* [fgkFctnArraySize]),
66 fDisplayMixingInfo(fgkDefaultMixingInfo),
68 fProcessOption(kSimulatedAndReconstructed),
78 /*************************************************************************************/
80 AliHBTAnalysis::AliHBTAnalysis(const AliHBTAnalysis& in):
85 fNParticleFunctions(0),
86 fNParticleAndTrackFunctions(0),
87 fNTrackMonitorFunctions(0),
88 fNParticleMonitorFunctions(0),
89 fNParticleAndTrackMonitorFunctions(0),
91 fParticleFunctions(0x0),
92 fParticleAndTrackFunctions(0x0),
93 fParticleMonitorFunctions(0x0),
94 fTrackMonitorFunctions(0x0),
95 fParticleAndTrackMonitorFunctions(0x0),
99 fBufferSize(fgkDefaultBufferSize),
100 fDisplayMixingInfo(fgkDefaultMixingInfo),
102 fProcessOption(kSimulatedAndReconstructed),
103 fNoCorrfctns(kFALSE),
104 fOutputFileName(0x0),
110 Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
112 /*************************************************************************************/
113 AliHBTAnalysis& AliHBTAnalysis::operator=(const AliHBTAnalysis& /*right*/)
116 Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
119 /*************************************************************************************/
120 AliHBTAnalysis::~AliHBTAnalysis()
123 //note that we do not delete functions itself
124 // they should be deleted by whom where created
125 //we only store pointers, and we use "new" only for pointers array
129 if (AliVAODParticle::GetDebug()>5)Info("~AliHBTAnalysis","Is Owner: Attempting to delete functions");
131 if (AliVAODParticle::GetDebug()>5)Info("~AliHBTAnalysis","Delete functions done");
133 delete [] fTrackFunctions;
134 delete [] fParticleFunctions;
135 delete [] fParticleAndTrackFunctions;
137 delete [] fParticleMonitorFunctions;
138 delete [] fTrackMonitorFunctions;
139 delete [] fParticleAndTrackMonitorFunctions;
142 delete fOutputFileName;
145 /*************************************************************************************/
147 Int_t AliHBTAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim)
149 //Processes one event
150 if (fProcEvent == 0x0)
152 Error("ProcessEvent","Analysis <<%s>> option not specified.",GetName());
155 if ( Rejected(aodrec,aodsim) )
157 // Double_t x = 0, y = 0, z = 0;
158 // aodrec->GetPrimaryVertex(x,y,z);
159 // Info("ProcessEvent","Event has vertex at %f %f %f",x,y,z);
160 Info("ProcessEvent","Nch is %d",aodsim->GetNumberOfCharged());
161 Info("ProcessEvent","%s: Event cut rejected this event",GetName());
165 //Move event to the apparent vertex -> must be after the event cut
166 //It is important for any cut that use any spacial coordiantes,
167 //f.g. anti-merging cut in order to preserve the same bahavior of variables (f.g. distance between tracks)
168 Double_t dvx = 0, dvy = 0, dvz = 0;
171 Double_t pvx,pvy,pvz;
172 aodrec->GetPrimaryVertex(pvx,pvy,pvz);
174 dvx = fVertexX - pvx;
175 dvy = fVertexY - pvy;
176 dvz = fVertexZ - pvz;
177 aodrec->Move(dvx,dvy,dvz);
180 Int_t result = (this->*fProcEvent)(aodrec,aodsim);
182 if (aodrec) aodrec->Move(-dvx,-dvy,-dvz);//move event back to the same spacial coordinates
186 /*************************************************************************************/
188 Int_t AliHBTAnalysis::Finish()
194 /*************************************************************************************/
196 void AliHBTAnalysis::DeleteFunctions()
198 //Deletes all functions added to analysis
201 for(ii = 0;ii<fNParticleFunctions;ii++)
203 if (AliVAODParticle::GetDebug()>5)
205 Info("DeleteFunctions","Deleting ParticleFunction %#x",fParticleFunctions[ii]);
206 Info("DeleteFunctions","Deleting ParticleFunction %s",fParticleFunctions[ii]->Name());
208 delete fParticleFunctions[ii];
210 fNParticleFunctions = 0;
212 for(ii = 0;ii<fNTrackFunctions;ii++)
214 if (AliVAODParticle::GetDebug()>5)
216 Info("DeleteFunctions","Deleting TrackFunction %#x",fTrackFunctions[ii]);
217 Info("DeleteFunctions","Deleting TrackFunction %s",fTrackFunctions[ii]->Name());
219 delete fTrackFunctions[ii];
221 fNTrackFunctions = 0;
223 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
225 if (AliVAODParticle::GetDebug()>5)
227 Info("DeleteFunctions","Deleting ParticleAndTrackFunction %#x",fParticleAndTrackFunctions[ii]);
228 Info("DeleteFunctions","Deleting ParticleAndTrackFunction %s",fParticleAndTrackFunctions[ii]->Name());
230 delete fParticleAndTrackFunctions[ii];
232 fNParticleAndTrackFunctions = 0;
234 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
236 if (AliVAODParticle::GetDebug()>5)
238 Info("DeleteFunctions","Deleting ParticleMonitorFunction %#x",fParticleMonitorFunctions[ii]);
239 Info("DeleteFunctions","Deleting ParticleMonitorFunction %s",fParticleMonitorFunctions[ii]->Name());
241 delete fParticleMonitorFunctions[ii];
243 fNParticleMonitorFunctions = 0;
245 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
247 if (AliVAODParticle::GetDebug()>5)
249 Info("DeleteFunctions","Deleting TrackMonitorFunction %#x",fTrackMonitorFunctions[ii]);
250 Info("DeleteFunctions","Deleting TrackMonitorFunction %s",fTrackMonitorFunctions[ii]->Name());
252 delete fTrackMonitorFunctions[ii];
254 fNTrackMonitorFunctions = 0;
256 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
258 if (AliVAODParticle::GetDebug()>5)
260 Info("DeleteFunctions","Deleting ParticleAndTrackMonitorFunction %#x",fParticleAndTrackMonitorFunctions[ii]);
261 Info("DeleteFunctions","Deleting ParticleAndTrackMonitorFunction %s",fParticleAndTrackMonitorFunctions[ii]->Name());
263 delete fParticleAndTrackMonitorFunctions[ii];
265 fNParticleAndTrackMonitorFunctions = 0;
268 /*************************************************************************************/
270 Int_t AliHBTAnalysis::Init()
272 //Initializeation method
273 //calls Init for all functions
275 //Depending on option and pair cut assigns proper analysis method
276 Bool_t nonid = IsNonIdentAnalysis();
277 switch(fProcessOption)
280 if (nonid) fProcEvent = &AliHBTAnalysis::ProcessRecNonId;
281 else fProcEvent = &AliHBTAnalysis::ProcessRec;
286 if (nonid) fProcEvent = &AliHBTAnalysis::ProcessSimNonId;
287 else fProcEvent = &AliHBTAnalysis::ProcessSim;
291 case kSimulatedAndReconstructed:
292 if (nonid) fProcEvent = &AliHBTAnalysis::ProcessRecAndSimNonId;
293 else fProcEvent = &AliHBTAnalysis::ProcessRecAndSim;
297 if (fPartBuffer == 0x0) fPartBuffer = new AliEventBuffer (fBufferSize);
298 else fPartBuffer->Reset();
300 if (fTrackBuffer == 0x0) fTrackBuffer = new AliEventBuffer (fBufferSize);
301 else fTrackBuffer->Reset();
304 fNoCorrfctns = (fNParticleFunctions == 0) && (fNTrackFunctions == 0) && (fNParticleAndTrackFunctions == 0);
307 for(ii = 0;ii<fNParticleFunctions;ii++)
308 fParticleFunctions[ii]->Init();
310 for(ii = 0;ii<fNTrackFunctions;ii++)
311 fTrackFunctions[ii]->Init();
313 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
314 fParticleAndTrackFunctions[ii]->Init();
316 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
317 fParticleMonitorFunctions[ii]->Init();
319 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
320 fTrackMonitorFunctions[ii]->Init();
322 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
323 fParticleAndTrackMonitorFunctions[ii]->Init();
327 /*************************************************************************************/
329 void AliHBTAnalysis::ResetFunctions()
331 //In case fOwner is true, deletes all functions
332 //in other case, just set number of analysis to 0
333 if (fIsOwner) DeleteFunctions();
336 fNParticleFunctions = 0;
337 fNTrackFunctions = 0;
338 fNParticleAndTrackFunctions = 0;
339 fNParticleMonitorFunctions = 0;
340 fNTrackMonitorFunctions = 0;
341 fNParticleAndTrackMonitorFunctions = 0;
344 /*************************************************************************************/
345 Int_t AliHBTAnalysis::ProcessRecAndSim(AliAOD* aodrec, AliAOD* aodsim)
347 //Does analysis for both tracks and particles
348 //mainly for resolution study and analysies with weighting algirithms
350 // cut on particles only -- why?
351 // - PID: when we make resolution analysis we want to take only tracks with correct PID
352 // We need cut on tracks because there are data characteristic
354 AliVAODParticle * part1, * part2;
355 AliVAODParticle * track1, * track2;
357 AliAOD * trackEvent = aodrec, *partEvent = aodsim;
358 AliAOD* trackEvent1 = new AliAOD();
359 AliAOD* partEvent1 = new AliAOD();
361 AliAOD * trackEvent2,*partEvent2;
363 // Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
365 // Int_t nev = fReader->GetNumberOfTrackEvents();
366 static AliHBTPair tpair;
367 static AliHBTPair ppair;
369 AliHBTPair* trackpair = &tpair;
370 AliHBTPair* partpair = &ppair;
372 AliHBTPair * tmptrackpair;//temprary pointers to pairs
373 AliHBTPair * tmppartpair;
379 if ( !partEvent || !trackEvent )
381 Error("ProcessRecAndSim","<<%s>> Can not get event",GetName());
385 if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
387 Error("ProcessRecAndSim",
388 "Number of simulated particles (%d) not equal to number of reconstructed tracks (%d). Skipping Event.",
389 partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
394 for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
396 /***************************************/
397 /****** Looping same events ********/
398 /****** filling numerators ********/
399 /***************************************/
400 if ( (j%fDisplayMixingInfo) == 0)
401 Info("ProcessRecAndSim",
402 "Mixing particle %d with particles from the same event",j);
404 part1= partEvent->GetParticle(j);
405 track1= trackEvent->GetParticle(j);
407 Bool_t firstcut = (this->*fkPass1)(part1,track1);
408 if (fBufferSize != 0)
409 if ( (firstcut == kFALSE) || ( (this->*fkPass2)(part1,track1) == kFALSE ) )
411 //accepted by any cut
412 // we have to copy because reader keeps only one event
414 partEvent1->AddParticle(part1);
415 trackEvent1->AddParticle(track1);
418 if (firstcut) continue;
420 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
421 fParticleMonitorFunctions[ii]->Process(part1);
422 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
423 fTrackMonitorFunctions[ii]->Process(track1);
424 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
425 fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
427 if (fNoCorrfctns) continue;
429 for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
431 part2= partEvent->GetParticle(k);
432 if (part1->GetUID() == part2->GetUID()) continue;
433 partpair->SetParticles(part1,part2);
435 track2= trackEvent->GetParticle(k);
436 trackpair->SetParticles(track1,track2);
438 if( (this->*fkPass)(partpair,trackpair) ) //check pair cut
439 { //do not meets crietria of the pair cut, try with swapped pairs
440 if( (this->*fkPass)((AliHBTPair*)partpair->GetSwappedPair(),(AliHBTPair*)trackpair->GetSwappedPair()) )
441 continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
443 { //swaped pair meets all the criteria
444 tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
445 tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
449 {//meets criteria of the pair cut
450 tmptrackpair = trackpair;
451 tmppartpair = partpair;
454 for(ii = 0;ii<fNParticleFunctions;ii++)
455 fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
457 for(ii = 0;ii<fNTrackFunctions;ii++)
458 fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
460 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
461 fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair,tmppartpair);
462 //end of 2nd loop over particles from the same event
463 }//for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
465 /***************************************/
466 /***** Filling denominators *********/
467 /***************************************/
468 if (fBufferSize == 0) continue;
470 fPartBuffer->ResetIter();
471 fTrackBuffer->ResetIter();
473 while (( partEvent2 = fPartBuffer->Next() ))
475 trackEvent2 = fTrackBuffer->Next();
478 if ( (j%fDisplayMixingInfo) == 0)
479 Info("ProcessRecAndSim",
480 "Mixing particle %d from current event with particles from event %d",j,-m);
482 for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
484 part2= partEvent2->GetParticle(l);
485 partpair->SetParticles(part1,part2);
487 track2= trackEvent2->GetParticle(l);
488 trackpair->SetParticles(track1,track2);
490 if( (this->*fkPass)(partpair,trackpair) ) //check pair cut
491 { //do not meets crietria of the
492 if( (this->*fkPass)((AliHBTPair*)partpair->GetSwappedPair(),(AliHBTPair*)trackpair->GetSwappedPair()) )
496 tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
497 tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
501 {//meets criteria of the pair cut
502 tmptrackpair = trackpair;
503 tmppartpair = partpair;
506 for(ii = 0;ii<fNParticleFunctions;ii++)
507 fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
509 for(ii = 0;ii<fNTrackFunctions;ii++)
510 fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
512 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
513 fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair);
514 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
517 //end of loop over particles from first event
518 }//for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
519 delete fPartBuffer->Push(partEvent1);
520 delete fTrackBuffer->Push(trackEvent1);
521 //end of loop over events
524 /*************************************************************************************/
526 Int_t AliHBTAnalysis::ProcessSim(AliAOD* /*aodrec*/, AliAOD* aodsim)
528 //Does analysis of simulated data
529 AliVAODParticle * part1, * part2;
531 AliAOD* partEvent = aodsim;
532 AliAOD* partEvent1 = new AliAOD();
538 AliHBTPair* partpair = &ppair;
540 AliHBTPair * tmppartpair;
546 Error("ProcessSim","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)
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)
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
644 delete fPartBuffer->Push(partEvent1);
645 //end of loop over events
648 /*************************************************************************************/
649 Int_t AliHBTAnalysis::ProcessRec(AliAOD* aodrec, AliAOD* /*aodsim*/)
651 //Does analysis of reconstructed data
652 AliVAODParticle * track1, * track2;
654 AliAOD* trackEvent = aodrec;
655 AliAOD* trackEvent1 = new AliAOD();
661 AliHBTPair* trackpair = &tpair;
663 AliHBTPair * tmptrackpair;
669 Error("ProcessRecAndSim","Can not get event");
674 for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
676 /***************************************/
677 /****** Looping same events ********/
678 /****** filling numerators ********/
679 /***************************************/
680 if ( (j%fDisplayMixingInfo) == 0)
682 "Mixing Particle %d with Particles from the same event",j);
684 track1= trackEvent->GetParticle(j);
686 Bool_t firstcut = fPairCut->GetFirstPartCut()->Rejected(track1);
688 if (fBufferSize != 0)
689 if ( (firstcut == kFALSE) || ( fPairCut->GetSecondPartCut()->Rejected(track1) == kFALSE ) )
691 //accepted by any cut
692 // we have to copy because reader keeps only one event
694 trackEvent1->AddParticle(track1);
697 if (firstcut) continue;
699 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
700 fParticleMonitorFunctions[ii]->Process(track1);
702 if ( fNTrackFunctions == 0 ) continue;
704 for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
706 track2= trackEvent->GetParticle(k);
707 if (track1->GetUID() == track2->GetUID()) continue;
708 trackpair->SetParticles(track1,track2);
710 if(fPairCut->Rejected(trackpair)) //check pair cut
711 { //do not meets crietria of the
712 if( fPairCut->Rejected((AliHBTPair*)trackpair->GetSwappedPair()) ) continue;
713 else tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
717 tmptrackpair = trackpair;
720 for(ii = 0;ii<fNTrackFunctions;ii++)
721 fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
723 //end of 2nd loop over Particles from the same event
724 }//for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
726 /***************************************/
727 /***** Filling denominators *********/
728 /***************************************/
729 if (fBufferSize == 0) continue;
731 fTrackBuffer->ResetIter();
733 while (( trackEvent2 = fTrackBuffer->Next() ))
736 if ( (j%fDisplayMixingInfo) == 0)
738 "Mixing Particle %d from current event with Particles from event %d",j,-m);
739 for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++) // ... on all Particles
742 track2= trackEvent2->GetParticle(l);
743 trackpair->SetParticles(track1,track2);
745 if( fPairCut->Rejected(trackpair) ) //check pair cut
746 { //do not meets crietria of the
747 if( fPairCut->Rejected((AliHBTPair*)trackpair->GetSwappedPair()) )
751 tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
755 {//meets criteria of the pair cut
756 tmptrackpair = trackpair;
759 for(ii = 0;ii<fNTrackFunctions;ii++)
760 fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
762 }//for(Int_t l = 0; l<N2;l++) // ... on all Particles
765 delete fTrackBuffer->Push(trackEvent1);
766 //end of loop over events
769 /*************************************************************************************/
771 Int_t AliHBTAnalysis::ProcessRecAndSimNonId(AliAOD* aodrec, AliAOD* aodsim)
773 //Analyzes both reconstructed and simulated data
776 Error("ProcessTracksAndParticlesNonIdentAnal","Reconstructed event is NULL");
782 Error("ProcessTracksAndParticlesNonIdentAnal","Simulated event is NULL");
786 if ( aodrec->GetNumberOfParticles() != aodsim->GetNumberOfParticles() )
788 Error("ProcessTracksAndParticlesNonIdentAnal",
789 "Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
790 aodsim->GetNumberOfParticles() , aodrec->GetNumberOfParticles());
795 AliVAODParticle * part1, * part2;
796 AliVAODParticle * track1, * track2;
798 static AliAOD aodrec1;
799 static AliAOD aodsim1;
801 AliAOD * trackEvent1=&aodrec1,*partEvent1=&aodsim1;//Particle that passes first particle cut, this event
802 trackEvent1->Reset();
804 AliAOD * trackEvent2=0x0,*partEvent2=0x0;//Particle that passes second particle cut, this event
805 AliAOD * trackEvent3=0x0,*partEvent3=0x0;//Particle that passes second particle cut, events from buffer
807 AliAOD* rawtrackEvent = aodrec;//this we get from Reader
808 AliAOD* rawpartEvent = aodsim;//this we get from Reader
810 static AliHBTPair tpair;
811 static AliHBTPair ppair;
813 AliHBTPair* trackpair = &tpair;
814 AliHBTPair* partpair = &ppair;
818 /********************************/
820 /********************************/
821 if ( ( (partEvent2==0x0) || (trackEvent2==0x0)) )//in case fBufferSize == 0 and pointers are created do not eneter
823 partEvent2 = new AliAOD();
824 trackEvent2 = new AliAOD();
827 FilterOut(partEvent1, partEvent2, rawpartEvent, trackEvent1, trackEvent2, rawtrackEvent);
829 for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
831 if ( (j%fDisplayMixingInfo) == 0)
832 Info("ProcessTracksAndParticlesNonIdentAnal",
833 "Mixing particle %d from current event with particles from current event",j);
835 part1= partEvent1->GetParticle(j);
836 track1= trackEvent1->GetParticle(j);
839 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
840 fParticleMonitorFunctions[ii]->Process(part1);
841 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
842 fTrackMonitorFunctions[ii]->Process(track1);
843 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
844 fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
846 if (fNoCorrfctns) continue;
848 /***************************************/
849 /****** filling numerators ********/
850 /****** (particles from event2) ********/
851 /***************************************/
853 for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups
855 part2= partEvent2->GetParticle(k);
856 if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
857 partpair->SetParticles(part1,part2);
859 track2= trackEvent2->GetParticle(k);
860 trackpair->SetParticles(track1,track2);
862 if( (this->*fkPassPairProp)(partpair,trackpair) ) //check pair cut
863 { //do not meets crietria of the pair cut
867 {//meets criteria of the pair cut
868 for(ii = 0;ii<fNParticleFunctions;ii++)
869 fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
871 for(ii = 0;ii<fNTrackFunctions;ii++)
872 fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
874 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
875 fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(trackpair,partpair);
879 if ( fBufferSize == 0) continue;//do not mix diff histograms
880 /***************************************/
881 /***** Filling denominators *********/
882 /***************************************/
883 fPartBuffer->ResetIter();
884 fTrackBuffer->ResetIter();
888 while ( (partEvent3 = fPartBuffer->Next() ) != 0x0)
890 trackEvent3 = fTrackBuffer->Next();
892 if ( (j%fDisplayMixingInfo) == 0)
893 Info("ProcessTracksAndParticlesNonIdentAnal",
894 "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
896 for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
898 part2= partEvent3->GetParticle(k);
899 partpair->SetParticles(part1,part2);
901 track2= trackEvent3->GetParticle(k);
902 trackpair->SetParticles(track1,track2);
904 if( (this->*fkPassPairProp)(partpair,trackpair) ) //check pair cut
905 { //do not meets crietria of the pair cut
909 {//meets criteria of the pair cut
911 for(ii = 0;ii<fNParticleFunctions;ii++)
912 fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
914 for(ii = 0;ii<fNTrackFunctions;ii++)
915 fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
917 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
918 fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(trackpair,partpair);
920 }// for particles event2
922 }//for over particles in event1
924 delete fPartBuffer->Push(partEvent2);
925 delete fTrackBuffer->Push(trackEvent2);
929 /*************************************************************************************/
930 Int_t AliHBTAnalysis::ProcessSimNonId(AliAOD* /*aodrec*/, AliAOD* aodsim)
932 //does analysis of simulated (MC) data in non-identical mode
933 //i.e. when particles selected by first part. cut are a disjunctive set than particles
934 //passed by the second part. cut
941 AliVAODParticle * part1, * part2;
943 static AliAOD aodsim1;
945 AliAOD* partEvent1=&aodsim1;//Particle that passes first particle cut, this event
947 AliAOD* partEvent2=0x0;//Particle that passes second particle cut, this event
948 AliAOD* partEvent3=0x0;//Particle that passes second particle cut, events from buffer
950 AliAOD* rawpartEvent = aodsim;//this we get from Reader
952 static AliHBTPair ppair;
954 AliHBTPair* partpair = &ppair;
958 /********************************/
960 /********************************/
961 if (partEvent2==0x0)//in case fBufferSize == 0 and pointers are created do not eneter
963 partEvent2 = new AliAOD();
966 FilterOut(partEvent1, partEvent2, rawpartEvent);
968 for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
970 if ( (j%fDisplayMixingInfo) == 0)
971 Info("ProcessParticlesNonIdentAnal",
972 "Mixing particle %d from current event with particles from current event",j);
974 part1= partEvent1->GetParticle(j);
977 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
978 fParticleMonitorFunctions[ii]->Process(part1);
980 if (fNParticleFunctions == 0) continue;
982 /***************************************/
983 /****** filling numerators ********/
984 /****** (particles from event2) ********/
985 /***************************************/
987 for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups
989 part2= partEvent2->GetParticle(k);
990 if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
991 partpair->SetParticles(part1,part2);
994 if(fPairCut->PassPairProp(partpair) ) //check pair cut
995 { //do not meets crietria of the pair cut
999 {//meets criteria of the pair cut
1000 for(ii = 0;ii<fNParticleFunctions;ii++)
1001 fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
1005 if ( fBufferSize == 0) continue;//do not mix diff histograms
1006 /***************************************/
1007 /***** Filling denominators *********/
1008 /***************************************/
1009 fPartBuffer->ResetIter();
1013 while ( (partEvent3 = fPartBuffer->Next() ) != 0x0)
1016 if ( (j%fDisplayMixingInfo) == 0)
1017 Info("ProcessParticlesNonIdentAnal",
1018 "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
1020 for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
1022 part2= partEvent3->GetParticle(k);
1023 partpair->SetParticles(part1,part2);
1026 if(fPairCut->PassPairProp(partpair) ) //check pair cut
1027 { //do not meets crietria of the pair cut
1031 {//meets criteria of the pair cut
1032 for(ii = 0;ii<fNParticleFunctions;ii++)
1034 fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
1037 }// for particles event2
1039 }//for over particles in event1
1041 delete fPartBuffer->Push(partEvent2);
1045 /*************************************************************************************/
1046 Int_t AliHBTAnalysis::ProcessRecNonId(AliAOD* aodrec, AliAOD* /*aodsim*/)
1048 //Analyzes both reconstructed and simulated data
1054 AliVAODParticle * track1, * track2;
1056 static AliAOD aodrec1;
1057 AliAOD * trackEvent1=&aodrec1;//Particle that passes first particle cut, this event
1058 trackEvent1->Reset();
1059 AliAOD * trackEvent2=0x0;//Particle that passes second particle cut, this event
1060 AliAOD * trackEvent3=0x0;//Particle that passes second particle cut, events from buffer
1061 AliAOD* rawtrackEvent = aodrec;//this we get from Reader
1063 static AliHBTPair tpair;
1065 AliHBTPair* trackpair = &tpair;
1070 /********************************/
1072 /********************************/
1073 if ( trackEvent2==0x0 )//in case fBufferSize == 0 and pointers are created do not eneter
1075 trackEvent2 = new AliAOD();
1078 FilterOut(trackEvent1, trackEvent2, rawtrackEvent);
1080 for (Int_t j = 0; j<trackEvent1->GetNumberOfParticles() ; j++)
1082 if ( (j%fDisplayMixingInfo) == 0)
1083 Info("ProcessTracksNonIdentAnal",
1084 "Mixing particle %d from current event with particles from current event",j);
1086 track1= trackEvent1->GetParticle(j);
1089 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
1090 fTrackMonitorFunctions[ii]->Process(track1);
1092 if (fNTrackFunctions == 0x0) continue;
1094 /***************************************/
1095 /****** filling numerators ********/
1096 /****** (particles from event2) ********/
1097 /***************************************/
1099 for (Int_t k = 0; k < trackEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups
1101 track2= trackEvent2->GetParticle(k);
1102 if (track1->GetUID() == track2->GetUID()) continue;//this is the same particle but with different PID
1103 trackpair->SetParticles(track1,track2);
1106 if( fPairCut->PassPairProp(trackpair)) //check pair cut
1107 { //do not meets crietria of the pair cut
1111 {//meets criteria of the pair cut
1113 for(ii = 0;ii<fNTrackFunctions;ii++)
1114 fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
1118 if ( fBufferSize == 0) continue;//do not mix diff histograms
1119 /***************************************/
1120 /***** Filling denominators *********/
1121 /***************************************/
1122 fTrackBuffer->ResetIter();
1126 while ( (trackEvent3 = fTrackBuffer->Next() ) != 0x0)
1128 if ( (j%fDisplayMixingInfo) == 0)
1129 Info("ProcessTracksNonIdentAnal",
1130 "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
1132 for (Int_t k = 0; k < trackEvent3->GetNumberOfParticles() ; k++)
1134 track2= trackEvent3->GetParticle(k);
1135 trackpair->SetParticles(track1,track2);
1137 if( fPairCut->PassPairProp(trackpair)) //check pair cut
1138 { //do not meets crietria of the pair cut
1142 {//meets criteria of the pair cut
1143 for(ii = 0;ii<fNTrackFunctions;ii++)
1144 fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
1146 }// for particles event2
1148 }//for over particles in event1
1150 delete fTrackBuffer->Push(trackEvent2);
1154 /*************************************************************************************/
1156 void AliHBTAnalysis::Process(Option_t* option)
1158 //default option = "TracksAndParticles"
1159 //Main method of the HBT Analysis Package
1160 //It triggers reading with the global cut (default is an empty cut)
1161 //Than it checks options and data which are read
1162 //if everything is OK, then it calls one of the looping methods
1163 //depending on tfReaderhe option
1164 //These methods differs on what thay are looping on
1167 //--------------------------------------------------------------------
1168 //ProcessTracksAndParticles - "TracksAndParticles"
1170 // it loops over both, tracks(reconstructed) and particles(simulated)
1171 // all function gethered in all 3 lists are called for each (double)pair
1173 //ProcessTracks - "Tracks"
1174 // it loops only on tracks(reconstructed),
1175 // functions ONLY from fTrackFunctions list are called
1177 //ProcessParticles - "Particles"
1178 // it loops only on particles(simulated),
1179 // functions ONLY from fParticleAndTrackFunctions list are called
1184 Error("Process","The reader is not set");
1188 const char *oT = strstr(option,"Tracks");
1189 const char *oP = strstr(option,"Particles");
1191 Bool_t nonid = IsNonIdentAnalysis();
1197 if (nonid) ProcessTracksAndParticlesNonIdentAnal();
1198 else ProcessTracksAndParticles();
1204 if (nonid) ProcessTracksNonIdentAnal();
1205 else ProcessTracks();
1211 if (nonid) ProcessParticlesNonIdentAnal();
1212 else ProcessParticles();
1217 /*************************************************************************************/
1219 void AliHBTAnalysis::ProcessTracksAndParticles()
1221 //Makes analysis for both tracks and particles
1222 //mainly for resolution study and analysies with weighting algirithms
1223 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1224 //the loops are splited
1226 // cut on particles only -- why?
1227 // - PID: when we make resolution analysis we want to take only tracks with correct PID
1228 // We need cut on tracks because there are data characteristic to
1230 AliAOD * trackEvent, *partEvent;
1233 while (fReader->Next() == kFALSE)
1235 partEvent = fReader->GetEventSim();
1236 trackEvent = fReader->GetEventRec();
1237 ProcessRecAndSim(trackEvent,partEvent);
1239 }//while (fReader->Next() == kFALSE)
1242 /*************************************************************************************/
1244 void AliHBTAnalysis::ProcessTracks()
1246 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1247 //the loops are splited
1248 AliAOD * trackEvent;
1250 while (fReader->Next() == kFALSE)
1252 trackEvent = fReader->GetEventRec();
1253 ProcessRec(trackEvent,0x0);
1254 }//while (fReader->Next() == kFALSE)
1257 /*************************************************************************************/
1259 void AliHBTAnalysis::ProcessParticles()
1261 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1262 //the loops are splited
1265 while (fReader->Next() == kFALSE)
1267 partEvent = fReader->GetEventSim();
1268 ProcessSim(0x0,partEvent);
1269 }//while (fReader->Next() == kFALSE)
1271 /*************************************************************************************/
1273 void AliHBTAnalysis::WriteFunctions()
1275 //Calls Write for all defined functions in analysis
1276 //== writes all results
1277 TFile* oututfile = 0x0;
1278 if (fOutputFileName)
1280 oututfile = TFile::Open(*fOutputFileName,"update");
1283 for(ii = 0;ii<fNParticleFunctions;ii++)
1285 if (AliVAODParticle::GetDebug()>5)
1287 Info("WriteFunctions","Writing ParticleFunction %#x",fParticleFunctions[ii]);
1288 Info("WriteFunctions","Writing ParticleFunction %s",fParticleFunctions[ii]->Name());
1290 fParticleFunctions[ii]->Write();
1293 for(ii = 0;ii<fNTrackFunctions;ii++)
1295 if (AliVAODParticle::GetDebug()>5)
1297 Info("WriteFunctions","Writing TrackFunction %#x",fTrackFunctions[ii]);
1298 Info("WriteFunctions","Writing TrackFunction %s",fTrackFunctions[ii]->Name());
1300 fTrackFunctions[ii]->Write();
1303 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
1305 if (AliVAODParticle::GetDebug()>5)
1307 Info("WriteFunctions","Writing ParticleAndTrackFunction %#x",fParticleAndTrackFunctions[ii]);
1308 Info("WriteFunctions","Writing ParticleAndTrackFunction %s",fParticleAndTrackFunctions[ii]->Name());
1310 fParticleAndTrackFunctions[ii]->Write();
1313 for(ii = 0;ii<fNParticleMonitorFunctions;ii++)
1315 if (AliVAODParticle::GetDebug()>5)
1317 Info("WriteFunctions","Writing ParticleMonitorFunction %#x",fParticleMonitorFunctions[ii]);
1318 Info("WriteFunctions","Writing ParticleMonitorFunction %s",fParticleMonitorFunctions[ii]->Name());
1320 fParticleMonitorFunctions[ii]->Write();
1323 for(ii = 0;ii<fNTrackMonitorFunctions;ii++)
1325 if (AliVAODParticle::GetDebug()>5)
1327 Info("WriteFunctions","Writing TrackMonitorFunction %#x",fTrackMonitorFunctions[ii]);
1328 Info("WriteFunctions","Writing TrackMonitorFunction %s",fTrackMonitorFunctions[ii]->Name());
1330 fTrackMonitorFunctions[ii]->Write();
1333 for(ii = 0;ii<fNParticleAndTrackMonitorFunctions;ii++)
1335 if (AliVAODParticle::GetDebug()>5)
1337 Info("WriteFunctions","Writing ParticleAndTrackMonitorFunction %#x",fParticleAndTrackMonitorFunctions[ii]);
1338 Info("WriteFunctions","Writing ParticleAndTrackMonitorFunction %s",fParticleAndTrackMonitorFunctions[ii]->Name());
1340 fParticleAndTrackMonitorFunctions[ii]->Write();
1344 /*************************************************************************************/
1346 void AliHBTAnalysis::SetOutputFileName(const char* fname)
1348 //Sets fiele name where to dump results,
1349 //if not specified reults are written to gDirectory
1352 delete fOutputFileName;
1353 fOutputFileName = 0x0;
1356 if ( strcmp(fname,"") == 0 )
1358 delete fOutputFileName;
1359 fOutputFileName = 0x0;
1362 if (fOutputFileName == 0x0) fOutputFileName = new TString(fname);
1363 else *fOutputFileName = fname;
1365 /*************************************************************************************/
1367 void AliHBTAnalysis::SetGlobalPairCut(AliAODPairCut* cut)
1369 //Sets the global cut
1372 Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
1375 fPairCut = (AliAODPairCut*)cut->Clone();
1378 /*************************************************************************************/
1380 void AliHBTAnalysis::AddTrackFunction(AliHBTOnePairFctn* f)
1382 //Adds track function
1383 if (f == 0x0) return;
1384 if (fNTrackFunctions == fgkFctnArraySize)
1386 Error("AliHBTAnalysis::AddTrackFunction","Can not add this function, not enough place in the array.");
1388 fTrackFunctions[fNTrackFunctions] = f;
1391 /*************************************************************************************/
1393 void AliHBTAnalysis::AddParticleFunction(AliHBTOnePairFctn* f)
1395 //adds particle function
1396 if (f == 0x0) return;
1398 if (fNParticleFunctions == fgkFctnArraySize)
1400 Error("AliHBTAnalysis::AddParticleFunction","Can not add this function, not enough place in the array.");
1402 fParticleFunctions[fNParticleFunctions] = f;
1403 fNParticleFunctions++;
1405 /*************************************************************************************/
1407 void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTTwoPairFctn* f)
1409 //add resolution function
1410 if (f == 0x0) return;
1411 if (fNParticleAndTrackFunctions == fgkFctnArraySize)
1413 Error("AliHBTAnalysis::AddParticleAndTrackFunction","Can not add this function, not enough place in the array.");
1415 fParticleAndTrackFunctions[fNParticleAndTrackFunctions] = f;
1416 fNParticleAndTrackFunctions++;
1418 /*************************************************************************************/
1420 void AliHBTAnalysis::AddParticleMonitorFunction(AliHBTMonOneParticleFctn* f)
1422 //add particle monitoring function
1423 if (f == 0x0) return;
1425 if (fNParticleMonitorFunctions == fgkFctnArraySize)
1427 Error("AliHBTAnalysis::AddParticleMonitorFunction","Can not add this function, not enough place in the array.");
1429 fParticleMonitorFunctions[fNParticleMonitorFunctions] = f;
1430 fNParticleMonitorFunctions++;
1432 /*************************************************************************************/
1434 void AliHBTAnalysis::AddTrackMonitorFunction(AliHBTMonOneParticleFctn* f)
1436 //add track monitoring function
1437 if (f == 0x0) return;
1439 if (fNTrackMonitorFunctions == fgkFctnArraySize)
1441 Error("AliHBTAnalysis::AddTrackMonitorFunction","Can not add this function, not enough place in the array.");
1443 fTrackMonitorFunctions[fNTrackMonitorFunctions] = f;
1444 fNTrackMonitorFunctions++;
1446 /*************************************************************************************/
1448 void AliHBTAnalysis::AddParticleAndTrackMonitorFunction(AliHBTMonTwoParticleFctn* f)
1450 //add resolution monitoring function
1451 if (f == 0x0) return;
1452 if (fNParticleAndTrackMonitorFunctions == fgkFctnArraySize)
1454 Error("AliHBTAnalysis::AddParticleAndTrackMonitorFunction","Can not add this function, not enough place in the array.");
1456 fParticleAndTrackMonitorFunctions[fNParticleAndTrackMonitorFunctions] = f;
1457 fNParticleAndTrackMonitorFunctions++;
1461 /*************************************************************************************/
1462 /*************************************************************************************/
1464 Bool_t AliHBTAnalysis::RunCoherencyCheck()
1466 //Checks if both HBTRuns are similar
1467 //return true if error found
1468 //if they seem to be OK return false
1472 Info("RunCoherencyCheck","Checking HBT Runs Coherency");
1474 //When we use non-buffering reader this is a big waste of time -> We need to read all data to check it
1475 //and reader is implemented safe in this case anyway
1476 // Info("RunCoherencyCheck","Number of events ...");
1477 // if (fReader->GetNumberOfPartEvents() == fReader->GetNumberOfTrackEvents() ) //check whether there is the same number of events
1479 // Info("RunCoherencyCheck","OK. %d found\n",fReader->GetNumberOfTrackEvents());
1482 // { //if not the same - ERROR
1483 // Error("RunCoherencyCheck",
1484 // "Number of simulated events (%d) is not equal to number of reconstructed events(%d)",
1485 // fReader->GetNumberOfPartEvents(),fReader->GetNumberOfTrackEvents());
1489 Info("RunCoherencyCheck","Checking number of Particles AND Particles Types in each event ...");
1493 for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
1495 partEvent= fReader->GetEventSim(i); //gets the "ith" event
1496 trackEvent = fReader->GetEventRec(i);
1498 if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
1499 if ( (partEvent == 0x0) || (partEvent == 0x0) )
1501 Error("RunCoherencyCheck",
1502 "One event is NULL and the other one not. Event Number %d",i);
1506 if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
1508 Error("RunCoherencyCheck",
1509 "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
1510 i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
1514 for (Int_t j = 0; j<partEvent->GetNumberOfParticles(); j++)
1516 if( partEvent->GetParticle(j)->GetPdgCode() != trackEvent->GetParticle(j)->GetPdgCode() )
1518 Error("RunCoherencyCheck",
1519 "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
1520 i,j, partEvent->GetParticle(j)->GetPdgCode(),trackEvent->GetParticle(j)->GetPdgCode() );
1526 Info("RunCoherencyCheck"," Done");
1527 Info("RunCoherencyCheck"," Everything looks OK");
1532 /*************************************************************************************/
1534 void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
1536 //Performs analysis for both, tracks and particles
1537 AliAOD* rawtrackEvent, * rawpartEvent;
1540 Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
1541 Info("ProcessTracksAndParticlesNonIdentAnal","***** NON IDENT MODE ****************");
1542 Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
1544 for (Int_t i = 0;;i++)//infinite loop
1546 if (fReader->Next()) break; //end when no more events available
1548 rawpartEvent = fReader->GetEventSim();
1549 rawtrackEvent = fReader->GetEventRec();
1551 ProcessRecAndSimNonId(rawtrackEvent,rawpartEvent);
1552 }//end of loop over events (1)
1554 /*************************************************************************************/
1556 void AliHBTAnalysis::ProcessTracksNonIdentAnal()
1558 //Process Tracks only with non identical mode
1559 AliAOD * rawtrackEvent;
1562 Info("ProcessTracksNonIdentAnal","**************************************");
1563 Info("ProcessTracksNonIdentAnal","***** NON IDENT MODE ****************");
1564 Info("ProcessTracksNonIdentAnal","**************************************");
1566 for (Int_t i = 0;;i++)//infinite loop
1568 if (fReader->Next()) break; //end when no more events available
1569 rawtrackEvent = fReader->GetEventRec();
1570 ProcessRecNonId(rawtrackEvent,0x0);
1571 }//end of loop over events (1)
1573 /*************************************************************************************/
1575 void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
1577 //process paricles only with non identical mode
1578 AliAOD * rawpartEvent = 0x0;
1581 Info("ProcessParticlesNonIdentAnal","**************************************");
1582 Info("ProcessParticlesNonIdentAnal","***** NON IDENT MODE ****************");
1583 Info("ProcessParticlesNonIdentAnal","**************************************");
1585 for (Int_t i = 0;;i++)//infinite loop
1587 if (fReader->Next()) break; //end when no more events available
1589 rawpartEvent = fReader->GetEventSim();
1590 ProcessSimNonId(0x0,rawpartEvent);
1591 }//end of loop over events (1)
1594 /*************************************************************************************/
1595 void AliHBTAnalysis::FilterOut(AliAOD* outpart1, AliAOD* outpart2, AliAOD* inpart,
1596 AliAOD* outtrack1, AliAOD* outtrack2, AliAOD* intrack) const
1598 //Puts particles accepted as a first particle by global cut in out1
1599 //and as a second particle in out2
1601 AliVAODParticle* part, *track;
1610 for (Int_t i = 0; i < inpart->GetNumberOfParticles(); i++)
1613 part = inpart->GetParticle(i);
1614 track = intrack->GetParticle(i);
1616 if ( ((this->*fkPass1)(part,track)) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
1617 if ( ((this->*fkPass2)(part,track)) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
1619 if (gDebug)//to be removed in real analysis
1620 if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1622 //Particle accpted by both cuts
1623 Error("FilterOut","Particle accepted by both cuts");
1629 outpart1->AddParticle(part);
1630 outtrack1->AddParticle(track);
1636 outpart2->AddParticle(part);
1637 outtrack2->AddParticle(track);
1642 /*************************************************************************************/
1644 void AliHBTAnalysis::FilterOut(AliAOD* out1, AliAOD* out2, AliAOD* in) const
1646 //Puts particles accepted as a first particle by global cut in out1
1647 //and as a second particle in out2
1648 AliVAODParticle* part;
1653 AliAODParticleCut *cut1 = fPairCut->GetFirstPartCut();
1654 AliAODParticleCut *cut2 = fPairCut->GetSecondPartCut();
1658 for (Int_t i = 0; i < in->GetNumberOfParticles(); i++)
1661 part = in->GetParticle(i);
1663 if ( cut1->Rejected(part) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
1664 if ( cut2->Rejected(part) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
1666 if (gDebug)//to be removed in real analysis
1667 if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1669 //Particle accpted by both cuts
1670 Error("FilterOut","Particle accepted by both cuts");
1676 out1->AddParticle(part);
1682 out2->AddParticle(part);
1687 /*************************************************************************************/
1689 Bool_t AliHBTAnalysis::IsNonIdentAnalysis()
1691 //checks if it is possible to use special analysis for non identical particles
1692 //it means - in global pair cut first particle id is different than second one
1693 //and both are different from 0
1694 //in the future is possible to perform more sophisticated check
1695 //if cuts have excluding requirements
1697 if (fPairCut->IsEmpty())
1700 if (fPairCut->GetFirstPartCut()->IsEmpty())
1703 if (fPairCut->GetSecondPartCut()->IsEmpty())
1706 Int_t id1 = fPairCut->GetFirstPartCut()->GetPID();
1707 Int_t id2 = fPairCut->GetSecondPartCut()->GetPID();
1709 if ( (id1==0) || (id2==0) || (id1==id2) )
1714 /*************************************************************************************/
1716 void AliHBTAnalysis::SetApparentVertex(Double_t x, Double_t y, Double_t z)
1718 //Sets apparent vertex
1719 // All events have to be moved to the same vertex position in order to
1720 // to be able to comare any space positions (f.g. anti-merging)
1721 // This method defines this position
1727 /*************************************************************************************/
1729 void AliHBTAnalysis::PressAnyKey()
1731 //small utility function that helps to make comfortable macros
1734 fcntl(0, F_SETFL, O_NONBLOCK);
1735 ::Info("","Press Any Key to continue ...");
1738 nread = read(0, &c, 1);
1739 gSystem->ProcessEvents();