-
#include "AliHBTAnalysis.h"
-
-#include <iostream.h>
-
-#include "AliHBTRun.h"
-#include "AliHBTEvent.h"
-#include "AliHBTReader.h"
-#include "AliHBTParticle.h"
-#include "AliHBTParticleCut.h"
+//_________________________________________________________
+////////////////////////////////////////////////////////////////////////////
+//
+// class AliHBTAnalysis
+//
+// Central Object Of HBTAnalyser:
+// This class performs main looping within HBT Analysis
+// User must plug a reader of Type AliReader
+// User plugs in coorelation and monitor functions
+// as well as monitor functions
+//
+// HBT Analysis Tool, which is integral part of AliRoot,
+// ALICE Off-Line framework:
+//
+// Piotr.Skowronski@cern.ch
+// more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
+//
+////////////////////////////////////////////////////////////////////////////
+//_________________________________________________________
+
+#include <Riostream.h>
+#include <AliESD.h>
+
+#include <TSystem.h>
+#include <TFile.h>
+
+#include "AliAOD.h"
+#include "AliAODParticle.h"
+#include "AliAODPairCut.h"
+#include "AliEventCut.h"
+
+#include "AliEventBuffer.h"
+
+#include "AliReader.h"
#include "AliHBTPair.h"
-#include "AliHBTPairCut.h"
#include "AliHBTFunction.h"
-
-#include <TBenchmark.h>
-#include <TList.h>
-
-
+#include "AliHBTMonitorFunction.h"
+
ClassImp(AliHBTAnalysis)
const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100;
-const Int_t AliHBTAnalysis::fgkHbtAnalyzeAll = 0;
-
-AliHBTAnalysis::AliHBTAnalysis()
+const UInt_t AliHBTAnalysis::fgkDefaultMixingInfo = 1000;
+const Int_t AliHBTAnalysis::fgkDefaultBufferSize = 5;
+
+AliHBTAnalysis::AliHBTAnalysis():
+ fProcEvent(0x0),
+ fReader(0x0),
+ fNTrackFunctions(0),
+ fNParticleFunctions(0),
+ fNParticleAndTrackFunctions(0),
+ fNTrackMonitorFunctions(0),
+ fNParticleMonitorFunctions(0),
+ fNParticleAndTrackMonitorFunctions(0),
+ fTrackFunctions ( new AliHBTOnePairFctn* [fgkFctnArraySize]),
+ fParticleFunctions ( new AliHBTOnePairFctn* [fgkFctnArraySize]),
+ fParticleAndTrackFunctions ( new AliHBTTwoPairFctn* [fgkFctnArraySize]),
+ fParticleMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),
+ fTrackMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),
+ fParticleAndTrackMonitorFunctions ( new AliHBTMonTwoParticleFctn* [fgkFctnArraySize]),
+ fBkgEventCut(0x0),
+ fPartBuffer(0x0),
+ fTrackBuffer(0x0),
+ fBufferSize(2),
+ fDisplayMixingInfo(fgkDefaultMixingInfo),
+ fIsOwner(kFALSE),
+ fProcessOption(kSimulatedAndReconstructed),
+ fNoCorrfctns(kFALSE),
+ fOutputFileName(0x0),
+ fVertexX(0.0),
+ fVertexY(0.0),
+ fVertexZ(0.0)
{
- fReader = 0x0;
+ //default constructor
- fTrackFunctions = new AliHBTOnePairFctn* [fgkFctnArraySize];
- fParticleFunctions = new AliHBTOnePairFctn* [fgkFctnArraySize];
- fParticleAndTrackFunctions = new AliHBTTwoPairFctn* [fgkFctnArraySize];
-
- fNTrackFunctions = 0;
- fNParticleFunctions = 0;
- fNParticleAndTrackFunctions = 0;
-
- fPairCut = new AliHBTEmptyPairCut();//empty cut - accepts all particles
-
- fBufferSize = 2;
}
/*************************************************************************************/
+AliHBTAnalysis::AliHBTAnalysis(const AliHBTAnalysis& in):
+ AliAnalysis(in),
+ fProcEvent(0x0),
+ fReader(0x0),
+ fNTrackFunctions(0),
+ fNParticleFunctions(0),
+ fNParticleAndTrackFunctions(0),
+ fNTrackMonitorFunctions(0),
+ fNParticleMonitorFunctions(0),
+ fNParticleAndTrackMonitorFunctions(0),
+ fTrackFunctions(0x0),
+ fParticleFunctions(0x0),
+ fParticleAndTrackFunctions(0x0),
+ fParticleMonitorFunctions(0x0),
+ fTrackMonitorFunctions(0x0),
+ fParticleAndTrackMonitorFunctions(0x0),
+ fBkgEventCut(0x0),
+ fPartBuffer(0x0),
+ fTrackBuffer(0x0),
+ fBufferSize(fgkDefaultBufferSize),
+ fDisplayMixingInfo(fgkDefaultMixingInfo),
+ fIsOwner(kFALSE),
+ fProcessOption(kSimulatedAndReconstructed),
+ fNoCorrfctns(kFALSE),
+ fOutputFileName(0x0),
+ fVertexX(0.0),
+ fVertexY(0.0),
+ fVertexZ(0.0)
+ {
+//copy constructor
+ Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
+ }
+/*************************************************************************************/
+AliHBTAnalysis& AliHBTAnalysis::operator=(const AliHBTAnalysis& /*right*/)
+ {
+//operator =
+ Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
+ return *this;
+ }
+/*************************************************************************************/
AliHBTAnalysis::~AliHBTAnalysis()
{
//destructor
//note that we do not delete functions itself
// they should be deleted by whom where created
//we only store pointers, and we use "new" only for pointers array
+
+ if (fIsOwner)
+ {
+ if (AliVAODParticle::GetDebug()>5)Info("~AliHBTAnalysis","Is Owner: Attempting to delete functions");
+ DeleteFunctions();
+ if (AliVAODParticle::GetDebug()>5)Info("~AliHBTAnalysis","Delete functions done");
+ }
delete [] fTrackFunctions;
delete [] fParticleFunctions;
delete [] fParticleAndTrackFunctions;
- delete fPairCut; // always have an copy of an object - we create we dstroy
+ delete [] fParticleMonitorFunctions;
+ delete [] fTrackMonitorFunctions;
+ delete [] fParticleAndTrackMonitorFunctions;
+
+ delete fBkgEventCut;
+ delete fOutputFileName;
}
/*************************************************************************************/
-void AliHBTAnalysis::Process(Option_t* option)
+Int_t AliHBTAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim)
{
- //default option = "TracksAndParticles"
- //Main method of the HBT Analysis Package
- //It triggers reading with the global cut (default is an empty cut)
- //Than it checks options and data which are read
- //if everything is OK, then it calls one of the looping methods
- //depending on tfReaderhe option
- //These methods differs on what thay are looping on
- //
- // METHOD OPTION
- //--------------------------------------------------------------------
- //ProcessTracksAndParticles - "TracksAndParticles"
- // DEFAULT
- // it loops over both, tracks(reconstructed) and particles(simulated)
- // all function gethered in all 3 lists are called for each (double)pair
- //
- //ProcessTracks - "Tracks"
- // it loops only on tracks(reconstructed),
- // functions ONLY from fTrackFunctions list are called
- //
- //ProcessParticles - "Particles"
- // it loops only on particles(simulated),
- // functions ONLY from fParticleAndTrackFunctions list are called
- //
- //
- if (!fReader)
- {
- Error("Process","The reader is not set");
- return;
- }
-
+ //Processes one event
+ if (fProcEvent == 0x0)
+ {
+ Error("ProcessEvent","Analysis <<%s>> option not specified.",GetName());
+ return 1;
+ }
+ if ( Rejected(aodrec,aodsim) )
+ {
+// Double_t x = 0, y = 0, z = 0;
+// aodrec->GetPrimaryVertex(x,y,z);
+// Info("ProcessEvent","Event has vertex at %f %f %f",x,y,z);
+ Info("ProcessEvent","Nch is %d",aodsim->GetNumberOfCharged());
+ Info("ProcessEvent","%s: Event cut rejected this event",GetName());
+ return 0;
+ }
+
+ //Move event to the apparent vertex -> must be after the event cut
+ //It is important for any cut that use any spacial coordiantes,
+ //f.g. anti-merging cut in order to preserve the same bahavior of variables (f.g. distance between tracks)
+ Double_t dvx = 0, dvy = 0, dvz = 0;
+ if (aodrec)
+ {
+ Double_t pvx,pvy,pvz;
+ aodrec->GetPrimaryVertex(pvx,pvy,pvz);
+
+ dvx = fVertexX - pvx;
+ dvy = fVertexY - pvy;
+ dvz = fVertexZ - pvz;
+ aodrec->Move(dvx,dvy,dvz);
+ }
+
+ Int_t result = (this->*fProcEvent)(aodrec,aodsim);
+
+ if (aodrec) aodrec->Move(-dvx,-dvy,-dvz);//move event back to the same spacial coordinates
+
+ return result;
+}
+/*************************************************************************************/
+
+Int_t AliHBTAnalysis::Finish()
+{
+//Finishes analysis
+ WriteFunctions();
+ return 0;
+}
+/*************************************************************************************/
+
+void AliHBTAnalysis::DeleteFunctions()
+{
+ //Deletes all functions added to analysis
- const char *oT = strstr(option,"Tracks");
- const char *oP = strstr(option,"Particles");
+ UInt_t ii;
+ for(ii = 0;ii<fNParticleFunctions;ii++)
+ {
+ if (AliVAODParticle::GetDebug()>5)
+ {
+ Info("DeleteFunctions","Deleting ParticleFunction %#x",fParticleFunctions[ii]);
+ Info("DeleteFunctions","Deleting ParticleFunction %s",fParticleFunctions[ii]->Name());
+ }
+ delete fParticleFunctions[ii];
+ }
+ fNParticleFunctions = 0;
+
+ for(ii = 0;ii<fNTrackFunctions;ii++)
+ {
+ if (AliVAODParticle::GetDebug()>5)
+ {
+ Info("DeleteFunctions","Deleting TrackFunction %#x",fTrackFunctions[ii]);
+ Info("DeleteFunctions","Deleting TrackFunction %s",fTrackFunctions[ii]->Name());
+ }
+ delete fTrackFunctions[ii];
+ }
+ fNTrackFunctions = 0;
- Bool_t nonid = IsNonIdentAnalysis();
+ for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
+ {
+ if (AliVAODParticle::GetDebug()>5)
+ {
+ Info("DeleteFunctions","Deleting ParticleAndTrackFunction %#x",fParticleAndTrackFunctions[ii]);
+ Info("DeleteFunctions","Deleting ParticleAndTrackFunction %s",fParticleAndTrackFunctions[ii]->Name());
+ }
+ delete fParticleAndTrackFunctions[ii];
+ }
+ fNParticleAndTrackFunctions = 0;
- if(oT && oP)
+ for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
{
- if (fReader->GetNumberOfPartEvents() <1)
+ if (AliVAODParticle::GetDebug()>5)
{
- Error("Process","There is no Particles. Maybe change the option?");
- return;
+ Info("DeleteFunctions","Deleting ParticleMonitorFunction %#x",fParticleMonitorFunctions[ii]);
+ Info("DeleteFunctions","Deleting ParticleMonitorFunction %s",fParticleMonitorFunctions[ii]->Name());
}
- if (fReader->GetNumberOfTrackEvents() <1)
+ delete fParticleMonitorFunctions[ii];
+ }
+ fNParticleMonitorFunctions = 0;
+
+ for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
+ {
+ if (AliVAODParticle::GetDebug()>5)
{
- Error("Process","There is no Tracks. Maybe change the option?");
- return;
+ Info("DeleteFunctions","Deleting TrackMonitorFunction %#x",fTrackMonitorFunctions[ii]);
+ Info("DeleteFunctions","Deleting TrackMonitorFunction %s",fTrackMonitorFunctions[ii]->Name());
}
-
- if ( RunCoherencyCheck() )
- {
- Error("Process",
- "Coherency check not passed. Maybe change the option?\n");
- return;
- }
- if (nonid) ProcessTracksAndParticlesNonIdentAnal();
- else ProcessTracksAndParticles();
- return;
- }
-
- if(oT)
- {
- if (fReader->GetNumberOfTrackEvents() <1)
+ delete fTrackMonitorFunctions[ii];
+ }
+ fNTrackMonitorFunctions = 0;
+
+ for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
+ {
+ if (AliVAODParticle::GetDebug()>5)
{
- Error("Process","There is no data to analyze.");
- return;
+ Info("DeleteFunctions","Deleting ParticleAndTrackMonitorFunction %#x",fParticleAndTrackMonitorFunctions[ii]);
+ Info("DeleteFunctions","Deleting ParticleAndTrackMonitorFunction %s",fParticleAndTrackMonitorFunctions[ii]->Name());
}
- if (nonid) ProcessTracksNonIdentAnal();
- else ProcessTracks();
- return;
- }
+ delete fParticleAndTrackMonitorFunctions[ii];
+ }
+ fNParticleAndTrackMonitorFunctions = 0;
- if(oP)
+}
+/*************************************************************************************/
+
+Int_t AliHBTAnalysis::Init()
+{
+//Initializeation method
+//calls Init for all functions
+
+ //Depending on option and pair cut assigns proper analysis method
+ Bool_t nonid = IsNonIdentAnalysis();
+ switch(fProcessOption)
{
- if (fReader->GetNumberOfPartEvents() <1)
- {
- Error("Process","There is no data to analyze.");
- return;
- }
- if (nonid) ProcessParticlesNonIdentAnal();
- else ProcessParticles();
-
-// cout<<"NON ID"<<endl;
-// ProcessParticlesNonIdentAnal();
-// cout<<"\n\n\n NORMAL"<<endl;
-// ProcessParticles();
-
- return;
+ case kReconstructed:
+ if (nonid) fProcEvent = &AliHBTAnalysis::ProcessRecNonId;
+ else fProcEvent = &AliHBTAnalysis::ProcessRec;
+ SetCutsOnRec();
+ break;
+
+ case kSimulated:
+ if (nonid) fProcEvent = &AliHBTAnalysis::ProcessSimNonId;
+ else fProcEvent = &AliHBTAnalysis::ProcessSim;
+ SetCutsOnSim();
+ break;
+
+ case kSimulatedAndReconstructed:
+ if (nonid) fProcEvent = &AliHBTAnalysis::ProcessRecAndSimNonId;
+ else fProcEvent = &AliHBTAnalysis::ProcessRecAndSim;
+ break;
}
+
+ if (fPartBuffer == 0x0) fPartBuffer = new AliEventBuffer (fBufferSize);
+ else fPartBuffer->Reset();
+
+ if (fTrackBuffer == 0x0) fTrackBuffer = new AliEventBuffer (fBufferSize);
+ else fTrackBuffer->Reset();
+
+
+ fNoCorrfctns = (fNParticleFunctions == 0) && (fNTrackFunctions == 0) && (fNParticleAndTrackFunctions == 0);
-}
+ UInt_t ii;
+ for(ii = 0;ii<fNParticleFunctions;ii++)
+ fParticleFunctions[ii]->Init();
+
+ for(ii = 0;ii<fNTrackFunctions;ii++)
+ fTrackFunctions[ii]->Init();
+ for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
+ fParticleAndTrackFunctions[ii]->Init();
+
+ for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
+ fParticleMonitorFunctions[ii]->Init();
+
+ for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
+ fTrackMonitorFunctions[ii]->Init();
+
+ for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
+ fParticleAndTrackMonitorFunctions[ii]->Init();
+
+ return 0;
+}
/*************************************************************************************/
-void AliHBTAnalysis::ProcessTracksAndParticles()
+void AliHBTAnalysis::ResetFunctions()
{
-
-//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
-//the loops are splited
+//In case fOwner is true, deletes all functions
+//in other case, just set number of analysis to 0
+ if (fIsOwner) DeleteFunctions();
+ else
+ {
+ fNParticleFunctions = 0;
+ fNTrackFunctions = 0;
+ fNParticleAndTrackFunctions = 0;
+ fNParticleMonitorFunctions = 0;
+ fNTrackMonitorFunctions = 0;
+ fNParticleAndTrackMonitorFunctions = 0;
+ }
+}
+/*************************************************************************************/
+Int_t AliHBTAnalysis::ProcessRecAndSim(AliAOD* aodrec, AliAOD* aodsim)
+{
+//Does analysis for both tracks and particles
+//mainly for resolution study and analysies with weighting algirithms
+// cut on particles only -- why?
+// - PID: when we make resolution analysis we want to take only tracks with correct PID
+// We need cut on tracks because there are data characteristic
- AliHBTParticle * part1, * part2;
- AliHBTParticle * track1, * track2;
+ AliVAODParticle * part1, * part2;
+ AliVAODParticle * track1, * track2;
- AliHBTEvent * trackEvent, *partEvent;
- AliHBTEvent * trackEvent2,*partEvent2;
+ AliAOD * trackEvent = aodrec, *partEvent = aodsim;
+ AliAOD* trackEvent1 = new AliAOD();
+ AliAOD* partEvent1 = new AliAOD();
+
+ AliAOD * trackEvent2,*partEvent2;
// Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
- Int_t Nev = fReader->GetNumberOfTrackEvents();
-
- /***************************************/
- /****** Looping same events ********/
- /****** filling numerators ********/
- /***************************************/
- AliHBTPair * trackpair = new AliHBTPair();
- AliHBTPair * partpair = new AliHBTPair();
+// Int_t nev = fReader->GetNumberOfTrackEvents();
+ static AliHBTPair tpair;
+ static AliHBTPair ppair;
+ AliHBTPair* trackpair = &tpair;
+ AliHBTPair* partpair = &ppair;
+
AliHBTPair * tmptrackpair;//temprary pointers to pairs
AliHBTPair * tmppartpair;
+ register UInt_t ii;
- for (Int_t i = 0;i<Nev;i++)
- {
- partEvent= fReader->GetParticleEvent(i);
- trackEvent = fReader->GetTrackEvent(i);
-
- if (!partEvent) continue;
-
- //N = 0;
-
- for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
- {
- if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
-
- part1= partEvent->GetParticle(j);
- track1= trackEvent->GetParticle(j);
-
- for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
+
+ if ( !partEvent || !trackEvent )
+ {
+ Error("ProcessRecAndSim","<<%s>> Can not get event",GetName());
+ return 1;
+ }
+
+ if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
+ {
+ Error("ProcessRecAndSim",
+ "Number of simulated particles (%d) not equal to number of reconstructed tracks (%d). Skipping Event.",
+ partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
+ return 2;
+ }
+
+
+ for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
+ {
+ /***************************************/
+ /****** Looping same events ********/
+ /****** filling numerators ********/
+ /***************************************/
+ if ( (j%fDisplayMixingInfo) == 0)
+ Info("ProcessRecAndSim",
+ "Mixing particle %d with particles from the same event",j);
+
+ part1= partEvent->GetParticle(j);
+ track1= trackEvent->GetParticle(j);
+
+ Bool_t firstcut = (this->*fkPass1)(part1,track1);
+ if (fBufferSize != 0)
+ if ( (firstcut == kFALSE) || ( (this->*fkPass2)(part1,track1) == kFALSE ) )
+ {
+ //accepted by any cut
+ // we have to copy because reader keeps only one event
+
+ partEvent1->AddParticle(part1);
+ trackEvent1->AddParticle(track1);
+ }
+
+ if (firstcut) continue;
+
+ for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
+ fParticleMonitorFunctions[ii]->Process(part1);
+ for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
+ fTrackMonitorFunctions[ii]->Process(track1);
+ for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
+ fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
+
+ if (fNoCorrfctns) continue;
+
+ for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
+ {
+ part2= partEvent->GetParticle(k);
+ if (part1->GetUID() == part2->GetUID()) continue;
+ partpair->SetParticles(part1,part2);
+
+ track2= trackEvent->GetParticle(k);
+ trackpair->SetParticles(track1,track2);
+
+ if( (this->*fkPass)(partpair,trackpair) ) //check pair cut
+ { //do not meets crietria of the pair cut, try with swapped pairs
+ if( (this->*fkPass)((AliHBTPair*)partpair->GetSwappedPair(),(AliHBTPair*)trackpair->GetSwappedPair()) )
+ continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
+ else
+ { //swaped pair meets all the criteria
+ tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
+ tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
+ }
+ }
+ else
+ {//meets criteria of the pair cut
+ tmptrackpair = trackpair;
+ tmppartpair = partpair;
+ }
+
+ for(ii = 0;ii<fNParticleFunctions;ii++)
+ fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
+
+ for(ii = 0;ii<fNTrackFunctions;ii++)
+ fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
+
+ for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
+ fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair,tmppartpair);
+ //end of 2nd loop over particles from the same event
+ }//for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
+
+ /***************************************/
+ /***** Filling denominators *********/
+ /***************************************/
+ if (fBufferSize == 0) continue;
+
+ fPartBuffer->ResetIter();
+ fTrackBuffer->ResetIter();
+ Int_t m = 0;
+ while (( partEvent2 = fPartBuffer->Next() ))
+ {
+ trackEvent2 = fTrackBuffer->Next();
+
+ m++;
+ if ( (j%fDisplayMixingInfo) == 0)
+ Info("ProcessRecAndSim",
+ "Mixing particle %d from current event with particles from event %d",j,-m);
+
+ for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
{
- part2= partEvent->GetParticle(k);
+ part2= partEvent2->GetParticle(l);
partpair->SetParticles(part1,part2);
-
- track2= trackEvent->GetParticle(k);
+
+ track2= trackEvent2->GetParticle(l);
trackpair->SetParticles(track1,track2);
- if(fPairCut->Pass(partpair) || (fPairCut->Pass(trackpair))) //check pair cut
- { //do not meets crietria of the pair cut, try with swaped pairs
- if( ( fPairCut->Pass(partpair->GetSwapedPair()) ) || ( fPairCut->Pass(trackpair->GetSwapedPair()) ) )
- continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
+ if( (this->*fkPass)(partpair,trackpair) ) //check pair cut
+ { //do not meets crietria of the
+ if( (this->*fkPass)((AliHBTPair*)partpair->GetSwappedPair(),(AliHBTPair*)trackpair->GetSwappedPair()) )
+ continue;
else
- { //swaped pair meets all the criteria
- tmppartpair = partpair->GetSwapedPair();
- tmptrackpair = trackpair->GetSwapedPair();
-
+ {
+ tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
+ tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
}
}
else
{//meets criteria of the pair cut
- tmptrackpair = trackpair;
- tmppartpair = partpair;
+ tmptrackpair = trackpair;
+ tmppartpair = partpair;
}
- UInt_t ii;
+
for(ii = 0;ii<fNParticleFunctions;ii++)
- fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
-
+ fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
+
for(ii = 0;ii<fNTrackFunctions;ii++)
- fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
-
+ fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
+
for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
- fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair,tmppartpair);
- }
- }
- }
+ fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair);
+ }//for(Int_t l = 0; l<N2;l++) // ... on all particles
+ }
+ //end of loop over particles from first event
+ }//for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
+ delete fPartBuffer->Push(partEvent1);
+ delete fTrackBuffer->Push(trackEvent1);
+ //end of loop over events
+ return 0;
+}
+/*************************************************************************************/
- /***** Filling denominators *********/
- /***************************************/
- for (Int_t i = 0;i<Nev-1;i++) //In each event (but last) ....
- {
-
- partEvent= fReader->GetParticleEvent(i);
- if (!partEvent) continue;
-
- trackEvent = fReader->GetTrackEvent(i);
-
-// N=0;
-
- for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
- {
-
- part1= partEvent->GetParticle(j);
+Int_t AliHBTAnalysis::ProcessSim(AliAOD* /*aodrec*/, AliAOD* aodsim)
+{
+ //Does analysis of simulated data
+ AliVAODParticle * part1, * part2;
+
+ AliAOD* partEvent = aodsim;
+ AliAOD* partEvent1 = new AliAOD();
- track1= trackEvent->GetParticle(j);
-
- Int_t NNN;
-
- if ( ((i+fBufferSize) >= Nev) ||( fBufferSize < 0) ) //if buffer size is negative
- //or current event+buffersize is greater
- //than max nuber of events
- {
- NNN = Nev; //set the max event number
- }
- else
- {
- NNN = i+fBufferSize; //set the current event number + fBufferSize
- }
+ AliAOD* partEvent2;
+
+ AliHBTPair ppair;
+
+ AliHBTPair* partpair = &ppair;
- for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
+ AliHBTPair * tmppartpair;
+
+ register UInt_t ii;
+
+ if ( !partEvent )
+ {
+ Error("ProcessSim","Can not get event");
+ return 1;
+ }
+
+
+ for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
+ {
+ /***************************************/
+ /****** Looping same events ********/
+ /****** filling numerators ********/
+ /***************************************/
+ if ( (j%fDisplayMixingInfo) == 0)
+ Info("ProcessSim",
+ "Mixing particle %d with particles from the same event",j);
+
+ part1= partEvent->GetParticle(j);
+
+ Bool_t firstcut = fPairCut->GetFirstPartCut()->Rejected(part1);
+
+ if (fBufferSize != 0)
+ if ( (firstcut == kFALSE) || ( fPairCut->GetSecondPartCut()->Rejected(part1) == kFALSE ) )
+ {
+ //accepted by any cut
+ // we have to copy because reader keeps only one event
+
+ partEvent1->AddParticle(part1);
+ }
+
+ if (firstcut) continue;
+
+ for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
+ fParticleMonitorFunctions[ii]->Process(part1);
+
+ if ( fNParticleFunctions == 0 ) continue;
+
+ for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
+ {
+ part2= partEvent->GetParticle(k);
+ if (part1->GetUID() == part2->GetUID()) continue;
+ partpair->SetParticles(part1,part2);
+
+ if(fPairCut->Rejected(partpair)) //check pair cut
+ { //do not meets crietria of the
+ if( fPairCut->Rejected((AliHBTPair*)partpair->GetSwappedPair()) ) continue;
+ else tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
+ }
+ else
+ {
+ tmppartpair = partpair;
+ }
+
+ for(ii = 0;ii<fNParticleFunctions;ii++)
+ fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
+
+ //end of 2nd loop over particles from the same event
+ }//for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
+
+ /***************************************/
+ /***** Filling denominators *********/
+ /***************************************/
+ if (fBufferSize == 0) continue;
+
+ fPartBuffer->ResetIter();
+ Int_t m = 0;
+ while (( partEvent2 = fPartBuffer->Next() ))
+ {
+ m++;
+ if ( (j%fDisplayMixingInfo) == 0)
+ Info("ProcessSim",
+ "Mixing particle %d from current event with particles from event %d",j,-m);
+ for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
{
-
- partEvent2= fReader->GetParticleEvent(k);
- if (!partEvent2) continue;
-
- trackEvent2 = fReader->GetTrackEvent(k);
-
- if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
-
- for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
- {
-
- // if (N>MAXCOMB) break;
-
- part2= partEvent2->GetParticle(l);
- partpair->SetParticles(part1,part2);
-
- track2= trackEvent2->GetParticle(l);
- trackpair->SetParticles(track1,track2);
-
- if( (fPairCut->Pass(partpair)) || (fPairCut->Pass(trackpair)) ) //check pair cut
- { //do not meets crietria of the
- if( ( fPairCut->Pass(partpair->GetSwapedPair()) ) || ( fPairCut->Pass(trackpair->GetSwapedPair()) ) )
- continue;
- else
- {
- tmppartpair = partpair->GetSwapedPair();
- tmptrackpair = trackpair->GetSwapedPair();
- }
- }
- else
- {//meets criteria of the pair cut
- tmptrackpair = trackpair;
- tmppartpair = partpair;
- }
- UInt_t ii;
- for(ii = 0;ii<fNParticleFunctions;ii++)
- fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
-
- for(ii = 0;ii<fNTrackFunctions;ii++)
- fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
-
- for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
- fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair);
-
-
- }//for(Int_t l = 0; l<N2;l++) // ... on all particles
- }//for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
- }
-
- }
- /***************************************/
-
-
-}
-/*************************************************************************************/
+ part2= partEvent2->GetParticle(l);
+ partpair->SetParticles(part1,part2);
+
+ if( fPairCut->Rejected(partpair) ) //check pair cut
+ { //do not meets crietria of the
+ if( fPairCut->Rejected((AliHBTPair*)partpair->GetSwappedPair()) )
+ continue;
+ else
+ {
+ tmppartpair = (AliHBTPair*)partpair->GetSwappedPair();
+ }
+ }
+ else
+ {//meets criteria of the pair cut
+ tmppartpair = partpair;
+ }
-void AliHBTAnalysis::ProcessTracks()
+ for(ii = 0;ii<fNParticleFunctions;ii++)
+ fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
+
+ }//for(Int_t l = 0; l<N2;l++) // ... on all particles
+ }
+ }
+
+
+ delete fPartBuffer->Push(partEvent1);
+ //end of loop over events
+ return 0;
+}
+/*************************************************************************************/
+Int_t AliHBTAnalysis::ProcessRec(AliAOD* aodrec, AliAOD* /*aodsim*/)
{
- //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
-//the loops are splited
- AliHBTParticle * track1, * track2;
-
- AliHBTEvent * trackEvent;
- AliHBTEvent * trackEvent2;
-
-// Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
+ //Does analysis of reconstructed data
+ AliVAODParticle * track1, * track2;
- Int_t Nev = fReader->GetNumberOfTrackEvents();
+ AliAOD* trackEvent = aodrec;
+ AliAOD* trackEvent1 = new AliAOD();
+
+ AliAOD* trackEvent2;
- /***************************************/
- /****** Looping same events ********/
- /****** filling numerators ********/
- /***************************************/
- AliHBTPair * trackpair = new AliHBTPair();
- AliHBTPair * tmptrackpair; //temporary pointer
+ AliHBTPair tpair;
+
+ AliHBTPair* trackpair = &tpair;
+
+ AliHBTPair * tmptrackpair;
- for (Int_t i = 0;i<Nev;i++)
- {
- trackEvent = fReader->GetTrackEvent(i);
- if (!trackEvent) continue;
- //N = 0;
-
- for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
- {
- if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
+ register UInt_t ii;
- track1= trackEvent->GetParticle(j);
-
- for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
- {
- track2= trackEvent->GetParticle(k);
- trackpair->SetParticles(track1,track2);
- if(fPairCut->Pass(trackpair)) //check pair cut
+ if ( !trackEvent )
+ {
+ Error("ProcessRecAndSim","Can not get event");
+ return 1;
+ }
+
+
+ for (Int_t j = 0; j<trackEvent->GetNumberOfParticles() ; j++)
+ {
+ /***************************************/
+ /****** Looping same events ********/
+ /****** filling numerators ********/
+ /***************************************/
+ if ( (j%fDisplayMixingInfo) == 0)
+ Info("ProcessRec",
+ "Mixing Particle %d with Particles from the same event",j);
+
+ track1= trackEvent->GetParticle(j);
+
+ Bool_t firstcut = fPairCut->GetFirstPartCut()->Rejected(track1);
+
+ if (fBufferSize != 0)
+ if ( (firstcut == kFALSE) || ( fPairCut->GetSecondPartCut()->Rejected(track1) == kFALSE ) )
+ {
+ //accepted by any cut
+ // we have to copy because reader keeps only one event
+
+ trackEvent1->AddParticle(track1);
+ }
+
+ if (firstcut) continue;
+
+ for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
+ fParticleMonitorFunctions[ii]->Process(track1);
+
+ if ( fNTrackFunctions == 0 ) continue;
+
+ for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
+ {
+ track2= trackEvent->GetParticle(k);
+ if (track1->GetUID() == track2->GetUID()) continue;
+ trackpair->SetParticles(track1,track2);
+
+ if(fPairCut->Rejected(trackpair)) //check pair cut
+ { //do not meets crietria of the
+ if( fPairCut->Rejected((AliHBTPair*)trackpair->GetSwappedPair()) ) continue;
+ else tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
+ }
+ else
+ {
+ tmptrackpair = trackpair;
+ }
+
+ for(ii = 0;ii<fNTrackFunctions;ii++)
+ fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
+
+ //end of 2nd loop over Particles from the same event
+ }//for (Int_t k =j+1; k < trackEvent->GetNumberOfParticles() ; k++)
+
+ /***************************************/
+ /***** Filling denominators *********/
+ /***************************************/
+ if (fBufferSize == 0) continue;
+
+ fTrackBuffer->ResetIter();
+ Int_t m = 0;
+ while (( trackEvent2 = fTrackBuffer->Next() ))
+ {
+ m++;
+ if ( (j%fDisplayMixingInfo) == 0)
+ Info("ProcessRec",
+ "Mixing Particle %d from current event with Particles from event %d",j,-m);
+ for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++) // ... on all Particles
+ {
+
+ track2= trackEvent2->GetParticle(l);
+ trackpair->SetParticles(track1,track2);
+
+ if( fPairCut->Rejected(trackpair) ) //check pair cut
{ //do not meets crietria of the
- if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
- else tmptrackpair = trackpair->GetSwapedPair();
+ if( fPairCut->Rejected((AliHBTPair*)trackpair->GetSwappedPair()) )
+ continue;
+ else
+ {
+ tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair();
+ }
}
else
- {
- tmptrackpair = trackpair;
- }
- UInt_t ii;
-
+ {//meets criteria of the pair cut
+ tmptrackpair = trackpair;
+ }
+
for(ii = 0;ii<fNTrackFunctions;ii++)
- fTrackFunctions[ii]->ProcessSameEventParticles(tmptrackpair);
-
-
- }
+ fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
+
+ }//for(Int_t l = 0; l<N2;l++) // ... on all Particles
+ }
+ }
+ delete fTrackBuffer->Push(trackEvent1);
+ //end of loop over events
+ return 0;
+}
+/*************************************************************************************/
+
+Int_t AliHBTAnalysis::ProcessRecAndSimNonId(AliAOD* aodrec, AliAOD* aodsim)
+{
+//Analyzes both reconstructed and simulated data
+ if (aodrec == 0x0)
+ {
+ Error("ProcessTracksAndParticlesNonIdentAnal","Reconstructed event is NULL");
+ return 1;
+ }
+
+ if (aodsim == 0x0)
+ {
+ Error("ProcessTracksAndParticlesNonIdentAnal","Simulated event is NULL");
+ return 1;
+ }
+
+ if ( aodrec->GetNumberOfParticles() != aodsim->GetNumberOfParticles() )
+ {
+ Error("ProcessTracksAndParticlesNonIdentAnal",
+ "Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
+ aodsim->GetNumberOfParticles() , aodrec->GetNumberOfParticles());
+ return 2;
+ }
+
+
+ AliVAODParticle * part1, * part2;
+ AliVAODParticle * track1, * track2;
+
+ static AliAOD aodrec1;
+ static AliAOD aodsim1;
+
+ AliAOD * trackEvent1=&aodrec1,*partEvent1=&aodsim1;//Particle that passes first particle cut, this event
+ trackEvent1->Reset();
+ partEvent1->Reset();
+ AliAOD * trackEvent2=0x0,*partEvent2=0x0;//Particle that passes second particle cut, this event
+ AliAOD * trackEvent3=0x0,*partEvent3=0x0;//Particle that passes second particle cut, events from buffer
+
+ AliAOD* rawtrackEvent = aodrec;//this we get from Reader
+ AliAOD* rawpartEvent = aodsim;//this we get from Reader
+
+ static AliHBTPair tpair;
+ static AliHBTPair ppair;
+
+ AliHBTPair* trackpair = &tpair;
+ AliHBTPair* partpair = &ppair;
+
+ register UInt_t ii;
+
+ /********************************/
+ /* Filtering out */
+ /********************************/
+ if ( ( (partEvent2==0x0) || (trackEvent2==0x0)) )//in case fBufferSize == 0 and pointers are created do not eneter
+ {
+ partEvent2 = new AliAOD();
+ trackEvent2 = new AliAOD();
+ }
+
+ FilterOut(partEvent1, partEvent2, rawpartEvent, trackEvent1, trackEvent2, rawtrackEvent);
+
+ for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
+ {
+ if ( (j%fDisplayMixingInfo) == 0)
+ Info("ProcessTracksAndParticlesNonIdentAnal",
+ "Mixing particle %d from current event with particles from current event",j);
+
+ part1= partEvent1->GetParticle(j);
+ track1= trackEvent1->GetParticle(j);
+
+
+ for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
+ fParticleMonitorFunctions[ii]->Process(part1);
+ for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
+ fTrackMonitorFunctions[ii]->Process(track1);
+ for(ii = 0; ii<fNParticleAndTrackMonitorFunctions; ii++)
+ fParticleAndTrackMonitorFunctions[ii]->Process(track1,part1);
+
+ if (fNoCorrfctns) continue;
+
+ /***************************************/
+ /****** filling numerators ********/
+ /****** (particles from event2) ********/
+ /***************************************/
+
+ for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups
+ {
+ part2= partEvent2->GetParticle(k);
+ if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
+ partpair->SetParticles(part1,part2);
+
+ track2= trackEvent2->GetParticle(k);
+ trackpair->SetParticles(track1,track2);
+
+ if( (this->*fkPassPairProp)(partpair,trackpair) ) //check pair cut
+ { //do not meets crietria of the pair cut
+ continue;
+ }
+ else
+ {//meets criteria of the pair cut
+ for(ii = 0;ii<fNParticleFunctions;ii++)
+ fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
+
+ for(ii = 0;ii<fNTrackFunctions;ii++)
+ fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
+
+ for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
+ fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(trackpair,partpair);
+ }
}
- }
- /***************************************/
- /***** Filling diff histogram *********/
- /***************************************/
- for (Int_t i = 0;i<Nev-1;i++) //In each event (but last) ....
+ if ( fBufferSize == 0) continue;//do not mix diff histograms
+ /***************************************/
+ /***** Filling denominators *********/
+ /***************************************/
+ fPartBuffer->ResetIter();
+ fTrackBuffer->ResetIter();
+
+ Int_t nmonitor = 0;
+
+ while ( (partEvent3 = fPartBuffer->Next() ) != 0x0)
{
- trackEvent = fReader->GetTrackEvent(i);
- if (!trackEvent) continue;
-// N=0;
-
- for (Int_t j = 0; j< trackEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
- {
-// if (N>MAXCOMB) break;
-
- track1= trackEvent->GetParticle(j);
+ trackEvent3 = fTrackBuffer->Next();
+
+ if ( (j%fDisplayMixingInfo) == 0)
+ Info("ProcessTracksAndParticlesNonIdentAnal",
+ "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
+
+ for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
+ {
+ part2= partEvent3->GetParticle(k);
+ partpair->SetParticles(part1,part2);
+
+ track2= trackEvent3->GetParticle(k);
+ trackpair->SetParticles(track1,track2);
+
+ if( (this->*fkPassPairProp)(partpair,trackpair) ) //check pair cut
+ { //do not meets crietria of the pair cut
+ continue;
+ }
+ else
+ {//meets criteria of the pair cut
+ UInt_t ii;
+ for(ii = 0;ii<fNParticleFunctions;ii++)
+ fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
+
+ for(ii = 0;ii<fNTrackFunctions;ii++)
+ fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
+
+ for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
+ fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(trackpair,partpair);
+ }
+ }// for particles event2
+ }//while event2
+ }//for over particles in event1
+
+ delete fPartBuffer->Push(partEvent2);
+ delete fTrackBuffer->Push(trackEvent2);
- Int_t NNN;
-
- if ( ((i+fBufferSize) >= Nev) ||( fBufferSize < 0) ) //if buffer size is negative
- //or current event+buffersize is greater
- //than max nuber of events
- {
- NNN = Nev; //set the max event number
- }
- else
- {
- NNN = i+fBufferSize; //set the current event number + fBufferSize
- }
+ return 0;
+}
+/*************************************************************************************/
+Int_t AliHBTAnalysis::ProcessSimNonId(AliAOD* /*aodrec*/, AliAOD* aodsim)
+{
+//does analysis of simulated (MC) data in non-identical mode
+//i.e. when particles selected by first part. cut are a disjunctive set than particles
+//passed by the second part. cut
+ if (aodsim == 0x0)
+ {
+ return 1;
+ }
- for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
- {
-
- trackEvent2 = fReader->GetTrackEvent(k);
- if (!trackEvent2) continue;
-
- if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
-
- for(Int_t l = 0; l<trackEvent2->GetNumberOfParticles();l++) // ... on all particles
- {
-
- // if (N>MAXCOMB) break;
- track2= trackEvent2->GetParticle(l);
- trackpair->SetParticles(track1,track2);
-
- if(fPairCut->Pass(trackpair)) //check pair cut
- { //do not meets crietria of the
- if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
- else tmptrackpair = trackpair->GetSwapedPair();
- }
- else
- {
- tmptrackpair = trackpair;
- }
- UInt_t ii;
- for(ii = 0;ii<fNTrackFunctions;ii++)
- fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
-
- }//for(Int_t l = 0; l<N2;l++) // ... on all particles
- }//for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
- }
-
- }
+
+ AliVAODParticle * part1, * part2;
- /***************************************/
-
+ static AliAOD aodsim1;
-}
+ AliAOD* partEvent1=&aodsim1;//Particle that passes first particle cut, this event
+ partEvent1->Reset();
+ AliAOD* partEvent2=0x0;//Particle that passes second particle cut, this event
+ AliAOD* partEvent3=0x0;//Particle that passes second particle cut, events from buffer
-/*************************************************************************************/
+ AliAOD* rawpartEvent = aodsim;//this we get from Reader
-void AliHBTAnalysis::ProcessParticles()
-{
- //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
-//the loops are splited
- AliHBTParticle * part1, * part2;
+ static AliHBTPair ppair;
- AliHBTEvent * partEvent;
- AliHBTEvent * partEvent2;
+ AliHBTPair* partpair = &ppair;
- AliHBTPair * partpair = new AliHBTPair();
- AliHBTPair * tmppartpair; //temporary pointer to the pair
+ register UInt_t ii;
-// Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
+ /********************************/
+ /* Filtering out */
+ /********************************/
+ if (partEvent2==0x0)//in case fBufferSize == 0 and pointers are created do not eneter
+ {
+ partEvent2 = new AliAOD();
+ }
+
+ FilterOut(partEvent1, partEvent2, rawpartEvent);
+
+ for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
+ {
+ if ( (j%fDisplayMixingInfo) == 0)
+ Info("ProcessParticlesNonIdentAnal",
+ "Mixing particle %d from current event with particles from current event",j);
+
+ part1= partEvent1->GetParticle(j);
+
+
+ for(ii = 0; ii<fNParticleMonitorFunctions; ii++)
+ fParticleMonitorFunctions[ii]->Process(part1);
+
+ if (fNParticleFunctions == 0) continue;
+
+ /***************************************/
+ /****** filling numerators ********/
+ /****** (particles from event2) ********/
+ /***************************************/
+
+ for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups
+ {
+ part2= partEvent2->GetParticle(k);
+ if (part1->GetUID() == part2->GetUID()) continue;//this is the same particle but with different PID
+ partpair->SetParticles(part1,part2);
+
+
+ if(fPairCut->PassPairProp(partpair) ) //check pair cut
+ { //do not meets crietria of the pair cut
+ continue;
+ }
+ else
+ {//meets criteria of the pair cut
+ for(ii = 0;ii<fNParticleFunctions;ii++)
+ fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
+ }
+ }
+
+ if ( fBufferSize == 0) continue;//do not mix diff histograms
+ /***************************************/
+ /***** Filling denominators *********/
+ /***************************************/
+ fPartBuffer->ResetIter();
+
+ Int_t nmonitor = 0;
+
+ while ( (partEvent3 = fPartBuffer->Next() ) != 0x0)
+ {
+
+ if ( (j%fDisplayMixingInfo) == 0)
+ Info("ProcessParticlesNonIdentAnal",
+ "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
+
+ for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
+ {
+ part2= partEvent3->GetParticle(k);
+ partpair->SetParticles(part1,part2);
+
+
+ if(fPairCut->PassPairProp(partpair) ) //check pair cut
+ { //do not meets crietria of the pair cut
+ continue;
+ }
+ else
+ {//meets criteria of the pair cut
+ for(ii = 0;ii<fNParticleFunctions;ii++)
+ {
+ fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
+ }
+ }
+ }// for particles event2
+ }//while event2
+ }//for over particles in event1
+
+ delete fPartBuffer->Push(partEvent2);
+
+ return 0;
+}
+/*************************************************************************************/
+Int_t AliHBTAnalysis::ProcessRecNonId(AliAOD* aodrec, AliAOD* /*aodsim*/)
+{
+//Analyzes both reconstructed and simulated data
+ if (aodrec == 0x0)
+ {
+ return 1;
+ }
+
+ AliVAODParticle * track1, * track2;
+
+ static AliAOD aodrec1;
+ AliAOD * trackEvent1=&aodrec1;//Particle that passes first particle cut, this event
+ trackEvent1->Reset();
+ AliAOD * trackEvent2=0x0;//Particle that passes second particle cut, this event
+ AliAOD * trackEvent3=0x0;//Particle that passes second particle cut, events from buffer
+ AliAOD* rawtrackEvent = aodrec;//this we get from Reader
+
+ static AliHBTPair tpair;
- Int_t Nev = fReader->GetNumberOfPartEvents();
+ AliHBTPair* trackpair = &tpair;
+
+ register UInt_t ii;
- // Nev = 1;
- /***************************************/
- /****** Looping same events ********/
- /****** filling numerators ********/
- /***************************************/
-// gBenchmark->Start("id");
- for (Int_t i = 0;i<Nev;i++)
- {
- partEvent= fReader->GetParticleEvent(i);
- if (!partEvent) continue;
- //N = 0;
-
- for (Int_t j = 0; j<partEvent->GetNumberOfParticles() ; j++)
- {
- if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
+ /********************************/
+ /* Filtering out */
+ /********************************/
+ if ( trackEvent2==0x0 )//in case fBufferSize == 0 and pointers are created do not eneter
+ {
+ trackEvent2 = new AliAOD();
+ }
- part1= partEvent->GetParticle(j);
-
- for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++)
- {
- part2= partEvent->GetParticle(k);
- partpair->SetParticles(part1,part2);
-
- if( fPairCut->Pass(partpair) ) //check pair cut
- { //do not meets crietria of the pair cut, try with swaped pairs
- if( fPairCut->Pass(partpair->GetSwapedPair() ) )
- continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
- else
- { //swaped pair meets all the criteria
- tmppartpair = partpair->GetSwapedPair();
- }
- }
- else
- {
- tmppartpair = partpair;
- }
+ FilterOut(trackEvent1, trackEvent2, rawtrackEvent);
+
+ for (Int_t j = 0; j<trackEvent1->GetNumberOfParticles() ; j++)
+ {
+ if ( (j%fDisplayMixingInfo) == 0)
+ Info("ProcessTracksNonIdentAnal",
+ "Mixing particle %d from current event with particles from current event",j);
+
+ track1= trackEvent1->GetParticle(j);
+
+
+ for(ii = 0; ii<fNTrackMonitorFunctions; ii++)
+ fTrackMonitorFunctions[ii]->Process(track1);
+
+ if (fNTrackFunctions == 0x0) continue;
+
+ /***************************************/
+ /****** filling numerators ********/
+ /****** (particles from event2) ********/
+ /***************************************/
+
+ for (Int_t k = 0; k < trackEvent2->GetNumberOfParticles() ; k++) //partEvent1 and partEvent2 are particles from the same event but separated to two groups
+ {
+ track2= trackEvent2->GetParticle(k);
+ if (track1->GetUID() == track2->GetUID()) continue;//this is the same particle but with different PID
+ trackpair->SetParticles(track1,track2);
+
+
+ if( fPairCut->PassPairProp(trackpair)) //check pair cut
+ { //do not meets crietria of the pair cut
+ continue;
+ }
+ else
+ {//meets criteria of the pair cut
+ UInt_t ii;
+ for(ii = 0;ii<fNTrackFunctions;ii++)
+ fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
+ }
+ }
+
+ if ( fBufferSize == 0) continue;//do not mix diff histograms
+ /***************************************/
+ /***** Filling denominators *********/
+ /***************************************/
+ fTrackBuffer->ResetIter();
+
+ Int_t nmonitor = 0;
+
+ while ( (trackEvent3 = fTrackBuffer->Next() ) != 0x0)
+ {
+ if ( (j%fDisplayMixingInfo) == 0)
+ Info("ProcessTracksNonIdentAnal",
+ "Mixing particle %d from current event with particles from event%d",j,-(++nmonitor));
+
+ for (Int_t k = 0; k < trackEvent3->GetNumberOfParticles() ; k++)
+ {
+ track2= trackEvent3->GetParticle(k);
+ trackpair->SetParticles(track1,track2);
+
+ if( fPairCut->PassPairProp(trackpair)) //check pair cut
+ { //do not meets crietria of the pair cut
+ continue;
+ }
+ else
+ {//meets criteria of the pair cut
+ for(ii = 0;ii<fNTrackFunctions;ii++)
+ fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
+ }
+ }// for particles event2
+ }//while event2
+ }//for over particles in event1
+
+ delete fTrackBuffer->Push(trackEvent2);
+
+ return 0;
+}
+/*************************************************************************************/
+
+void AliHBTAnalysis::Process(Option_t* option)
+{
+ //default option = "TracksAndParticles"
+ //Main method of the HBT Analysis Package
+ //It triggers reading with the global cut (default is an empty cut)
+ //Than it checks options and data which are read
+ //if everything is OK, then it calls one of the looping methods
+ //depending on tfReaderhe option
+ //These methods differs on what thay are looping on
+ //
+ // METHOD OPTION
+ //--------------------------------------------------------------------
+ //ProcessTracksAndParticles - "TracksAndParticles"
+ // DEFAULT
+ // it loops over both, tracks(reconstructed) and particles(simulated)
+ // all function gethered in all 3 lists are called for each (double)pair
+ //
+ //ProcessTracks - "Tracks"
+ // it loops only on tracks(reconstructed),
+ // functions ONLY from fTrackFunctions list are called
+ //
+ //ProcessParticles - "Particles"
+ // it loops only on particles(simulated),
+ // functions ONLY from fParticleAndTrackFunctions list are called
+ //
+ //
+ if (!fReader)
+ {
+ Error("Process","The reader is not set");
+ return;
+ }
+
+ const char *oT = strstr(option,"Tracks");
+ const char *oP = strstr(option,"Particles");
+
+ Bool_t nonid = IsNonIdentAnalysis();
+
+ Init();
+
+ if(oT && oP)
+ {
+ if (nonid) ProcessTracksAndParticlesNonIdentAnal();
+ else ProcessTracksAndParticles();
+ return;
+ }
+
+ if(oT)
+ {
+ if (nonid) ProcessTracksNonIdentAnal();
+ else ProcessTracks();
+ return;
+ }
+
+ if(oP)
+ {
+ if (nonid) ProcessParticlesNonIdentAnal();
+ else ProcessParticles();
+ return;
+ }
+
+}
+/*************************************************************************************/
- UInt_t ii;
-
- for(ii = 0;ii<fNParticleFunctions;ii++)
- fParticleFunctions[ii]->ProcessSameEventParticles(tmppartpair);
- }
- }
- }
+void AliHBTAnalysis::ProcessTracksAndParticles()
+{
+//Makes analysis for both tracks and particles
+//mainly for resolution study and analysies with weighting algirithms
+//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
+//the loops are splited
+
+// cut on particles only -- why?
+// - PID: when we make resolution analysis we want to take only tracks with correct PID
+// We need cut on tracks because there are data characteristic to
+
+ AliAOD * trackEvent, *partEvent;
- /***************************************/
- /***** Filling diff histogram *********/
- /***************************************/
- for (Int_t i = 0;i<Nev-1;i++) //In each event (but last)....
+ fReader->Rewind();
+ while (fReader->Next() == kFALSE)
{
- partEvent= fReader->GetParticleEvent(i);
- if (!partEvent) continue;
+ partEvent = fReader->GetEventSim();
+ trackEvent = fReader->GetEventRec();
+ ProcessRecAndSim(trackEvent,partEvent);
-// N=0;
-
- for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
- {
-// if (N>MAXCOMB) break;
-
- part1= partEvent->GetParticle(j);
-
- Int_t NNN;
-
- if ( ((i+fBufferSize) >= Nev) ||( fBufferSize < 0) ) //if buffer size is negative
- //or current event+buffersize is greater
- //than max nuber of events
- {
- NNN = Nev; //set the max event number
- }
- else
- {
- NNN = i+fBufferSize; //set the current event number + fBufferSize
- }
-
-// cout<<"NNN = "<<NNN<<endl;
- for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
- {
-
- partEvent2= fReader->GetParticleEvent(k);
- if (!partEvent2) continue;
-
- if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<k<<endl;
-
- for(Int_t l = 0; l<partEvent2->GetNumberOfParticles();l++) // ... on all particles
- {
-
- // if (N>MAXCOMB) break;
- part2= partEvent2->GetParticle(l);
- partpair->SetParticles(part1,part2);
-
- if(fPairCut->Pass(partpair)) //check pair cut
- { //do not meets crietria of the
- if( fPairCut->Pass(partpair->GetSwapedPair()) ) continue;
- else tmppartpair = partpair->GetSwapedPair();
- }
- else
- {
- tmppartpair = partpair;
- }
- UInt_t ii;
- for(ii = 0;ii<fNParticleFunctions;ii++)
- fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
-
- }//for(Int_t l = 0; l<N2;l++) // ... on all particles
- }//for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
- }
+ }//while (fReader->Next() == kFALSE)
- }
-
- /***************************************/
-
-// gBenchmark->Stop("id");
-// gBenchmark->Show("id");
+}
+/*************************************************************************************/
+void AliHBTAnalysis::ProcessTracks()
+{
+//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
+//the loops are splited
+ AliAOD * trackEvent;
+ fReader->Rewind();
+ while (fReader->Next() == kFALSE)
+ {
+ trackEvent = fReader->GetEventRec();
+ ProcessRec(trackEvent,0x0);
+ }//while (fReader->Next() == kFALSE)
}
/*************************************************************************************/
+void AliHBTAnalysis::ProcessParticles()
+{
+//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
+//the loops are splited
+ AliAOD * partEvent;
+ fReader->Rewind();
+ while (fReader->Next() == kFALSE)
+ {
+ partEvent = fReader->GetEventSim();
+ ProcessSim(0x0,partEvent);
+ }//while (fReader->Next() == kFALSE)
+}
+/*************************************************************************************/
void AliHBTAnalysis::WriteFunctions()
{
+//Calls Write for all defined functions in analysis
+//== writes all results
+ TFile* oututfile = 0x0;
+ if (fOutputFileName)
+ {
+ oututfile = TFile::Open(*fOutputFileName,"update");
+ }
UInt_t ii;
for(ii = 0;ii<fNParticleFunctions;ii++)
- fParticleFunctions[ii]->Write();
+ {
+ if (AliVAODParticle::GetDebug()>5)
+ {
+ Info("WriteFunctions","Writing ParticleFunction %#x",fParticleFunctions[ii]);
+ Info("WriteFunctions","Writing ParticleFunction %s",fParticleFunctions[ii]->Name());
+ }
+ fParticleFunctions[ii]->Write();
+ }
for(ii = 0;ii<fNTrackFunctions;ii++)
- fTrackFunctions[ii]->Write();
+ {
+ if (AliVAODParticle::GetDebug()>5)
+ {
+ Info("WriteFunctions","Writing TrackFunction %#x",fTrackFunctions[ii]);
+ Info("WriteFunctions","Writing TrackFunction %s",fTrackFunctions[ii]->Name());
+ }
+ fTrackFunctions[ii]->Write();
+ }
for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
- fParticleAndTrackFunctions[ii]->Write();
+ {
+ if (AliVAODParticle::GetDebug()>5)
+ {
+ Info("WriteFunctions","Writing ParticleAndTrackFunction %#x",fParticleAndTrackFunctions[ii]);
+ Info("WriteFunctions","Writing ParticleAndTrackFunction %s",fParticleAndTrackFunctions[ii]->Name());
+ }
+ fParticleAndTrackFunctions[ii]->Write();
+ }
+
+ for(ii = 0;ii<fNParticleMonitorFunctions;ii++)
+ {
+ if (AliVAODParticle::GetDebug()>5)
+ {
+ Info("WriteFunctions","Writing ParticleMonitorFunction %#x",fParticleMonitorFunctions[ii]);
+ Info("WriteFunctions","Writing ParticleMonitorFunction %s",fParticleMonitorFunctions[ii]->Name());
+ }
+ fParticleMonitorFunctions[ii]->Write();
+ }
+
+ for(ii = 0;ii<fNTrackMonitorFunctions;ii++)
+ {
+ if (AliVAODParticle::GetDebug()>5)
+ {
+ Info("WriteFunctions","Writing TrackMonitorFunction %#x",fTrackMonitorFunctions[ii]);
+ Info("WriteFunctions","Writing TrackMonitorFunction %s",fTrackMonitorFunctions[ii]->Name());
+ }
+ fTrackMonitorFunctions[ii]->Write();
+ }
+
+ for(ii = 0;ii<fNParticleAndTrackMonitorFunctions;ii++)
+ {
+ if (AliVAODParticle::GetDebug()>5)
+ {
+ Info("WriteFunctions","Writing ParticleAndTrackMonitorFunction %#x",fParticleAndTrackMonitorFunctions[ii]);
+ Info("WriteFunctions","Writing ParticleAndTrackMonitorFunction %s",fParticleAndTrackMonitorFunctions[ii]->Name());
+ }
+ fParticleAndTrackMonitorFunctions[ii]->Write();
+ }
+ delete oututfile;
+}
+/*************************************************************************************/
+
+void AliHBTAnalysis::SetOutputFileName(const char* fname)
+{
+ //Sets fiele name where to dump results,
+ //if not specified reults are written to gDirectory
+ if (fname == 0x0)
+ {
+ delete fOutputFileName;
+ fOutputFileName = 0x0;
+ return;
+ }
+ if ( strcmp(fname,"") == 0 )
+ {
+ delete fOutputFileName;
+ fOutputFileName = 0x0;
+ return;
+ }
+ if (fOutputFileName == 0x0) fOutputFileName = new TString(fname);
+ else *fOutputFileName = fname;
}
/*************************************************************************************/
-void AliHBTAnalysis::SetGlobalPairCut(AliHBTPairCut* cut)
+void AliHBTAnalysis::SetGlobalPairCut(AliAODPairCut* cut)
{
+//Sets the global cut
if (cut == 0x0)
{
Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
}
delete fPairCut;
- fPairCut = (AliHBTPairCut*)cut->Clone();
+ fPairCut = (AliAODPairCut*)cut->Clone();
}
/*************************************************************************************/
void AliHBTAnalysis::AddTrackFunction(AliHBTOnePairFctn* f)
{
+//Adds track function
if (f == 0x0) return;
if (fNTrackFunctions == fgkFctnArraySize)
{
void AliHBTAnalysis::AddParticleFunction(AliHBTOnePairFctn* f)
{
+//adds particle function
if (f == 0x0) return;
if (fNParticleFunctions == fgkFctnArraySize)
}
fParticleFunctions[fNParticleFunctions] = f;
fNParticleFunctions++;
-
-
}
+/*************************************************************************************/
+
void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTTwoPairFctn* f)
{
+//add resolution function
if (f == 0x0) return;
if (fNParticleAndTrackFunctions == fgkFctnArraySize)
{
fParticleAndTrackFunctions[fNParticleAndTrackFunctions] = f;
fNParticleAndTrackFunctions++;
}
+/*************************************************************************************/
+
+void AliHBTAnalysis::AddParticleMonitorFunction(AliHBTMonOneParticleFctn* f)
+{
+//add particle monitoring function
+ if (f == 0x0) return;
+
+ if (fNParticleMonitorFunctions == fgkFctnArraySize)
+ {
+ Error("AliHBTAnalysis::AddParticleMonitorFunction","Can not add this function, not enough place in the array.");
+ }
+ fParticleMonitorFunctions[fNParticleMonitorFunctions] = f;
+ fNParticleMonitorFunctions++;
+}
+/*************************************************************************************/
+void AliHBTAnalysis::AddTrackMonitorFunction(AliHBTMonOneParticleFctn* f)
+{
+//add track monitoring function
+ if (f == 0x0) return;
+ if (fNTrackMonitorFunctions == fgkFctnArraySize)
+ {
+ Error("AliHBTAnalysis::AddTrackMonitorFunction","Can not add this function, not enough place in the array.");
+ }
+ fTrackMonitorFunctions[fNTrackMonitorFunctions] = f;
+ fNTrackMonitorFunctions++;
+}
/*************************************************************************************/
+void AliHBTAnalysis::AddParticleAndTrackMonitorFunction(AliHBTMonTwoParticleFctn* f)
+{
+//add resolution monitoring function
+ if (f == 0x0) return;
+ if (fNParticleAndTrackMonitorFunctions == fgkFctnArraySize)
+ {
+ Error("AliHBTAnalysis::AddParticleAndTrackMonitorFunction","Can not add this function, not enough place in the array.");
+ }
+ fParticleAndTrackMonitorFunctions[fNParticleAndTrackMonitorFunctions] = f;
+ fNParticleAndTrackMonitorFunctions++;
+}
+
+/*************************************************************************************/
/*************************************************************************************/
Bool_t AliHBTAnalysis::RunCoherencyCheck()
//Checks if both HBTRuns are similar
//return true if error found
//if they seem to be OK return false
+/*
+
Int_t i;
- cout<<"Checking HBT Runs Coherency"<<endl;
+ Info("RunCoherencyCheck","Checking HBT Runs Coherency");
+
+//When we use non-buffering reader this is a big waste of time -> We need to read all data to check it
+//and reader is implemented safe in this case anyway
+// Info("RunCoherencyCheck","Number of events ...");
+// if (fReader->GetNumberOfPartEvents() == fReader->GetNumberOfTrackEvents() ) //check whether there is the same number of events
+// {
+// Info("RunCoherencyCheck","OK. %d found\n",fReader->GetNumberOfTrackEvents());
+// }
+// else
+// { //if not the same - ERROR
+// Error("RunCoherencyCheck",
+// "Number of simulated events (%d) is not equal to number of reconstructed events(%d)",
+// fReader->GetNumberOfPartEvents(),fReader->GetNumberOfTrackEvents());
+// return kTRUE;
+// }
- cout<<"Number of events ..."; fflush(stdout);
+ Info("RunCoherencyCheck","Checking number of Particles AND Particles Types in each event ...");
- if (fReader->GetNumberOfPartEvents() == fReader->GetNumberOfTrackEvents() ) //check whether there is the same number of events
- {
- cout<<"OK. "<<fReader->GetNumberOfTrackEvents()<<" found"<<endl;
- }
- else
- { //if not the same - ERROR
- Error("AliHBTAnalysis::RunCoherencyCheck()",
- "Number of simulated events (%d) is not equal to number of reconstructed events(%d)",
- fReader->GetNumberOfPartEvents(),fReader->GetNumberOfTrackEvents());
- return kTRUE;
- }
-
- cout<<"Checking number of Particles AND Particles Types in each event ...";fflush(stdout);
-
- AliHBTEvent *partEvent;
- AliHBTEvent *trackEvent;
+ AliAOD *partEvent;
+ AliAOD *trackEvent;
for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
{
- partEvent= fReader->GetParticleEvent(i); //gets the "ith" event
- trackEvent = fReader->GetTrackEvent(i);
+ partEvent= fReader->GetEventSim(i); //gets the "ith" event
+ trackEvent = fReader->GetEventRec(i);
if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
if ( (partEvent == 0x0) || (partEvent == 0x0) )
{
- Error("AliHBTAnalysis::RunCoherencyCheck()",
+ Error("RunCoherencyCheck",
"One event is NULL and the other one not. Event Number %d",i);
return kTRUE;
}
if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
{
- Error("AliHBTAnalysis::RunCoherencyCheck()",
+ Error("RunCoherencyCheck",
"Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
return kTRUE;
{
if( partEvent->GetParticle(j)->GetPdgCode() != trackEvent->GetParticle(j)->GetPdgCode() )
{
- Error("AliHBTAnalysis::RunCoherencyCheck()",
+ Error("RunCoherencyCheck",
"Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)",
i,j, partEvent->GetParticle(j)->GetPdgCode(),trackEvent->GetParticle(j)->GetPdgCode() );
return kTRUE;
}
}
}
- cout<<" Done"<<endl;
- cout<<" Everything looks OK"<<endl;
+ Info("RunCoherencyCheck"," Done");
+ Info("RunCoherencyCheck"," Everything looks OK");
+*/
return kFALSE;
}
void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
{
- AliHBTParticle * part1, * part2;
- AliHBTParticle * track1, * track2;
-
- AliHBTEvent * trackEvent1,*partEvent1;
- AliHBTEvent * trackEvent2,*partEvent2;
- AliHBTEvent * trackEvent3,*partEvent3;
-
- AliHBTEvent * rawtrackEvent, *rawpartEvent;
-// AliHBTEvent * rawtrackEvent2,*rawpartEvent2;
-
-
- Int_t Nev = fReader->GetNumberOfTrackEvents();
-
- AliHBTPair * trackpair = new AliHBTPair();
- AliHBTPair * partpair = new AliHBTPair();
-
- TList tbuffer;
- TList pbuffer;
- Int_t ninbuffer = 0;
+//Performs analysis for both, tracks and particles
+ AliAOD* rawtrackEvent, * rawpartEvent;
+ fReader->Rewind();
+
+ Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
+ Info("ProcessTracksAndParticlesNonIdentAnal","***** NON IDENT MODE ****************");
+ Info("ProcessTracksAndParticlesNonIdentAnal","**************************************");
- trackEvent1 = new AliHBTEvent();
- partEvent1 = new AliHBTEvent();
- trackEvent1->SetOwner(kFALSE);
- partEvent1->SetOwner(kFALSE);;
-
- cout<<"**************************************"<<endl;
- cout<<"***** NON IDENT MODE ****************"<<endl;
- cout<<"**************************************"<<endl;
- for (Int_t i = 0;i<Nev;i++)
+ for (Int_t i = 0;;i++)//infinite loop
{
- rawpartEvent = fReader->GetParticleEvent(i);
- rawtrackEvent = fReader->GetTrackEvent(i);
- if ( (rawpartEvent == 0x0) || (rawtrackEvent == 0x0) ) continue;//in case of any error
-
- /********************************/
- /* Filtering out */
- /********************************/
- if (ninbuffer > fBufferSize)
- {
- partEvent2 = (AliHBTEvent*)pbuffer.Remove(pbuffer.Last()); //remove from the end to be reset, filled and put on the beginning
- trackEvent2 = (AliHBTEvent*)tbuffer.Remove(tbuffer.Last());
- ninbuffer--;
- }
- else
- {
- partEvent2 = new AliHBTEvent();
- trackEvent2 = new AliHBTEvent();
- partEvent2->SetOwner(kFALSE);
- trackEvent2->SetOwner(kFALSE);
- }
- FilterOut(partEvent1, partEvent2, rawpartEvent, trackEvent1, trackEvent2, rawtrackEvent);
+ if (fReader->Next()) break; //end when no more events available
- for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
- {
- if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
-
- part1= partEvent1->GetParticle(j);
- track1= trackEvent1->GetParticle(j);
-
-
- /***************************************/
- /****** filling numerators ********/
- /****** (particles from event2) ********/
- /***************************************/
- for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++)
- {
- part2= partEvent2->GetParticle(k);
- partpair->SetParticles(part1,part2);
-
- track2= trackEvent2->GetParticle(k);
- trackpair->SetParticles(track1,track2);
-
-
- if( (fPairCut->PassPairProp(partpair)) || (fPairCut->PassPairProp(trackpair))) //check pair cut
- { //do not meets crietria of the pair cut
- continue;
- }
- else
- {//meets criteria of the pair cut
- UInt_t ii;
- for(ii = 0;ii<fNParticleFunctions;ii++)
- fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
-
- for(ii = 0;ii<fNTrackFunctions;ii++)
- fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
-
- for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
- fParticleAndTrackFunctions[ii]->ProcessSameEventParticles(trackpair,partpair);
- }
- }
- /***************************************/
- /***** Filling denominators *********/
- /***************************************/
- TIter piter(&pbuffer);
- TIter titer(&tbuffer);
- Int_t nmonitor = 0;
-
- while ( (partEvent3 = (AliHBTEvent*)piter.Next()) != 0x0)
- {
- trackEvent3 = (AliHBTEvent*)titer.Next();
- if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i - (++nmonitor)<<endl;
-
- for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
- {
- part2= partEvent3->GetParticle(k);
- partpair->SetParticles(part1,part2);
-
- track2= trackEvent3->GetParticle(k);
- trackpair->SetParticles(track1,track2);
-
-
- if( (fPairCut->PassPairProp(partpair)) || (fPairCut->PassPairProp(trackpair))) //check pair cut
- { //do not meets crietria of the pair cut
- continue;
- }
- else
- {//meets criteria of the pair cut
- UInt_t ii;
- for(ii = 0;ii<fNParticleFunctions;ii++)
- fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
-
- for(ii = 0;ii<fNTrackFunctions;ii++)
- fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
-
- for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
- fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(trackpair,partpair);
- }
- }// for particles event2
- }//while event2
- }//for over particles in event1
-
- pbuffer.AddFirst(partEvent2);
- tbuffer.AddFirst(trackEvent2);
- ninbuffer++;
-
+ rawpartEvent = fReader->GetEventSim();
+ rawtrackEvent = fReader->GetEventRec();
+
+ ProcessRecAndSimNonId(rawtrackEvent,rawpartEvent);
}//end of loop over events (1)
-
- pbuffer.SetOwner(); //to clean stored events by the way of buffer destruction
- tbuffer.SetOwner();
}
/*************************************************************************************/
void AliHBTAnalysis::ProcessTracksNonIdentAnal()
{
- AliHBTParticle * track1, * track2;
-
- AliHBTEvent * trackEvent1;
- AliHBTEvent * trackEvent2;
- AliHBTEvent * trackEvent3;
-
- AliHBTEvent * rawtrackEvent;
-// AliHBTEvent * rawtrackEvent2,*rawpartEvent2;
-
+//Process Tracks only with non identical mode
+ AliAOD * rawtrackEvent;
+ fReader->Rewind();
- Int_t Nev = fReader->GetNumberOfTrackEvents();
-
- AliHBTPair * trackpair = new AliHBTPair();
-
- TList tbuffer;
- Int_t ninbuffer = 0;
-
- trackEvent1 = new AliHBTEvent();
- trackEvent1->SetOwner(kFALSE);
+ Info("ProcessTracksNonIdentAnal","**************************************");
+ Info("ProcessTracksNonIdentAnal","***** NON IDENT MODE ****************");
+ Info("ProcessTracksNonIdentAnal","**************************************");
- cout<<"**************************************"<<endl;
- cout<<"***** NON IDENT MODE ****************"<<endl;
- cout<<"**************************************"<<endl;
- for (Int_t i = 0;i<Nev;i++)
+ for (Int_t i = 0;;i++)//infinite loop
{
- rawtrackEvent = fReader->GetTrackEvent(i);
- if (rawtrackEvent == 0x0) continue;//in case of any error
-
- /********************************/
- /* Filtering out */
- /********************************/
- if (ninbuffer > fBufferSize)
- {
- trackEvent2 = (AliHBTEvent*)tbuffer.Remove(tbuffer.Last());
- ninbuffer--;
- }
- else
- {
- trackEvent2 = new AliHBTEvent();
- trackEvent2->SetOwner(kFALSE);
- }
- FilterOut(trackEvent1, trackEvent2, rawtrackEvent);
-
- for (Int_t j = 0; j<trackEvent1->GetNumberOfParticles() ; j++)
- {
- if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
-
- track1= trackEvent1->GetParticle(j);
-
-
- /***************************************/
- /****** filling numerators ********/
- /****** (particles from event2) ********/
- /***************************************/
- for (Int_t k = 0; k < trackEvent2->GetNumberOfParticles() ; k++)
- {
- track2= trackEvent2->GetParticle(k);
- trackpair->SetParticles(track1,track2);
-
-
- if( fPairCut->PassPairProp(trackpair)) //check pair cut
- { //do not meets crietria of the pair cut
- continue;
- }
- else
- {//meets criteria of the pair cut
- UInt_t ii;
- for(ii = 0;ii<fNTrackFunctions;ii++)
- fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
- }
- }
- /***************************************/
- /***** Filling denominators *********/
- /***************************************/
- TIter titer(&tbuffer);
- Int_t nmonitor = 0;
-
- while ( (trackEvent3 = (AliHBTEvent*)titer.Next()) != 0x0)
- {
-
- if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i - (++nmonitor)<<endl;
-
- for (Int_t k = 0; k < trackEvent3->GetNumberOfParticles() ; k++)
- {
-
- track2= trackEvent3->GetParticle(k);
- trackpair->SetParticles(track1,track2);
-
-
- if( fPairCut->PassPairProp(trackpair)) //check pair cut
- { //do not meets crietria of the pair cut
- continue;
- }
- else
- {//meets criteria of the pair cut
- UInt_t ii;
- for(ii = 0;ii<fNTrackFunctions;ii++)
- fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
-
- }
- }// for particles event2
- }//while event2
- }//for over particles in event1
-
- tbuffer.AddFirst(trackEvent2);
- ninbuffer++;
-
+ if (fReader->Next()) break; //end when no more events available
+ rawtrackEvent = fReader->GetEventRec();
+ ProcessRecNonId(rawtrackEvent,0x0);
}//end of loop over events (1)
-
- tbuffer.SetOwner();
}
/*************************************************************************************/
void AliHBTAnalysis::ProcessParticlesNonIdentAnal()
{
- AliHBTParticle * part1 = 0x0, * part2 = 0x0;
-
- AliHBTEvent * partEvent1 = 0x0;
- AliHBTEvent * partEvent2 = 0x0;
- AliHBTEvent * partEvent3 = 0x0;
-
- AliHBTEvent * rawpartEvent = 0x0;
-
- Int_t Nev = fReader->GetNumberOfPartEvents();
-
- AliHBTPair * partpair = new AliHBTPair();
-
- TList pbuffer;
- Int_t ninbuffer = 0;
-
- partEvent1 = new AliHBTEvent();
- partEvent1->SetOwner(kFALSE);;
+//process paricles only with non identical mode
+ AliAOD * rawpartEvent = 0x0;
+ fReader->Rewind();
- cout<<"**************************************"<<endl;
- cout<<"***** PART NON IDENT MODE ******"<<endl;
- cout<<"**************************************"<<endl;
+ Info("ProcessParticlesNonIdentAnal","**************************************");
+ Info("ProcessParticlesNonIdentAnal","***** NON IDENT MODE ****************");
+ Info("ProcessParticlesNonIdentAnal","**************************************");
-// gBenchmark->Start("non_id");
- for (Int_t i = 0;i<Nev;i++)
+ for (Int_t i = 0;;i++)//infinite loop
{
- rawpartEvent = fReader->GetParticleEvent(i);
- if ( rawpartEvent == 0x0 ) continue;//in case of any error
-
- /********************************/
- /* Filtering out */
- /********************************/
- if (ninbuffer > fBufferSize)
- {
- partEvent2 = (AliHBTEvent*)pbuffer.Remove(pbuffer.Last()); //remove from the end to be reset, filled and put on the beginning
- ninbuffer--;
- }
- else
- {
- partEvent2 = new AliHBTEvent();
- partEvent2->SetOwner(kFALSE);
- }
- FilterOut(partEvent1, partEvent2, rawpartEvent);
+ if (fReader->Next()) break; //end when no more events available
- for (Int_t j = 0; j<partEvent1->GetNumberOfParticles() ; j++)
- {
- if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
-
- part1= partEvent1->GetParticle(j);
-
-
- /***************************************/
- /****** filling numerators ********/
- /****** (particles from event2) ********/
- /***************************************/
- for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++)
- {
- part2= partEvent2->GetParticle(k);
- partpair->SetParticles(part1,part2);
-
- if(fPairCut->PassPairProp(partpair) ) //check pair cut
- { //do not meets crietria of the pair cut
- continue;
- }
- else
- {//meets criteria of the pair cut
- UInt_t ii;
- for(ii = 0;ii<fNParticleFunctions;ii++)
- fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
-
- }
- }
- /***************************************/
- /***** Filling denominators *********/
- /***************************************/
- TIter piter(&pbuffer);
- Int_t nmonitor = 0;
-
- while ( (partEvent3 = (AliHBTEvent*)piter.Next()) != 0x0)
- {
- if ( (j%100) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i - (++nmonitor)<<endl;
-
- for (Int_t k = 0; k < partEvent3->GetNumberOfParticles() ; k++)
- {
- part2= partEvent3->GetParticle(k);
- partpair->SetParticles(part1,part2);
-
-
- if(fPairCut->PassPairProp(partpair) ) //check pair cut
- { //do not meets crietria of the pair cut
- continue;
- }
- else
- {//meets criteria of the pair cut
- UInt_t ii;
- for(ii = 0;ii<fNParticleFunctions;ii++)
- fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
-
- }
- }// for particles event2
- }//while event2
- }//for over particles in event1
-
- pbuffer.AddFirst(partEvent2);
- ninbuffer++;
-
+ rawpartEvent = fReader->GetEventSim();
+ ProcessSimNonId(0x0,rawpartEvent);
}//end of loop over events (1)
-
-// gBenchmark->Stop("non_id");
-// gBenchmark->Show("non_id");
- pbuffer.SetOwner();//to delete stered events.
}
/*************************************************************************************/
-void AliHBTAnalysis::FilterOut(AliHBTEvent* outpart1, AliHBTEvent* outpart2, AliHBTEvent* inpart,
- AliHBTEvent* outtrack1, AliHBTEvent* outtrack2, AliHBTEvent* intrack)
+void AliHBTAnalysis::FilterOut(AliAOD* outpart1, AliAOD* outpart2, AliAOD* inpart,
+ AliAOD* outtrack1, AliAOD* outtrack2, AliAOD* intrack) const
{
//Puts particles accepted as a first particle by global cut in out1
//and as a second particle in out2
- AliHBTParticle* part, *track;
+ AliVAODParticle* part, *track;
outpart1->Reset();
outpart2->Reset();
outtrack1->Reset();
outtrack2->Reset();
- AliHBTParticleCut *cut1 = fPairCut->GetFirstPartCut();
- AliHBTParticleCut *cut2 = fPairCut->GetSecondPartCut();
-
Bool_t in1, in2;
for (Int_t i = 0; i < inpart->GetNumberOfParticles(); i++)
part = inpart->GetParticle(i);
track = intrack->GetParticle(i);
- if ( (cut1->Pass(part)) || ( cut1->Pass(track) ) ) in1 = kFALSE; //if any, part or track, is rejected by cut1, in1 is false
- if ( (cut2->Pass(part)) || ( cut2->Pass(track) ) ) in2 = kFALSE; //if any, part or track, is rejected by cut2, in2 is false
+ if ( ((this->*fkPass1)(part,track)) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
+ if ( ((this->*fkPass2)(part,track)) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
if (gDebug)//to be removed in real analysis
if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
continue;
}
}
-
}
-
-
/*************************************************************************************/
-void AliHBTAnalysis::FilterOut(AliHBTEvent* out1, AliHBTEvent* out2, AliHBTEvent* in)
+
+void AliHBTAnalysis::FilterOut(AliAOD* out1, AliAOD* out2, AliAOD* in) const
{
//Puts particles accepted as a first particle by global cut in out1
//and as a second particle in out2
- AliHBTParticle* part;
+ AliVAODParticle* part;
out1->Reset();
out2->Reset();
- AliHBTParticleCut *cut1 = fPairCut->GetFirstPartCut();
- AliHBTParticleCut *cut2 = fPairCut->GetSecondPartCut();
+ AliAODParticleCut *cut1 = fPairCut->GetFirstPartCut();
+ AliAODParticleCut *cut2 = fPairCut->GetSecondPartCut();
Bool_t in1, in2;
in1 = in2 = kTRUE;
part = in->GetParticle(i);
- if ( cut1->Pass(part) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
- if ( cut2->Pass(part) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
+ if ( cut1->Rejected(part) ) in1 = kFALSE; //if part is rejected by cut1, in1 is false
+ if ( cut2->Rejected(part) ) in2 = kFALSE; //if part is rejected by cut2, in2 is false
if (gDebug)//to be removed in real analysis
if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
//in the future is possible to perform more sophisticated check
//if cuts have excluding requirements
- if (fPairCut->IsEmpty()) return kFALSE;
- if (fPairCut->GetFirstPartCut()->IsEmpty()) return kFALSE;
- if (fPairCut->GetSecondPartCut()->IsEmpty()) return kFALSE;
+ if (fPairCut->IsEmpty())
+ return kFALSE;
+
+ if (fPairCut->GetFirstPartCut()->IsEmpty())
+ return kFALSE;
+
+ if (fPairCut->GetSecondPartCut()->IsEmpty())
+ return kFALSE;
Int_t id1 = fPairCut->GetFirstPartCut()->GetPID();
Int_t id2 = fPairCut->GetSecondPartCut()->GetPID();
-
- if ( (id1==0) || (id2==0) || (id1==id2) ) return kFALSE;
+ if ( (id1==0) || (id2==0) || (id1==id2) )
+ return kFALSE;
+
return kTRUE;
}
+/*************************************************************************************/
+
+void AliHBTAnalysis::SetApparentVertex(Double_t x, Double_t y, Double_t z)
+{
+ //Sets apparent vertex
+ // All events have to be moved to the same vertex position in order to
+ // to be able to comare any space positions (f.g. anti-merging)
+ // This method defines this position
+
+ fVertexX = x;
+ fVertexY = y;
+ fVertexZ = z;
+}
+/*************************************************************************************/
+
+void AliHBTAnalysis::PressAnyKey()
+{
+ //small utility function that helps to make comfortable macros
+ char c;
+ int nread = -1;
+ fcntl(0, F_SETFL, O_NONBLOCK);
+ ::Info("","Press Any Key to continue ...");
+ while (nread<1)
+ {
+ nread = read(0, &c, 1);
+ gSystem->ProcessEvents();
+ }
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff