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://aliweb.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"
42 ClassImp(AliHBTAnalysis)
44 const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100;
45 const UInt_t AliHBTAnalysis::fgkDefaultMixingInfo = 1000;
46 const Int_t AliHBTAnalysis::fgkDefaultBufferSize = 5;
48 AliHBTAnalysis::AliHBTAnalysis():
52 fNParticleFunctions(0),
53 fNParticleAndTrackFunctions(0),
54 fNTrackMonitorFunctions(0),
55 fNParticleMonitorFunctions(0),
56 fNParticleAndTrackMonitorFunctions(0),
57 fTrackFunctions ( new AliHBTOnePairFctn* [fgkFctnArraySize]),
58 fParticleFunctions ( new AliHBTOnePairFctn* [fgkFctnArraySize]),
59 fParticleAndTrackFunctions ( new AliHBTTwoPairFctn* [fgkFctnArraySize]),
60 fParticleMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),
61 fTrackMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),
62 fParticleAndTrackMonitorFunctions ( new AliHBTMonTwoParticleFctn* [fgkFctnArraySize]),
67 fDisplayMixingInfo(fgkDefaultMixingInfo),
69 fProcessOption(kSimulatedAndReconstructed),
79 /*************************************************************************************/
81 AliHBTAnalysis::AliHBTAnalysis(const AliHBTAnalysis& in):
86 fNParticleFunctions(0),
87 fNParticleAndTrackFunctions(0),
88 fNTrackMonitorFunctions(0),
89 fNParticleMonitorFunctions(0),
90 fNParticleAndTrackMonitorFunctions(0),
92 fParticleFunctions(0x0),
93 fParticleAndTrackFunctions(0x0),
94 fParticleMonitorFunctions(0x0),
95 fTrackMonitorFunctions(0x0),
96 fParticleAndTrackMonitorFunctions(0x0),
100 fBufferSize(fgkDefaultBufferSize),
101 fDisplayMixingInfo(fgkDefaultMixingInfo),
103 fProcessOption(kSimulatedAndReconstructed),
104 fNoCorrfctns(kFALSE),
105 fOutputFileName(0x0),
111 Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
113 /*************************************************************************************/
114 AliHBTAnalysis& AliHBTAnalysis::operator=(const AliHBTAnalysis& /*right*/)
117 Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
120 /*************************************************************************************/
121 AliHBTAnalysis::~AliHBTAnalysis()
124 //note that we do not delete functions itself
125 // they should be deleted by whom where created
126 //we only store pointers, and we use "new" only for pointers array
130 AliDebug(5,"Is Owner: Attempting to delete functions");
132 if (AliVAODParticle::GetDebug()>5)Info("~AliHBTAnalysis","Delete functions done");
134 delete [] fTrackFunctions;
135 delete [] fParticleFunctions;
136 delete [] fParticleAndTrackFunctions;
138 delete [] fParticleMonitorFunctions;
139 delete [] fTrackMonitorFunctions;
140 delete [] fParticleAndTrackMonitorFunctions;
143 delete fOutputFileName;
146 /*************************************************************************************/
148 Int_t AliHBTAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim)
150 //Processes one event
151 if (fProcEvent == 0x0)
153 Error("ProcessEvent","Analysis <<%s>> option not specified.",GetName());
156 if ( Rejected(aodrec,aodsim) )
158 // Double_t x = 0, y = 0, z = 0;
159 // aodrec->GetPrimaryVertex(x,y,z);
160 // Info("ProcessEvent","Event has vertex at %f %f %f",x,y,z);
161 Info("ProcessEvent","Nch is %d",aodsim->GetNumberOfCharged());
162 Info("ProcessEvent","%s: Event cut rejected this event",GetName());
166 //Move event to the apparent vertex -> must be after the event cut
167 //It is important for any cut that use any spacial coordiantes,
168 //f.g. anti-merging cut in order to preserve the same bahavior of variables (f.g. distance between tracks)
169 Double_t dvx = 0, dvy = 0, dvz = 0;
172 Double_t pvx,pvy,pvz;
173 aodrec->GetPrimaryVertex(pvx,pvy,pvz);
175 dvx = fVertexX - pvx;
176 dvy = fVertexY - pvy;
177 dvz = fVertexZ - pvz;
178 aodrec->Move(dvx,dvy,dvz);
181 Int_t result = (this->*fProcEvent)(aodrec,aodsim);
183 if (aodrec) aodrec->Move(-dvx,-dvy,-dvz);//move event back to the same spacial coordinates
187 /*************************************************************************************/
189 Int_t AliHBTAnalysis::Finish()
195 /*************************************************************************************/
197 void AliHBTAnalysis::DeleteFunctions()
199 //Deletes all functions added to analysis
202 for(ii = 0;ii<fNParticleFunctions;ii++)
204 if (AliVAODParticle::GetDebug()>5)
206 Info("DeleteFunctions","Deleting ParticleFunction %#x",fParticleFunctions[ii]);
207 Info("DeleteFunctions","Deleting ParticleFunction %s",fParticleFunctions[ii]->Name());
209 delete fParticleFunctions[ii];
211 fNParticleFunctions = 0;
213 for(ii = 0;ii<fNTrackFunctions;ii++)
215 if (AliVAODParticle::GetDebug()>5)
217 Info("DeleteFunctions","Deleting TrackFunction %#x",fTrackFunctions[ii]);
218 Info("DeleteFunctions","Deleting TrackFunction %s",fTrackFunctions[ii]->Name());
220 delete fTrackFunctions[ii];
222 fNTrackFunctions = 0;
224 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
226 if (AliVAODParticle::GetDebug()>5)
228 Info("DeleteFunctions","Deleting ParticleAndTrackFunction %#x",fParticleAndTrackFunctions[ii]);
229 Info("DeleteFunctions","Deleting ParticleAndTrackFunction %s",fParticleAndTrackFunctions[ii]->Name());
231 delete fParticleAndTrackFunctions[ii];
233 fNParticleAndTrackFunctions = 0;
235 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
237 if (AliVAODParticle::GetDebug()>5)
239 Info("DeleteFunctions","Deleting ParticleMonitorFunction %#x",fParticleMonitorFunctions[ii]);
240 Info("DeleteFunctions","Deleting ParticleMonitorFunction %s",fParticleMonitorFunctions[ii]->Name());
242 delete fParticleMonitorFunctions[ii];
244 fNParticleMonitorFunctions = 0;
246 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
248 if (AliVAODParticle::GetDebug()>5)
250 Info("DeleteFunctions","Deleting TrackMonitorFunction %#x",fTrackMonitorFunctions[ii]);
251 Info("DeleteFunctions","Deleting TrackMonitorFunction %s",fTrackMonitorFunctions[ii]->Name());
253 delete fTrackMonitorFunctions[ii];
255 fNTrackMonitorFunctions = 0;
257 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
259 if (AliVAODParticle::GetDebug()>5)
261 Info("DeleteFunctions","Deleting ParticleAndTrackMonitorFunction %#x",fParticleAndTrackMonitorFunctions[ii]);
262 Info("DeleteFunctions","Deleting ParticleAndTrackMonitorFunction %s",fParticleAndTrackMonitorFunctions[ii]->Name());
264 delete fParticleAndTrackMonitorFunctions[ii];
266 fNParticleAndTrackMonitorFunctions = 0;
269 /*************************************************************************************/
271 Int_t AliHBTAnalysis::Init()
273 //Initializeation method
274 //calls Init for all functions
276 //Depending on option and pair cut assigns proper analysis method
277 Bool_t nonid = IsNonIdentAnalysis();
278 switch(fProcessOption)
281 if (nonid) fProcEvent = &AliHBTAnalysis::ProcessRecNonId;
282 else fProcEvent = &AliHBTAnalysis::ProcessRec;
287 if (nonid) fProcEvent = &AliHBTAnalysis::ProcessSimNonId;
288 else fProcEvent = &AliHBTAnalysis::ProcessSim;
292 case kSimulatedAndReconstructed:
293 if (nonid) fProcEvent = &AliHBTAnalysis::ProcessRecAndSimNonId;
294 else fProcEvent = &AliHBTAnalysis::ProcessRecAndSim;
298 if (fPartBuffer == 0x0) fPartBuffer = new AliEventBuffer (fBufferSize);
299 else fPartBuffer->Reset();
301 if (fTrackBuffer == 0x0) fTrackBuffer = new AliEventBuffer (fBufferSize);
302 else fTrackBuffer->Reset();
305 fNoCorrfctns = (fNParticleFunctions == 0) && (fNTrackFunctions == 0) && (fNParticleAndTrackFunctions == 0);
308 for(ii = 0;ii<fNParticleFunctions;ii++)
309 fParticleFunctions[ii]->Init();
311 for(ii = 0;ii<fNTrackFunctions;ii++)
312 fTrackFunctions[ii]->Init();
314 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
315 fParticleAndTrackFunctions[ii]->Init();
317 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
318 fParticleMonitorFunctions[ii]->Init();
320 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
321 fTrackMonitorFunctions[ii]->Init();
323 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
324 fParticleAndTrackMonitorFunctions[ii]->Init();
328 /*************************************************************************************/
330 void AliHBTAnalysis::ResetFunctions()
332 //In case fOwner is true, deletes all functions
333 //in other case, just set number of analysis to 0
334 if (fIsOwner) DeleteFunctions();
337 fNParticleFunctions = 0;
338 fNTrackFunctions = 0;
339 fNParticleAndTrackFunctions = 0;
340 fNParticleMonitorFunctions = 0;
341 fNTrackMonitorFunctions = 0;
342 fNParticleAndTrackMonitorFunctions = 0;
345 /*************************************************************************************/
346 Int_t AliHBTAnalysis::ProcessRecAndSim(AliAOD* aodrec, AliAOD* aodsim)
348 //Does analysis for both tracks and particles
349 //mainly for resolution study and analysies with weighting algirithms
351 // cut on particles only -- why?
352 // - PID: when we make resolution analysis we want to take only tracks with correct PID
353 // We need cut on tracks because there are data characteristic
355 AliVAODParticle * part1, * part2;
356 AliVAODParticle * track1, * track2;
358 AliAOD * trackEvent = aodrec, *partEvent = aodsim;
359 AliAOD* trackEvent1 = new AliAOD();
360 AliAOD* partEvent1 = new AliAOD();
362 AliAOD * trackEvent2,*partEvent2;
364 // Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
366 // Int_t nev = fReader->GetNumberOfTrackEvents();
367 static AliHBTPair tpair;
368 static AliHBTPair ppair;
370 AliHBTPair* trackpair = &tpair;
371 AliHBTPair* partpair = &ppair;
373 AliHBTPair * tmptrackpair;//temprary pointers to pairs
374 AliHBTPair * tmppartpair;
380 if ( !partEvent || !trackEvent )
382 Error("ProcessRecAndSim","<<%s>> Can not get event",GetName());
386 if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
388 Error("ProcessRecAndSim",
389 "Number of simulated particles (%d) not equal to number of reconstructed tracks (%d). Skipping Event.",
390 partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
395 for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
397 /***************************************/
398 /****** Looping same events ********/
399 /****** filling numerators ********/
400 /***************************************/
401 if ( (j%fDisplayMixingInfo) == 0)
402 Info("ProcessRecAndSim",
403 "Mixing particle %d with particles from the same event",j);
405 part1= partEvent->GetParticle(j);
406 track1= trackEvent->GetParticle(j);
408 Bool_t firstcut = (this->*fkPass1)(part1,track1);
409 if (fBufferSize != 0)
410 if ( (firstcut == kFALSE) || ( (this->*fkPass2)(part1,track1) == kFALSE ) )
412 //accepted by any cut
413 // we have to copy because reader keeps only one event
415 partEvent1->AddParticle(part1);
416 trackEvent1->AddParticle(track1);
419 if (firstcut) continue;
421 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
422 fParticleMonitorFunctions[ii]->Process(part1);
423 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
424 fTrackMonitorFunctions[ii]->Process(track1);
425 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
426 fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
428 if (fNoCorrfctns) continue;
430 for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
432 part2= partEvent->GetParticle(k);
433 if (part1->GetUID() == part2->GetUID()) continue;
434 partpair->SetParticles(part1,part2);
436 track2= trackEvent->GetParticle(k);
437 trackpair->SetParticles(track1,track2);
439 if( (this->*fkPass)(partpair,trackpair) ) //check pair cut
440 { //do not meets crietria of the pair cut, try with swapped pairs
441 if( (this->*fkPass)((AliHBTPair*)partpair->GetSwappedPair(),(AliHBTPair*)trackpair->GetSwappedPair()) )
442 continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
444 { //swaped pair meets all the criteria
445 tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
446 tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
450 {//meets criteria of the pair cut
451 tmptrackpair = trackpair;
452 tmppartpair = partpair;
455 for(ii = 0;ii<fNParticleFunctions;ii++)
456 fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
458 for(ii = 0;ii<fNTrackFunctions;ii++)
459 fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
461 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
462 fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair,tmppartpair);
463 //end of 2nd loop over particles from the same event
464 }//for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
466 /***************************************/
467 /***** Filling denominators *********/
468 /***************************************/
469 if (fBufferSize == 0) continue;
471 fPartBuffer->ResetIter();
472 fTrackBuffer->ResetIter();
474 while (( partEvent2 = fPartBuffer->Next() ))
476 trackEvent2 = fTrackBuffer->Next();
479 if ( (j%fDisplayMixingInfo) == 0)
480 Info("ProcessRecAndSim",
481 "Mixing particle %d from current event with particles from event %d",j,-m);
483 for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
485 part2= partEvent2->GetParticle(l);
486 partpair->SetParticles(part1,part2);
488 track2= trackEvent2->GetParticle(l);
489 trackpair->SetParticles(track1,track2);
491 if( (this->*fkPass)(partpair,trackpair) ) //check pair cut
492 { //do not meets crietria of the
493 if( (this->*fkPass)((AliHBTPair*)partpair->GetSwappedPair(),(AliHBTPair*)trackpair->GetSwappedPair()) )
497 tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
498 tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
502 {//meets criteria of the pair cut
503 tmptrackpair = trackpair;
504 tmppartpair = partpair;
507 for(ii = 0;ii<fNParticleFunctions;ii++)
508 fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
510 for(ii = 0;ii<fNTrackFunctions;ii++)
511 fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
513 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
514 fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair);
515 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
518 //end of loop over particles from first event
519 }//for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
520 delete fPartBuffer->Push(partEvent1);
521 delete fTrackBuffer->Push(trackEvent1);
522 //end of loop over events
525 /*************************************************************************************/
527 Int_t AliHBTAnalysis::ProcessSim(AliAOD* /*aodrec*/, AliAOD* aodsim)
529 //Does analysis of simulated data
530 AliVAODParticle * part1, * part2;
532 AliAOD* partEvent = aodsim;
533 AliAOD* partEvent1 = new AliAOD();
539 AliHBTPair* partpair = &ppair;
541 AliHBTPair * tmppartpair;
547 Error("ProcessSim","Can not get event");
552 for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
554 /***************************************/
555 /****** Looping same events ********/
556 /****** filling numerators ********/
557 /***************************************/
558 if ( (j%fDisplayMixingInfo) == 0)
560 "Mixing particle %d with particles from the same event",j);
562 part1= partEvent->GetParticle(j);
564 Bool_t firstcut = fPairCut->GetFirstPartCut()->Rejected(part1);
566 if (fBufferSize != 0)
567 if ( (firstcut == kFALSE) || ( fPairCut->GetSecondPartCut()->Rejected(part1) == kFALSE ) )
569 //accepted by any cut
570 // we have to copy because reader keeps only one event
572 partEvent1->AddParticle(part1);
575 if (firstcut) continue;
577 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
578 fParticleMonitorFunctions[ii]->Process(part1);
580 if ( fNParticleFunctions == 0 ) continue;
582 for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
584 part2= partEvent->GetParticle(k);
585 if (part1->GetUID() == part2->GetUID()) continue;
586 partpair->SetParticles(part1,part2);
588 if(fPairCut->Rejected(partpair)) //check pair cut
589 { //do not meets crietria of the
590 if( fPairCut->Rejected((AliHBTPair*)partpair->GetSwappedPair()) ) continue;
591 else tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
595 tmppartpair = partpair;
598 for(ii = 0;ii<fNParticleFunctions;ii++)
599 fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
601 //end of 2nd loop over particles from the same event
602 }//for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
604 /***************************************/
605 /***** Filling denominators *********/
606 /***************************************/
607 if (fBufferSize == 0) continue;
609 fPartBuffer->ResetIter();
611 while (( partEvent2 = fPartBuffer->Next() ))
614 if ( (j%fDisplayMixingInfo) == 0)
616 "Mixing particle %d from current event with particles from event %d",j,-m);
617 for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
620 part2= partEvent2->GetParticle(l);
621 partpair->SetParticles(part1,part2);
623 if( fPairCut->Rejected(partpair) ) //check pair cut
624 { //do not meets crietria of the
625 if( fPairCut->Rejected((AliHBTPair*)partpair->GetSwappedPair()) )
629 tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
633 {//meets criteria of the pair cut
634 tmppartpair = partpair;
637 for(ii = 0;ii<fNParticleFunctions;ii++)
638 fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
640 }//for(Int_t l = 0; l<N2;l++) // ... on all particles
645 delete fPartBuffer->Push(partEvent1);
646 //end of loop over events
649 /*************************************************************************************/
650 Int_t AliHBTAnalysis::ProcessRec(AliAOD* aodrec, AliAOD* /*aodsim*/)
652 //Does analysis of reconstructed data
653 AliVAODParticle * track1, * track2;
655 AliAOD* trackEvent = aodrec;
656 AliAOD* trackEvent1 = new AliAOD();
662 AliHBTPair* trackpair = &tpair;
664 AliHBTPair * tmptrackpair;
670 Error("ProcessRecAndSim","Can not get event");
675 for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
677 /***************************************/
678 /****** Looping same events ********/
679 /****** filling numerators ********/
680 /***************************************/
681 if ( (j%fDisplayMixingInfo) == 0)
683 "Mixing Particle %d with Particles from the same event",j);
685 track1= trackEvent->GetParticle(j);
687 Bool_t firstcut = fPairCut->GetFirstPartCut()->Rejected(track1);
689 if (fBufferSize != 0)
690 if ( (firstcut == kFALSE) || ( fPairCut->GetSecondPartCut()->Rejected(track1) == kFALSE ) )
692 //accepted by any cut
693 // we have to copy because reader keeps only one event
695 trackEvent1->AddParticle(track1);
698 if (firstcut) continue;
700 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
701 fParticleMonitorFunctions[ii]->Process(track1);
703 if ( fNTrackFunctions == 0 ) continue;
705 for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
707 track2= trackEvent->GetParticle(k);
708 if (track1->GetUID() == track2->GetUID()) continue;
709 trackpair->SetParticles(track1,track2);
711 if(fPairCut->Rejected(trackpair)) //check pair cut
712 { //do not meets crietria of the
713 if( fPairCut->Rejected((AliHBTPair*)trackpair->GetSwappedPair()) ) continue;
714 else tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
718 tmptrackpair = trackpair;
721 for(ii = 0;ii<fNTrackFunctions;ii++)
722 fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
724 //end of 2nd loop over Particles from the same event
725 }//for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
727 /***************************************/
728 /***** Filling denominators *********/
729 /***************************************/
730 if (fBufferSize == 0) continue;
732 fTrackBuffer->ResetIter();
734 while (( trackEvent2 = fTrackBuffer->Next() ))
737 if ( (j%fDisplayMixingInfo) == 0)
739 "Mixing Particle %d from current event with Particles from event %d",j,-m);
740 for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++) // ... on all Particles
743 track2= trackEvent2->GetParticle(l);
744 trackpair->SetParticles(track1,track2);
746 if( fPairCut->Rejected(trackpair) ) //check pair cut
747 { //do not meets crietria of the
748 if( fPairCut->Rejected((AliHBTPair*)trackpair->GetSwappedPair()) )
752 tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
756 {//meets criteria of the pair cut
757 tmptrackpair = trackpair;
760 for(ii = 0;ii<fNTrackFunctions;ii++)
761 fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
763 }//for(Int_t l = 0; l<N2;l++) // ... on all Particles
766 delete fTrackBuffer->Push(trackEvent1);
767 //end of loop over events
770 /*************************************************************************************/
772 Int_t AliHBTAnalysis::ProcessRecAndSimNonId(AliAOD* aodrec, AliAOD* aodsim)
774 //Analyzes both reconstructed and simulated data
777 Error("ProcessTracksAndParticlesNonIdentAnal","Reconstructed event is NULL");
783 Error("ProcessTracksAndParticlesNonIdentAnal","Simulated event is NULL");
787 if ( aodrec->GetNumberOfParticles() != aodsim->GetNumberOfParticles() )
789 Error("ProcessTracksAndParticlesNonIdentAnal",
790 "Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
791 aodsim->GetNumberOfParticles() , aodrec->GetNumberOfParticles());
796 AliVAODParticle * part1, * part2;
797 AliVAODParticle * track1, * track2;
799 static AliAOD aodrec1;
800 static AliAOD aodsim1;
802 AliAOD * trackEvent1=&aodrec1,*partEvent1=&aodsim1;//Particle that passes first particle cut, this event
803 trackEvent1->Reset();
805 AliAOD * trackEvent2=0x0,*partEvent2=0x0;//Particle that passes second particle cut, this event
806 AliAOD * trackEvent3=0x0,*partEvent3=0x0;//Particle that passes second particle cut, events from buffer
808 AliAOD* rawtrackEvent = aodrec;//this we get from Reader
809 AliAOD* rawpartEvent = aodsim;//this we get from Reader
811 static AliHBTPair tpair;
812 static AliHBTPair ppair;
814 AliHBTPair* trackpair = &tpair;
815 AliHBTPair* partpair = &ppair;
819 /********************************/
821 /********************************/
822 if ( ( (partEvent2==0x0) || (trackEvent2==0x0)) )//in case fBufferSize == 0 and pointers are created do not eneter
824 partEvent2 = new AliAOD();
825 trackEvent2 = new AliAOD();
828 FilterOut(partEvent1, partEvent2, rawpartEvent, trackEvent1, trackEvent2, rawtrackEvent);
830 for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
832 if ( (j%fDisplayMixingInfo) == 0)
833 Info("ProcessTracksAndParticlesNonIdentAnal",
834 "Mixing particle %d from current event with particles from current event",j);
836 part1= partEvent1->GetParticle(j);
837 track1= trackEvent1->GetParticle(j);
840 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
841 fParticleMonitorFunctions[ii]->Process(part1);
842 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
843 fTrackMonitorFunctions[ii]->Process(track1);
844 for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
845 fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
847 if (fNoCorrfctns) continue;
849 /***************************************/
850 /****** filling numerators ********/
851 /****** (particles from event2) ********/
852 /***************************************/
854 for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups
856 part2= partEvent2->GetParticle(k);
857 if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
858 partpair->SetParticles(part1,part2);
860 track2= trackEvent2->GetParticle(k);
861 trackpair->SetParticles(track1,track2);
863 if( (this->*fkPassPairProp)(partpair,trackpair) ) //check pair cut
864 { //do not meets crietria of the pair cut
868 {//meets criteria of the pair cut
869 for(ii = 0;ii<fNParticleFunctions;ii++)
870 fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
872 for(ii = 0;ii<fNTrackFunctions;ii++)
873 fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
875 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
876 fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(trackpair,partpair);
880 if ( fBufferSize == 0) continue;//do not mix diff histograms
881 /***************************************/
882 /***** Filling denominators *********/
883 /***************************************/
884 fPartBuffer->ResetIter();
885 fTrackBuffer->ResetIter();
889 while ( (partEvent3 = fPartBuffer->Next() ) != 0x0)
891 trackEvent3 = fTrackBuffer->Next();
893 if ( (j%fDisplayMixingInfo) == 0)
894 Info("ProcessTracksAndParticlesNonIdentAnal",
895 "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
897 for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
899 part2= partEvent3->GetParticle(k);
900 partpair->SetParticles(part1,part2);
902 track2= trackEvent3->GetParticle(k);
903 trackpair->SetParticles(track1,track2);
905 if( (this->*fkPassPairProp)(partpair,trackpair) ) //check pair cut
906 { //do not meets crietria of the pair cut
910 {//meets criteria of the pair cut
912 for(iPF = 0;iPF<fNParticleFunctions;iPF++)
913 fParticleFunctions[iPF]->ProcessDiffEventParticles(partpair);
915 for(iPF = 0;iPF<fNTrackFunctions;iPF++)
916 fTrackFunctions[iPF]->ProcessDiffEventParticles(trackpair);
918 for(iPF = 0;iPF<fNParticleAndTrackFunctions;iPF++)
919 fParticleAndTrackFunctions[iPF]->ProcessDiffEventParticles(trackpair,partpair);
921 }// for particles event2
923 }//for over particles in event1
925 delete fPartBuffer->Push(partEvent2);
926 delete fTrackBuffer->Push(trackEvent2);
930 /*************************************************************************************/
931 Int_t AliHBTAnalysis::ProcessSimNonId(AliAOD* /*aodrec*/, AliAOD* aodsim)
933 //does analysis of simulated (MC) data in non-identical mode
934 //i.e. when particles selected by first part. cut are a disjunctive set than particles
935 //passed by the second part. cut
942 AliVAODParticle * part1, * part2;
944 static AliAOD aodsim1;
946 AliAOD* partEvent1=&aodsim1;//Particle that passes first particle cut, this event
948 AliAOD* partEvent2=0x0;//Particle that passes second particle cut, this event
949 AliAOD* partEvent3=0x0;//Particle that passes second particle cut, events from buffer
951 AliAOD* rawpartEvent = aodsim;//this we get from Reader
953 static AliHBTPair ppair;
955 AliHBTPair* partpair = &ppair;
959 /********************************/
961 /********************************/
962 if (partEvent2==0x0)//in case fBufferSize == 0 and pointers are created do not eneter
964 partEvent2 = new AliAOD();
967 FilterOut(partEvent1, partEvent2, rawpartEvent);
969 for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
971 if ( (j%fDisplayMixingInfo) == 0)
972 Info("ProcessParticlesNonIdentAnal",
973 "Mixing particle %d from current event with particles from current event",j);
975 part1= partEvent1->GetParticle(j);
978 for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
979 fParticleMonitorFunctions[ii]->Process(part1);
981 if (fNParticleFunctions == 0) continue;
983 /***************************************/
984 /****** filling numerators ********/
985 /****** (particles from event2) ********/
986 /***************************************/
988 for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups
990 part2= partEvent2->GetParticle(k);
991 if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
992 partpair->SetParticles(part1,part2);
995 if(fPairCut->PassPairProp(partpair) ) //check pair cut
996 { //do not meets crietria of the pair cut
1000 {//meets criteria of the pair cut
1001 for(ii = 0;ii<fNParticleFunctions;ii++)
1002 fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
1006 if ( fBufferSize == 0) continue;//do not mix diff histograms
1007 /***************************************/
1008 /***** Filling denominators *********/
1009 /***************************************/
1010 fPartBuffer->ResetIter();
1014 while ( (partEvent3 = fPartBuffer->Next() ) != 0x0)
1017 if ( (j%fDisplayMixingInfo) == 0)
1018 Info("ProcessParticlesNonIdentAnal",
1019 "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
1021 for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
1023 part2= partEvent3->GetParticle(k);
1024 partpair->SetParticles(part1,part2);
1027 if(fPairCut->PassPairProp(partpair) ) //check pair cut
1028 { //do not meets crietria of the pair cut
1032 {//meets criteria of the pair cut
1033 for(ii = 0;ii<fNParticleFunctions;ii++)
1035 fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
1038 }// for particles event2
1040 }//for over particles in event1
1042 delete fPartBuffer->Push(partEvent2);
1046 /*************************************************************************************/
1047 Int_t AliHBTAnalysis::ProcessRecNonId(AliAOD* aodrec, AliAOD* /*aodsim*/)
1049 //Analyzes both reconstructed and simulated data
1055 AliVAODParticle * track1, * track2;
1057 static AliAOD aodrec1;
1058 AliAOD * trackEvent1=&aodrec1;//Particle that passes first particle cut, this event
1059 trackEvent1->Reset();
1060 AliAOD * trackEvent2=0x0;//Particle that passes second particle cut, this event
1061 AliAOD * trackEvent3=0x0;//Particle that passes second particle cut, events from buffer
1062 AliAOD* rawtrackEvent = aodrec;//this we get from Reader
1064 static AliHBTPair tpair;
1066 AliHBTPair* trackpair = &tpair;
1071 /********************************/
1073 /********************************/
1074 if ( trackEvent2==0x0 )//in case fBufferSize == 0 and pointers are created do not eneter
1076 trackEvent2 = new AliAOD();
1079 FilterOut(trackEvent1, trackEvent2, rawtrackEvent);
1081 for (Int_t j = 0; j<trackEvent1->GetNumberOfParticles() ; j++)
1083 if ( (j%fDisplayMixingInfo) == 0)
1084 Info("ProcessTracksNonIdentAnal",
1085 "Mixing particle %d from current event with particles from current event",j);
1087 track1= trackEvent1->GetParticle(j);
1090 for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
1091 fTrackMonitorFunctions[ii]->Process(track1);
1093 if (fNTrackFunctions == 0x0) continue;
1095 /***************************************/
1096 /****** filling numerators ********/
1097 /****** (particles from event2) ********/
1098 /***************************************/
1100 for (Int_t k = 0; k < trackEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups
1102 track2= trackEvent2->GetParticle(k);
1103 if (track1->GetUID() == track2->GetUID()) continue;//this is the same particle but with different PID
1104 trackpair->SetParticles(track1,track2);
1107 if( fPairCut->PassPairProp(trackpair)) //check pair cut
1108 { //do not meets crietria of the pair cut
1112 {//meets criteria of the pair cut
1114 for(iTF = 0;iTF<fNTrackFunctions;iTF++)
1115 fTrackFunctions[iTF]->ProcessSameEventParticles(trackpair);
1119 if ( fBufferSize == 0) continue;//do not mix diff histograms
1120 /***************************************/
1121 /***** Filling denominators *********/
1122 /***************************************/
1123 fTrackBuffer->ResetIter();
1127 while ( (trackEvent3 = fTrackBuffer->Next() ) != 0x0)
1129 if ( (j%fDisplayMixingInfo) == 0)
1130 Info("ProcessTracksNonIdentAnal",
1131 "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
1133 for (Int_t k = 0; k < trackEvent3->GetNumberOfParticles() ; k++)
1135 track2= trackEvent3->GetParticle(k);
1136 trackpair->SetParticles(track1,track2);
1138 if( fPairCut->PassPairProp(trackpair)) //check pair cut
1139 { //do not meets crietria of the pair cut
1143 {//meets criteria of the pair cut
1144 for(ii = 0;ii<fNTrackFunctions;ii++)
1145 fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
1147 }// for particles event2
1149 }//for over particles in event1
1151 delete fTrackBuffer->Push(trackEvent2);
1155 /*************************************************************************************/
1157 void AliHBTAnalysis::Process(Option_t* option)
1159 //default option = "TracksAndParticles"
1160 //Main method of the HBT Analysis Package
1161 //It triggers reading with the global cut (default is an empty cut)
1162 //Than it checks options and data which are read
1163 //if everything is OK, then it calls one of the looping methods
1164 //depending on tfReaderhe option
1165 //These methods differs on what thay are looping on
1168 //--------------------------------------------------------------------
1169 //ProcessTracksAndParticles - "TracksAndParticles"
1171 // it loops over both, tracks(reconstructed) and particles(simulated)
1172 // all function gethered in all 3 lists are called for each (double)pair
1174 //ProcessTracks - "Tracks"
1175 // it loops only on tracks(reconstructed),
1176 // functions ONLY from fTrackFunctions list are called
1178 //ProcessParticles - "Particles"
1179 // it loops only on particles(simulated),
1180 // functions ONLY from fParticleAndTrackFunctions list are called
1185 Error("Process","The reader is not set");
1189 const char *oT = strstr(option,"Tracks");
1190 const char *oP = strstr(option,"Particles");
1192 Bool_t nonid = IsNonIdentAnalysis();
1198 if (nonid) ProcessTracksAndParticlesNonIdentAnal();
1199 else ProcessTracksAndParticles();
1205 if (nonid) ProcessTracksNonIdentAnal();
1206 else ProcessTracks();
1212 if (nonid) ProcessParticlesNonIdentAnal();
1213 else ProcessParticles();
1218 /*************************************************************************************/
1220 void AliHBTAnalysis::ProcessTracksAndParticles()
1222 //Makes analysis for both tracks and particles
1223 //mainly for resolution study and analysies with weighting algirithms
1224 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1225 //the loops are splited
1227 // cut on particles only -- why?
1228 // - PID: when we make resolution analysis we want to take only tracks with correct PID
1229 // We need cut on tracks because there are data characteristic to
1231 AliAOD * trackEvent, *partEvent;
1234 while (fReader->Next() == kFALSE)
1236 partEvent = fReader->GetEventSim();
1237 trackEvent = fReader->GetEventRec();
1238 ProcessRecAndSim(trackEvent,partEvent);
1240 }//while (fReader->Next() == kFALSE)
1243 /*************************************************************************************/
1245 void AliHBTAnalysis::ProcessTracks()
1247 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1248 //the loops are splited
1249 AliAOD * trackEvent;
1251 while (fReader->Next() == kFALSE)
1253 trackEvent = fReader->GetEventRec();
1254 ProcessRec(trackEvent,0x0);
1255 }//while (fReader->Next() == kFALSE)
1258 /*************************************************************************************/
1260 void AliHBTAnalysis::ProcessParticles()
1262 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
1263 //the loops are splited
1266 while (fReader->Next() == kFALSE)
1268 partEvent = fReader->GetEventSim();
1269 ProcessSim(0x0,partEvent);
1270 }//while (fReader->Next() == kFALSE)
1272 /*************************************************************************************/
1274 void AliHBTAnalysis::WriteFunctions()
1276 //Calls Write for all defined functions in analysis
1277 //== writes all results
1278 TFile* oututfile = 0x0;
1279 if (fOutputFileName)
1281 oututfile = TFile::Open(*fOutputFileName,"update");
1284 for(ii = 0;ii<fNParticleFunctions;ii++)
1286 if (AliVAODParticle::GetDebug()>5)
1288 Info("WriteFunctions","Writing ParticleFunction %#x",fParticleFunctions[ii]);
1289 Info("WriteFunctions","Writing ParticleFunction %s",fParticleFunctions[ii]->Name());
1291 fParticleFunctions[ii]->Write();
1294 for(ii = 0;ii<fNTrackFunctions;ii++)
1296 if (AliVAODParticle::GetDebug()>5)
1298 Info("WriteFunctions","Writing TrackFunction %#x",fTrackFunctions[ii]);
1299 Info("WriteFunctions","Writing TrackFunction %s",fTrackFunctions[ii]->Name());
1301 fTrackFunctions[ii]->Write();
1304 for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
1306 if (AliVAODParticle::GetDebug()>5)
1308 Info("WriteFunctions","Writing ParticleAndTrackFunction %#x",fParticleAndTrackFunctions[ii]);
1309 Info("WriteFunctions","Writing ParticleAndTrackFunction %s",fParticleAndTrackFunctions[ii]->Name());
1311 fParticleAndTrackFunctions[ii]->Write();
1314 for(ii = 0;ii<fNParticleMonitorFunctions;ii++)
1316 if (AliVAODParticle::GetDebug()>5)
1318 Info("WriteFunctions","Writing ParticleMonitorFunction %#x",fParticleMonitorFunctions[ii]);
1319 Info("WriteFunctions","Writing ParticleMonitorFunction %s",fParticleMonitorFunctions[ii]->Name());
1321 fParticleMonitorFunctions[ii]->Write();
1324 for(ii = 0;ii<fNTrackMonitorFunctions;ii++)
1326 if (AliVAODParticle::GetDebug()>5)
1328 Info("WriteFunctions","Writing TrackMonitorFunction %#x",fTrackMonitorFunctions[ii]);
1329 Info("WriteFunctions","Writing TrackMonitorFunction %s",fTrackMonitorFunctions[ii]->Name());
1331 fTrackMonitorFunctions[ii]->Write();
1334 for(ii = 0;ii<fNParticleAndTrackMonitorFunctions;ii++)
1336 if (AliVAODParticle::GetDebug()>5)
1338 Info("WriteFunctions","Writing ParticleAndTrackMonitorFunction %#x",fParticleAndTrackMonitorFunctions[ii]);
1339 Info("WriteFunctions","Writing ParticleAndTrackMonitorFunction %s",fParticleAndTrackMonitorFunctions[ii]->Name());
1341 fParticleAndTrackMonitorFunctions[ii]->Write();
1345 /*************************************************************************************/
1347 void AliHBTAnalysis::SetOutputFileName(const char* fname)
1349 //Sets fiele name where to dump results,
1350 //if not specified reults are written to gDirectory
1353 delete fOutputFileName;
1354 fOutputFileName = 0x0;
1357 if ( strcmp(fname,"") == 0 )
1359 delete fOutputFileName;
1360 fOutputFileName = 0x0;
1363 if (fOutputFileName == 0x0) fOutputFileName = new TString(fname);
1364 else *fOutputFileName = fname;
1366 /*************************************************************************************/
1368 void AliHBTAnalysis::SetGlobalPairCut(AliAODPairCut* cut)
1370 //Sets the global cut
1373 Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
1376 fPairCut = (AliAODPairCut*)cut->Clone();
1379 /*************************************************************************************/
1381 void AliHBTAnalysis::AddTrackFunction(AliHBTOnePairFctn* f)
1383 //Adds track function
1384 if (f == 0x0) return;
1385 if (fNTrackFunctions == fgkFctnArraySize)
1387 Error("AliHBTAnalysis::AddTrackFunction","Can not add this function, not enough place in the array.");
1389 fTrackFunctions[fNTrackFunctions] = f;
1392 /*************************************************************************************/
1394 void AliHBTAnalysis::AddParticleFunction(AliHBTOnePairFctn* f)
1396 //adds particle function
1397 if (f == 0x0) return;
1399 if (fNParticleFunctions == fgkFctnArraySize)
1401 Error("AliHBTAnalysis::AddParticleFunction","Can not add this function, not enough place in the array.");
1403 fParticleFunctions[fNParticleFunctions] = f;
1404 fNParticleFunctions++;
1406 /*************************************************************************************/
1408 void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTTwoPairFctn* f)
1410 //add resolution function
1411 if (f == 0x0) return;
1412 if (fNParticleAndTrackFunctions == fgkFctnArraySize)
1414 Error("AliHBTAnalysis::AddParticleAndTrackFunction","Can not add this function, not enough place in the array.");
1416 fParticleAndTrackFunctions[fNParticleAndTrackFunctions] = f;
1417 fNParticleAndTrackFunctions++;
1419 /*************************************************************************************/
1421 void AliHBTAnalysis::AddParticleMonitorFunction(AliHBTMonOneParticleFctn* f)
1423 //add particle monitoring function
1424 if (f == 0x0) return;
1426 if (fNParticleMonitorFunctions == fgkFctnArraySize)
1428 Error("AliHBTAnalysis::AddParticleMonitorFunction","Can not add this function, not enough place in the array.");
1430 fParticleMonitorFunctions[fNParticleMonitorFunctions] = f;
1431 fNParticleMonitorFunctions++;
1433 /*************************************************************************************/
1435 void AliHBTAnalysis::AddTrackMonitorFunction(AliHBTMonOneParticleFctn* f)
1437 //add track monitoring function
1438 if (f == 0x0) return;
1440 if (fNTrackMonitorFunctions == fgkFctnArraySize)
1442 Error("AliHBTAnalysis::AddTrackMonitorFunction","Can not add this function, not enough place in the array.");
1444 fTrackMonitorFunctions[fNTrackMonitorFunctions] = f;
1445 fNTrackMonitorFunctions++;
1447 /*************************************************************************************/
1449 void AliHBTAnalysis::AddParticleAndTrackMonitorFunction(AliHBTMonTwoParticleFctn* f)
1451 //add resolution monitoring function
1452 if (f == 0x0) return;
1453 if (fNParticleAndTrackMonitorFunctions == fgkFctnArraySize)
1455 Error("AliHBTAnalysis::AddParticleAndTrackMonitorFunction","Can not add this function, not enough place in the array.");
1457 fParticleAndTrackMonitorFunctions[fNParticleAndTrackMonitorFunctions] = f;
1458 fNParticleAndTrackMonitorFunctions++;
1462 /*************************************************************************************/
1463 /*************************************************************************************/
1465 Bool_t AliHBTAnalysis::RunCoherencyCheck()
1467 //Checks if both HBTRuns are similar
1468 //return true if error found
1469 //if they seem to be OK return false
1473 Info("RunCoherencyCheck","Checking HBT Runs Coherency");
1475 //When we use non-buffering reader this is a big waste of time -> We need to read all data to check it
1476 //and reader is implemented safe in this case anyway
1477 // Info("RunCoherencyCheck","Number of events ...");
1478 // if (fReader->GetNumberOfPartEvents() == fReader->GetNumberOfTrackEvents() ) //check whether there is the same number of events
1480 // Info("RunCoherencyCheck","OK. %d found\n",fReader->GetNumberOfTrackEvents());
1483 // { //if not the same - ERROR
1484 // Error("RunCoherencyCheck",
1485 // "Number of simulated events (%d) is not equal to number of reconstructed events(%d)",
1486 // fReader->GetNumberOfPartEvents(),fReader->GetNumberOfTrackEvents());
1490 Info("RunCoherencyCheck","Checking number of Particles AND Particles Types in each event ...");
1494 for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
1496 partEvent= fReader->GetEventSim(i); //gets the "ith" event
1497 trackEvent = fReader->GetEventRec(i);
1499 if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
1500 if ( (partEvent == 0x0) || (partEvent == 0x0) )
1502 Error("RunCoherencyCheck",
1503 "One event is NULL and the other one not. Event Number %d",i);
1507 if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
1509 Error("RunCoherencyCheck",
1510 "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
1511 i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
1515 for (Int_t j = 0; j<partEvent->GetNumberOfParticles(); j++)
1517 if( partEvent->GetParticle(j)->GetPdgCode() != trackEvent->GetParticle(j)->GetPdgCode() )
1519 Error("RunCoherencyCheck",
1520 "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
1521 i,j, partEvent->GetParticle(j)->GetPdgCode(),trackEvent->GetParticle(j)->GetPdgCode() );
1527 Info("RunCoherencyCheck"," Done");
1528 Info("RunCoherencyCheck"," Everything looks OK");
1533 /*************************************************************************************/
1535 void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
1537 //Performs analysis for both, tracks and particles
1538 AliAOD* rawtrackEvent, * rawpartEvent;
1541 Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
1542 Info("ProcessTracksAndParticlesNonIdentAnal","***** NON IDENT MODE ****************");
1543 Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
1545 for (Int_t i = 0;;i++)//infinite loop
1547 if (fReader->Next()) break; //end when no more events available
1549 rawpartEvent = fReader->GetEventSim();
1550 rawtrackEvent = fReader->GetEventRec();
1552 ProcessRecAndSimNonId(rawtrackEvent,rawpartEvent);
1553 }//end of loop over events (1)
1555 /*************************************************************************************/
1557 void AliHBTAnalysis::ProcessTracksNonIdentAnal()
1559 //Process Tracks only with non identical mode
1560 AliAOD * rawtrackEvent;
1563 Info("ProcessTracksNonIdentAnal","**************************************");
1564 Info("ProcessTracksNonIdentAnal","***** NON IDENT MODE ****************");
1565 Info("ProcessTracksNonIdentAnal","**************************************");
1567 for (Int_t i = 0;;i++)//infinite loop
1569 if (fReader->Next()) break; //end when no more events available
1570 rawtrackEvent = fReader->GetEventRec();
1571 ProcessRecNonId(rawtrackEvent,0x0);
1572 }//end of loop over events (1)
1574 /*************************************************************************************/
1576 void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
1578 //process paricles only with non identical mode
1579 AliAOD * rawpartEvent = 0x0;
1582 Info("ProcessParticlesNonIdentAnal","**************************************");
1583 Info("ProcessParticlesNonIdentAnal","***** NON IDENT MODE ****************");
1584 Info("ProcessParticlesNonIdentAnal","**************************************");
1586 for (Int_t i = 0;;i++)//infinite loop
1588 if (fReader->Next()) break; //end when no more events available
1590 rawpartEvent = fReader->GetEventSim();
1591 ProcessSimNonId(0x0,rawpartEvent);
1592 }//end of loop over events (1)
1595 /*************************************************************************************/
1596 void AliHBTAnalysis::FilterOut(AliAOD* outpart1, AliAOD* outpart2, AliAOD* inpart,
1597 AliAOD* outtrack1, AliAOD* outtrack2, AliAOD* intrack) const
1599 //Puts particles accepted as a first particle by global cut in out1
1600 //and as a second particle in out2
1602 AliVAODParticle* part, *track;
1611 for (Int_t i = 0; i < inpart->GetNumberOfParticles(); i++)
1614 part = inpart->GetParticle(i);
1615 track = intrack->GetParticle(i);
1617 if ( ((this->*fkPass1)(part,track)) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
1618 if ( ((this->*fkPass2)(part,track)) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
1620 if (gDebug)//to be removed in real analysis
1621 if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1623 //Particle accpted by both cuts
1624 Error("FilterOut","Particle accepted by both cuts");
1630 outpart1->AddParticle(part);
1631 outtrack1->AddParticle(track);
1637 outpart2->AddParticle(part);
1638 outtrack2->AddParticle(track);
1643 /*************************************************************************************/
1645 void AliHBTAnalysis::FilterOut(AliAOD* out1, AliAOD* out2, AliAOD* in) const
1647 //Puts particles accepted as a first particle by global cut in out1
1648 //and as a second particle in out2
1649 AliVAODParticle* part;
1654 AliAODParticleCut *cut1 = fPairCut->GetFirstPartCut();
1655 AliAODParticleCut *cut2 = fPairCut->GetSecondPartCut();
1659 for (Int_t i = 0; i < in->GetNumberOfParticles(); i++)
1662 part = in->GetParticle(i);
1664 if ( cut1->Rejected(part) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
1665 if ( cut2->Rejected(part) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
1667 if (gDebug)//to be removed in real analysis
1668 if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
1670 //Particle accpted by both cuts
1671 Error("FilterOut","Particle accepted by both cuts");
1677 out1->AddParticle(part);
1683 out2->AddParticle(part);
1688 /*************************************************************************************/
1690 Bool_t AliHBTAnalysis::IsNonIdentAnalysis()
1692 //checks if it is possible to use special analysis for non identical particles
1693 //it means - in global pair cut first particle id is different than second one
1694 //and both are different from 0
1695 //in the future is possible to perform more sophisticated check
1696 //if cuts have excluding requirements
1698 if (fPairCut->IsEmpty())
1701 if (fPairCut->GetFirstPartCut()->IsEmpty())
1704 if (fPairCut->GetSecondPartCut()->IsEmpty())
1707 Int_t id1 = fPairCut->GetFirstPartCut()->GetPID();
1708 Int_t id2 = fPairCut->GetSecondPartCut()->GetPID();
1710 if ( (id1==0) || (id2==0) || (id1==id2) )
1715 /*************************************************************************************/
1717 void AliHBTAnalysis::SetApparentVertex(Double_t x, Double_t y, Double_t z)
1719 //Sets apparent vertex
1720 // All events have to be moved to the same vertex position in order to
1721 // to be able to comare any space positions (f.g. anti-merging)
1722 // This method defines this position
1728 /*************************************************************************************/
1730 void AliHBTAnalysis::PressAnyKey()
1732 //small utility function that helps to make comfortable macros
1735 fcntl(0, F_SETFL, O_NONBLOCK);
1736 ::Info("","Press Any Key to continue ...");
1739 nread = read(0, &c, 1);
1740 gSystem->ProcessEvents();