X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=HBTAN%2FAliHBTAnalysis.cxx;h=a216f8b2554282f84e2de96b61b35d2762920e23;hb=55d460483b6227b2554af960eb2672b4a03c43fd;hp=943a40d67a5e81fc486070708a81a0cbbd18280d;hpb=bed069a4fecf2a9734ae6be09bee7f5d2bbb852b;p=u%2Fmrichter%2FAliRoot.git diff --git a/HBTAN/AliHBTAnalysis.cxx b/HBTAN/AliHBTAnalysis.cxx index 943a40d67a5..a216f8b2554 100644 --- a/HBTAN/AliHBTAnalysis.cxx +++ b/HBTAN/AliHBTAnalysis.cxx @@ -1,35 +1,41 @@ #include "AliHBTAnalysis.h" -#include "AliHBTRun.h" -#include "AliHBTEvent.h" -#include "AliHBTReader.h" -#include "AliHBTParticle.h" -#include "AliHBTParticleCut.h" -#include "AliHBTPair.h" -#include "AliHBTPairCut.h" -#include "AliHBTFunction.h" -#include "AliHBTMonitorFunction.h" -#include "AliHBTEventBuffer.h" - -#include -#include - //_________________________________________________________ -/////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////// +// +// class AliHBTAnalysis // -//Central Object Of HBTAnalyser: -//This class performs main looping within HBT Analysis -//User must plug a reader of Type AliHBTReader -//User plugs in coorelation and monitor functions -//as well as monitor functions +// 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: +// 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 +// Piotr.Skowronski@cern.ch +// more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html // +//////////////////////////////////////////////////////////////////////////// //_________________________________________________________ + +#include +#include + +#include "AliAOD.h" +#include "AliAODParticle.h" +#include "AliAODPairCut.h" +#include "AliEventCut.h" + +#include "AliEventBuffer.h" + +#include "AliReader.h" +#include "AliHBTPair.h" +#include "AliHBTFunction.h" +#include "AliHBTMonitorFunction.h" + + ClassImp(AliHBTAnalysis) const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100; @@ -37,6 +43,7 @@ const UInt_t AliHBTAnalysis::fgkDefaultMixingInfo = 1000; const Int_t AliHBTAnalysis::fgkDefaultBufferSize = 5; AliHBTAnalysis::AliHBTAnalysis(): + fProcEvent(0x0), fReader(0x0), fNTrackFunctions(0), fNParticleFunctions(0), @@ -44,25 +51,33 @@ AliHBTAnalysis::AliHBTAnalysis(): 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) + fIsOwner(kFALSE), + fProcessOption(kSimulatedAndReconstructed), + fNoCorrfctns(kFALSE), + fOutputFileName(0x0), + fVertexX(0.0), + fVertexY(0.0), + fVertexZ(0.0) { -//default constructor - fTrackFunctions = new AliHBTOnePairFctn* [fgkFctnArraySize]; - fParticleFunctions = new AliHBTOnePairFctn* [fgkFctnArraySize]; - fParticleAndTrackFunctions = new AliHBTTwoPairFctn* [fgkFctnArraySize]; + //default constructor - fParticleMonitorFunctions = new AliHBTMonOneParticleFctn* [fgkFctnArraySize]; - fTrackMonitorFunctions = new AliHBTMonOneParticleFctn* [fgkFctnArraySize]; - fParticleAndTrackMonitorFunctions = new AliHBTMonTwoParticleFctn* [fgkFctnArraySize]; - - fPairCut = new AliHBTEmptyPairCut();//empty cut - accepts all particles } /*************************************************************************************/ AliHBTAnalysis::AliHBTAnalysis(const AliHBTAnalysis& in): - TObject(in), + AliAnalysis(in), + fProcEvent(0x0), fReader(0x0), fNTrackFunctions(0), fNParticleFunctions(0), @@ -76,16 +91,24 @@ AliHBTAnalysis::AliHBTAnalysis(const AliHBTAnalysis& in): fParticleMonitorFunctions(0x0), fTrackMonitorFunctions(0x0), fParticleAndTrackMonitorFunctions(0x0), - fPairCut(0x0), + fBkgEventCut(0x0), + fPartBuffer(0x0), + fTrackBuffer(0x0), fBufferSize(fgkDefaultBufferSize), fDisplayMixingInfo(fgkDefaultMixingInfo), - fIsOwner(kFALSE) + fIsOwner(kFALSE), + fProcessOption(kSimulatedAndReconstructed), + fNoCorrfctns(kFALSE), + fOutputFileName(0x0), + fVertexX(0.0), + fVertexY(0.0), + fVertexZ(0.0) { //copy constructor Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless"); } /*************************************************************************************/ -const AliHBTAnalysis& AliHBTAnalysis::operator=(const AliHBTAnalysis& /*right*/) +AliHBTAnalysis& AliHBTAnalysis::operator=(const AliHBTAnalysis& /*right*/) { //operator = Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless"); @@ -98,55 +121,186 @@ AliHBTAnalysis::~AliHBTAnalysis() //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) - DeleteFunctions(); + + 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 [] fParticleMonitorFunctions; delete [] fTrackMonitorFunctions; - delete [] fParticleAndTrackMonitorFunctions; + delete [] fParticleAndTrackMonitorFunctions; - delete fPairCut; // always have an copy of an object - we create we dstroy + delete fBkgEventCut; + delete fOutputFileName; } /*************************************************************************************/ +Int_t AliHBTAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim) +{ + //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 + UInt_t ii; for(ii = 0;ii5) + { + Info("DeleteFunctions","Deleting ParticleFunction %#x",fParticleFunctions[ii]); + Info("DeleteFunctions","Deleting ParticleFunction %s",fParticleFunctions[ii]->Name()); + } + delete fParticleFunctions[ii]; + } fNParticleFunctions = 0; for(ii = 0;ii5) + { + Info("DeleteFunctions","Deleting TrackFunction %#x",fTrackFunctions[ii]); + Info("DeleteFunctions","Deleting TrackFunction %s",fTrackFunctions[ii]->Name()); + } + delete fTrackFunctions[ii]; + } fNTrackFunctions = 0; for(ii = 0;ii5) + { + Info("DeleteFunctions","Deleting ParticleAndTrackFunction %#x",fParticleAndTrackFunctions[ii]); + Info("DeleteFunctions","Deleting ParticleAndTrackFunction %s",fParticleAndTrackFunctions[ii]->Name()); + } + delete fParticleAndTrackFunctions[ii]; + } fNParticleAndTrackFunctions = 0; for(ii = 0; ii5) + { + Info("DeleteFunctions","Deleting ParticleMonitorFunction %#x",fParticleMonitorFunctions[ii]); + Info("DeleteFunctions","Deleting ParticleMonitorFunction %s",fParticleMonitorFunctions[ii]->Name()); + } + delete fParticleMonitorFunctions[ii]; + } fNParticleMonitorFunctions = 0; for(ii = 0; ii5) + { + Info("DeleteFunctions","Deleting TrackMonitorFunction %#x",fTrackMonitorFunctions[ii]); + Info("DeleteFunctions","Deleting TrackMonitorFunction %s",fTrackMonitorFunctions[ii]->Name()); + } + delete fTrackMonitorFunctions[ii]; + } fNTrackMonitorFunctions = 0; for(ii = 0; ii5) + { + Info("DeleteFunctions","Deleting ParticleAndTrackMonitorFunction %#x",fParticleAndTrackMonitorFunctions[ii]); + Info("DeleteFunctions","Deleting ParticleAndTrackMonitorFunction %s",fParticleAndTrackMonitorFunctions[ii]->Name()); + } + delete fParticleAndTrackMonitorFunctions[ii]; + } fNParticleAndTrackMonitorFunctions = 0; + } /*************************************************************************************/ -void AliHBTAnalysis::Init() +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) + { + 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;iiInit(); @@ -165,6 +319,8 @@ void AliHBTAnalysis::Init() for(ii = 0; iiInit(); + + return 0; } /*************************************************************************************/ @@ -184,376 +340,418 @@ void AliHBTAnalysis::ResetFunctions() } } /*************************************************************************************/ - -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; - } - -} -/*************************************************************************************/ - -void AliHBTAnalysis::ProcessTracksAndParticles() +Int_t AliHBTAnalysis::ProcessRecAndSim(AliAOD* aodrec, AliAOD* aodsim) { -//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events), -//the loops are splited +//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 * trackEvent1 = 0x0,*partEvent1 = 0x0; - AliHBTEvent * trackEvent2,*partEvent2; + AliAOD * trackEvent = aodrec, *partEvent = aodsim; + AliAOD* trackEvent1 = new AliAOD(); + AliAOD* partEvent1 = new AliAOD(); + partEvent1->SetOwner(kTRUE); + trackEvent1->SetOwner(kTRUE); + + 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(); - AliHBTPair * trackpair = new AliHBTPair(); - AliHBTPair * partpair = new AliHBTPair(); + static AliHBTPair tpair; + static AliHBTPair ppair; + AliHBTPair* trackpair = &tpair; + AliHBTPair* partpair = &ppair; + AliHBTPair * tmptrackpair;//temprary pointers to pairs AliHBTPair * tmppartpair; - - AliHBTEventBuffer partbuffer(fBufferSize); - AliHBTEventBuffer trackbuffer(fBufferSize); - register UInt_t ii; - Bool_t nocorrfctns = (fNParticleFunctions == 0) && (fNTrackFunctions == 0) && (fNParticleAndTrackFunctions == 0); - fReader->Rewind(); - Int_t i = -1; - while (fReader->Next() == kFALSE) - { - i++; - partEvent= fReader->GetParticleEvent(); - trackEvent = fReader->GetTrackEvent(); - - if ( !partEvent || !trackEvent ) - { - Error("ProcessTracksAndParticles","Can not get event"); - return; - } - - if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() ) - { - Fatal("ProcessTracksAndParticles", - "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)", - i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles()); - return; - } - - if(partEvent1 == 0x0) - { - partEvent1 = new AliHBTEvent(); - partEvent1->SetOwner(kTRUE); - trackEvent1 = new AliHBTEvent(); - trackEvent1->SetOwner(kTRUE); - } - else - { - partEvent1->Reset(); - trackEvent1->Reset(); - } - - for (Int_t j = 0; jGetNumberOfParticles() ; j++) - { - /***************************************/ - /****** Looping same events ********/ - /****** filling numerators ********/ - /***************************************/ - if ( (j%fDisplayMixingInfo) == 0) - Info("ProcessTracksAndParticles", - "Mixing particle %d from event %d with particles from event %d",j,i,i); - - part1= partEvent->GetParticle(j); - track1= trackEvent->GetParticle(j); - - if( part1->GetPdgCode() != track1->GetPdgCode() ) - { - Fatal("ProcessTracksAndParticles", - "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)", - i,j, part1->GetPdgCode(),track1->GetPdgCode() ); - return; - } - - Bool_t firstcut = fPairCut->GetFirstPartCut()->Pass(part1); - - if ( (firstcut == kFALSE) || (fPairCut->GetSecondPartCut()->Pass(part1) == kFALSE) ) - { - //accepted by any cut - // we have to copy because reader keeps only one event - - partEvent1->AddParticle(new AliHBTParticle(*part1)); - trackEvent1->AddParticle(new AliHBTParticle(*track1)); + 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; jGetNumberOfParticles() ; j++) + { + /***************************************/ + /****** Looping same events ********/ + /****** filling numerators ********/ + /***************************************/ + if ( (j%fDisplayMixingInfo) == 0) + Info("ProcessTracksAndParticles", + "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; iiProcess(part1); + for(ii = 0; iiProcess(track1); + for(ii = 0; iiProcess(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;iiProcessSameEventParticles(tmppartpair); + + for(ii = 0;iiProcessSameEventParticles(tmptrackpair); + + for(ii = 0;iiProcessSameEventParticles(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(); - if (firstcut) continue; - - for(ii = 0; iiProcess(part1); - for(ii = 0; iiProcess(track1); - for(ii = 0; iiProcess(track1,part1); - - if (nocorrfctns) continue; - - for (Int_t k =j+1; k < partEvent->GetNumberOfParticles() ; k++) + m++; + if ( (j%fDisplayMixingInfo) == 0) + Info("ProcessTracksAndParticles", + "Mixing particle %d from current event with particles from event %d",j,-m); + + for(Int_t l = 0; lGetNumberOfParticles();l++) // ... on all particles { - part2= partEvent->GetParticle(k); - if (part1->GetUID() == part2->GetUID()) continue; + part2= partEvent2->GetParticle(l); partpair->SetParticles(part1,part2); - - track2= trackEvent->GetParticle(k); + + track2= trackEvent2->GetParticle(l); trackpair->SetParticles(track1,track2); - if(fPairCut->Pass(partpair) ) //check pair cut - { //do not meets crietria of the pair cut, try with swapped pairs - if( fPairCut->Pass(partpair->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; } + for(ii = 0;iiProcessSameEventParticles(tmppartpair); - + fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair); + for(ii = 0;iiProcessSameEventParticles(tmptrackpair); - + fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair); + for(ii = 0;iiProcessSameEventParticles(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; - - partbuffer.ResetIter(); - trackbuffer.ResetIter(); + fParticleAndTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair,tmppartpair); + }//for(Int_t l = 0; lGetNumberOfParticles() ; j++) + delete fPartBuffer->Push(partEvent1); + delete fTrackBuffer->Push(trackEvent1); + //end of loop over events + return 0; +} +/*************************************************************************************/ + +Int_t AliHBTAnalysis::ProcessSim(AliAOD* /*aodrec*/, AliAOD* aodsim) +{ + //Does analysis of simulated data + AliVAODParticle * part1, * part2; + + AliAOD* partEvent = aodsim; + AliAOD* partEvent1 = new AliAOD(); + partEvent1->SetOwner(kTRUE); + + AliAOD* partEvent2; + + AliHBTPair ppair; + + AliHBTPair* partpair = &ppair; + + AliHBTPair * tmppartpair; + + register UInt_t ii; + + + + if ( !partEvent ) + { + Error("ProcessRecAndSim","Can not get event"); + return 1; + } + + + for (Int_t j = 0; jGetNumberOfParticles() ; j++) + { + /***************************************/ + /****** Looping same events ********/ + /****** filling numerators ********/ + /***************************************/ + if ( (j%fDisplayMixingInfo) == 0) + Info("ProcessTracksAndParticles", + "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; iiProcess(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;iiProcessSameEventParticles(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 = partbuffer.Next() )) + while (( partEvent2 = fPartBuffer->Next() )) { - trackEvent2 = trackbuffer.Next(); - m++; if ( (j%fDisplayMixingInfo) == 0) - Info("ProcessTracksAndParticles", - "Mixing particle %d from event %d with particles from event %d",j,i,i-m); - + Info("ProcessParticles", + "Mixing particle %d from current event with particles from event %d",j,-m); for(Int_t l = 0; lGetNumberOfParticles();l++) // ... on all particles { + part2= partEvent2->GetParticle(l); partpair->SetParticles(part1,part2); - track2= trackEvent2->GetParticle(l); - trackpair->SetParticles(track1,track2); - - if( fPairCut->Pass(partpair) ) //check pair cut + if( fPairCut->Rejected(partpair) ) //check pair cut { //do not meets crietria of the - if( fPairCut->Pass(partpair->GetSwapedPair()) ) + if( fPairCut->Rejected((AliHBTPair*)partpair->GetSwappedPair()) ) continue; else { - tmppartpair = partpair->GetSwapedPair(); - tmptrackpair = trackpair->GetSwapedPair(); + tmppartpair = (AliHBTPair*)partpair->GetSwappedPair(); } } else {//meets criteria of the pair cut - tmptrackpair = trackpair; tmppartpair = partpair; } + for(ii = 0;iiProcessDiffEventParticles(tmppartpair); - for(ii = 0;iiProcessDiffEventParticles(tmptrackpair); - - for(ii = 0;iiProcessDiffEventParticles(tmptrackpair,tmppartpair); - }//for(Int_t l = 0; lGetNumberOfParticles() ; j++) - partEvent1 = partbuffer.Push(partEvent1); - trackEvent1 = trackbuffer.Push(trackEvent1); - //end of loop over events - }//while (fReader->Next() == kFALSE) -} -/*************************************************************************************/ - -void AliHBTAnalysis::ProcessTracks() -{ -//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 * trackEvent1 = 0x0; - AliHBTEvent * trackEvent2; + } + delete fPartBuffer->Push(partEvent1); + //end of loop over events + return 0; +} +/*************************************************************************************/ +Int_t AliHBTAnalysis::ProcessRec(AliAOD* aodrec, AliAOD* /*aodsim*/) +{ + //Does analysis of reconstructed data + AliVAODParticle * track1, * track2; + + AliAOD* trackEvent = aodrec; + AliAOD* trackEvent1 = new AliAOD(); + trackEvent1->SetOwner(kTRUE); - register UInt_t ii; + AliAOD* trackEvent2; + + AliHBTPair tpair; - AliHBTPair * trackpair = new AliHBTPair(); - AliHBTPair * tmptrackpair; //temporary pointer + AliHBTPair* trackpair = &tpair; + + AliHBTPair * tmptrackpair; - AliHBTEventBuffer trackbuffer(fBufferSize); + register UInt_t ii; - fReader->Rewind(); - Int_t i = -1; - while (fReader->Next() == kFALSE) - { - i++; - trackEvent = fReader->GetTrackEvent(); - if (!trackEvent) continue; - - if(trackEvent1 == 0x0) - { - trackEvent1 = new AliHBTEvent(); - trackEvent1->SetOwner(kTRUE); - } - else - { - trackEvent1->Reset(); - } - - for (Int_t j = 0; jGetNumberOfParticles() ; j++) - { - /***************************************/ - /****** Looping same events ********/ - /****** filling numerators ********/ - /***************************************/ - if ( (j%fDisplayMixingInfo) == 0) - Info("ProcessTracks", - "Mixing particle %d from event %d with particles from event %d",j,i,i); - - track1= trackEvent->GetParticle(j); - Bool_t firstcut = fPairCut->GetFirstPartCut()->Pass(track1); - - if ( (firstcut == kFALSE) || (fPairCut->GetSecondPartCut()->Pass(track1) == kFALSE) ) - { - //accepted by any cut - // we have to copy because reader keeps only one event - trackEvent1->AddParticle(new AliHBTParticle(*track1)); - } - if (firstcut) continue; + if ( !trackEvent ) + { + Error("ProcessRecAndSim","Can not get event"); + return 1; + } - for(ii = 0; iiProcess(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; + for (Int_t j = 0; jGetNumberOfParticles() ; j++) + { + /***************************************/ + /****** Looping same events ********/ + /****** filling numerators ********/ + /***************************************/ + if ( (j%fDisplayMixingInfo) == 0) + Info("ProcessTracksAndParticles", + "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; iiProcess(track1); + + if ( fNParticleFunctions == 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); - trackpair->SetParticles(track1,track2); - if(fPairCut->Pass(trackpair)) //check pair cut + 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; } - for(ii = 0;iiProcessSameEventParticles(tmptrackpair); - } - /***************************************/ - /***** Filling denominators *********/ - /***************************************/ - - if (fBufferSize == 0) continue; - - trackbuffer.ResetIter(); - - while (( trackEvent2 = trackbuffer.Next() )) + + for(ii = 0;iiProcessSameEventParticles(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() )) { - for(Int_t l = 0; lGetNumberOfParticles();l++) // ... on all particles + m++; + if ( (j%fDisplayMixingInfo) == 0) + Info("ProcessParticles", + "Mixing Particle %d from current event with Particles from event %d",j,-m); + for(Int_t l = 0; lGetNumberOfParticles();l++) // ... on all Particles { track2= trackEvent2->GetParticle(l); trackpair->SetParticles(track1,track2); - if( fPairCut->Pass(trackpair) ) //check pair cut + if( fPairCut->Rejected(trackpair) ) //check pair cut { //do not meets crietria of the - if( fPairCut->Pass(trackpair->GetSwapedPair()) ) + if( fPairCut->Rejected((AliHBTPair*)trackpair->GetSwappedPair()) ) continue; else { - tmptrackpair = trackpair->GetSwapedPair(); + tmptrackpair = (AliHBTPair*)trackpair->GetSwappedPair(); } } else @@ -564,132 +762,529 @@ void AliHBTAnalysis::ProcessTracks() for(ii = 0;iiProcessDiffEventParticles(tmptrackpair); - }//for(Int_t l = 0; lNext() == kFALSE) + 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 + 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; + + + trackEvent1 = new AliAOD(); + partEvent1 = new AliAOD(); + trackEvent1->SetOwner(kFALSE); + partEvent1->SetOwner(kFALSE);; + + /********************************/ + /* Filtering out */ + /********************************/ + if ( ( (partEvent2==0x0) || (trackEvent2==0x0)) )//in case fBufferSize == 0 and pointers are created do not eneter + { + partEvent2 = new AliAOD(); + trackEvent2 = new AliAOD(); + partEvent2->SetOwner(kTRUE); + trackEvent2->SetOwner(kTRUE); + } + + FilterOut(partEvent1, partEvent2, rawpartEvent, trackEvent1, trackEvent2, rawtrackEvent); + + for (Int_t j = 0; jGetNumberOfParticles() ; 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; iiProcess(part1); + for(ii = 0; iiProcess(track1); + for(ii = 0; iiProcess(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;iiProcessSameEventParticles(partpair); + + for(ii = 0;iiProcessSameEventParticles(trackpair); + + for(ii = 0;iiProcessSameEventParticles(trackpair,partpair); + } + } + + if ( fBufferSize == 0) continue;//do not mix diff histograms + /***************************************/ + /***** Filling denominators *********/ + /***************************************/ + fPartBuffer->ResetIter(); + fTrackBuffer->ResetIter(); + + Int_t nmonitor = 0; + + while ( (partEvent3 = fPartBuffer->Next() ) != 0x0) + { + 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;iiProcessDiffEventParticles(partpair); + + for(ii = 0;iiProcessDiffEventParticles(trackpair); + + for(ii = 0;iiProcessDiffEventParticles(trackpair,partpair); + } + }// for particles event2 + }//while event2 + }//for over particles in event1 + + delete fPartBuffer->Push(partEvent2); + delete fTrackBuffer->Push(trackEvent2); + + 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; + } + + + AliVAODParticle * part1, * part2; + + static AliAOD aodsim1; + + AliAOD* partEvent1=&aodsim1;//Particle that passes first particle cut, this event + 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 + + static AliHBTPair ppair; + + AliHBTPair* partpair = &ppair; + + register UInt_t ii; + + + partEvent1 = new AliAOD(); + partEvent1->SetOwner(kFALSE);; + + + /********************************/ + /* Filtering out */ + /********************************/ + if (partEvent2==0x0)//in case fBufferSize == 0 and pointers are created do not eneter + { + partEvent2 = new AliAOD(); + partEvent2->SetOwner(kTRUE); + } + + FilterOut(partEvent1, partEvent2, rawpartEvent); + + for (Int_t j = 0; jGetNumberOfParticles() ; 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; iiProcess(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;iiProcessSameEventParticles(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;iiProcessDiffEventParticles(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 + 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; + + AliHBTPair* trackpair = &tpair; + + register UInt_t ii; + + + trackEvent1 = new AliAOD(); + trackEvent1->SetOwner(kFALSE); + + /********************************/ + /* Filtering out */ + /********************************/ + if ( trackEvent2==0x0 )//in case fBufferSize == 0 and pointers are created do not eneter + { + trackEvent2 = new AliAOD(); + trackEvent2->SetOwner(kTRUE); + } + + FilterOut(trackEvent1, trackEvent2, rawtrackEvent); + + for (Int_t j = 0; jGetNumberOfParticles() ; 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; iiProcess(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;iiProcessSameEventParticles(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;iiProcessDiffEventParticles(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; + } + } +/*************************************************************************************/ -/*************************************************************************************/ - -void AliHBTAnalysis::ProcessParticles() +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 - AliHBTParticle * part1, * part2; - AliHBTEvent * partEvent; - AliHBTEvent * partEvent1 = 0x0; - AliHBTEvent * partEvent2; - - register UInt_t ii; - - AliHBTPair * partpair = new AliHBTPair(); - AliHBTPair * tmppartpair; //temporary pointer - AliHBTEventBuffer partbuffer(fBufferSize); +// 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; + fReader->Rewind(); - Int_t i = -1; while (fReader->Next() == kFALSE) { - i++; - partEvent = fReader->GetParticleEvent(); - if (!partEvent) continue; - - if(partEvent1 == 0x0) - { - partEvent1 = new AliHBTEvent(); - partEvent1->SetOwner(kTRUE); - } - else - { - partEvent1->Reset(); - } - - for (Int_t j = 0; jGetNumberOfParticles() ; j++) - { - /***************************************/ - /****** Looping same events ********/ - /****** filling numerators ********/ - /***************************************/ - if ( (j%fDisplayMixingInfo) == 0) - Info("ProcessParts", - "Mixing particle %d from event %d with particles from event %d",j,i,i); - - part1= partEvent->GetParticle(j); - Bool_t firstcut = fPairCut->GetFirstPartCut()->Pass(part1); - - if ( (firstcut == kFALSE) || (fPairCut->GetSecondPartCut()->Pass(part1) == kFALSE) ) - { - //accepted by any cut - // we have to copy because reader keeps only one event - partEvent1->AddParticle(new AliHBTParticle(*part1)); - } - - if (firstcut) continue; + partEvent = fReader->GetEventSim(); + trackEvent = fReader->GetEventRec(); + ProcessRecAndSim(trackEvent,partEvent); - for(ii = 0; iiProcess(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; + }//while (fReader->Next() == kFALSE) + +} +/*************************************************************************************/ - 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; - } - for(ii = 0;iiProcessSameEventParticles(tmppartpair); - } - /***************************************/ - /***** Filling denominators *********/ - /***************************************/ - - if (fBufferSize == 0) continue; - - partbuffer.ResetIter(); - - while (( partEvent2 = partbuffer.Next() )) - { - for(Int_t l = 0; lGetNumberOfParticles();l++) // ... on all particles - { +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) +} - 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 - {//meets criteria of the pair cut - tmppartpair = partpair; - } - - for(ii = 0;iiProcessDiffEventParticles(tmppartpair); - - }//for(Int_t l = 0; lRewind(); + while (fReader->Next() == kFALSE) + { + partEvent = fReader->GetEventSim(); + ProcessSim(0x0,partEvent); }//while (fReader->Next() == kFALSE) } /*************************************************************************************/ @@ -698,28 +1293,97 @@ 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;iiWrite(); + { + 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;iiWrite(); + { + 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;iiWrite(); + { + 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;iiWrite(); + { + 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;iiWrite(); + { + 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;iiWrite(); + { + 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) @@ -727,7 +1391,7 @@ void AliHBTAnalysis::SetGlobalPairCut(AliHBTPairCut* cut) Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring"); } delete fPairCut; - fPairCut = (AliHBTPairCut*)cut->Clone(); + fPairCut = (AliAODPairCut*)cut->Clone(); } /*************************************************************************************/ @@ -821,42 +1485,46 @@ 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; Info("RunCoherencyCheck","Checking HBT Runs Coherency"); - 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("AliHBTAnalysis::RunCoherencyCheck()", - "Number of simulated events (%d) is not equal to number of reconstructed events(%d)", - fReader->GetNumberOfPartEvents(),fReader->GetNumberOfTrackEvents()); - return kTRUE; - } +//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; +// } Info("RunCoherencyCheck","Checking number of Particles AND Particles Types in each event ..."); - AliHBTEvent *partEvent; - AliHBTEvent *trackEvent; + AliAOD *partEvent; + AliAOD *trackEvent; for( i = 0; iGetNumberOfTrackEvents();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; @@ -866,7 +1534,7 @@ Bool_t AliHBTAnalysis::RunCoherencyCheck() { 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; @@ -876,6 +1544,7 @@ Bool_t AliHBTAnalysis::RunCoherencyCheck() } Info("RunCoherencyCheck"," Done"); Info("RunCoherencyCheck"," Everything looks OK"); +*/ return kFALSE; } @@ -884,329 +1553,48 @@ Bool_t AliHBTAnalysis::RunCoherencyCheck() void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal() { //Performs analysis for both, tracks and particles - - AliHBTParticle * part1, * part2; - AliHBTParticle * track1, * track2; - - AliHBTEvent * trackEvent1=0x0,*partEvent1=0x0; - AliHBTEvent * trackEvent2=0x0,*partEvent2=0x0; - AliHBTEvent * trackEvent3=0x0,*partEvent3=0x0; - - AliHBTEvent * rawtrackEvent, *rawpartEvent;//this we get from Reader - - AliHBTPair * trackpair = new AliHBTPair(); - AliHBTPair * partpair = new AliHBTPair(); - - AliHBTEventBuffer partbuffer(fBufferSize); - AliHBTEventBuffer trackbuffer(fBufferSize); - - register UInt_t ii; - - trackEvent1 = new AliHBTEvent(); - partEvent1 = new AliHBTEvent(); - trackEvent1->SetOwner(kFALSE); - partEvent1->SetOwner(kFALSE);; - - Bool_t nocorrfctns = (fNParticleFunctions == 0) && (fNTrackFunctions == 0) && (fNParticleAndTrackFunctions == 0); - + AliAOD* rawtrackEvent, * rawpartEvent; fReader->Rewind(); - + Info("ProcessTracksAndParticlesNonIdentAnal","**************************************"); Info("ProcessTracksAndParticlesNonIdentAnal","***** NON IDENT MODE ****************"); Info("ProcessTracksAndParticlesNonIdentAnal","**************************************"); - + for (Int_t i = 0;;i++)//infinite loop { if (fReader->Next()) break; //end when no more events available - rawpartEvent = fReader->GetParticleEvent(); - rawtrackEvent = fReader->GetTrackEvent(); - if ( (rawpartEvent == 0x0) || (rawtrackEvent == 0x0) ) continue;//in case of any error - - if ( rawpartEvent->GetNumberOfParticles() != rawtrackEvent->GetNumberOfParticles() ) - { - Fatal("ProcessTracksAndParticlesNonIdentAnal", - "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)", - i,rawpartEvent->GetNumberOfParticles() , rawtrackEvent->GetNumberOfParticles()); - return; - } - - /********************************/ - /* Filtering out */ - /********************************/ - if ( ( (partEvent2==0x0) || (trackEvent2==0x0)) )//in case fBufferSize == 0 and pointers are created do not eneter - { - partEvent2 = new AliHBTEvent(); - trackEvent2 = new AliHBTEvent(); - partEvent2->SetOwner(kTRUE); - trackEvent2->SetOwner(kTRUE); - } - - FilterOut(partEvent1, partEvent2, rawpartEvent, trackEvent1, trackEvent2, rawtrackEvent); + rawpartEvent = fReader->GetEventSim(); + rawtrackEvent = fReader->GetEventRec(); - for (Int_t j = 0; jGetNumberOfParticles() ; j++) - { - if ( (j%fDisplayMixingInfo) == 0) - Info("ProcessTracksAndParticlesNonIdentAnal", - "Mixing particle %d from event %d with particles from event %d",j,i,i); - - part1= partEvent1->GetParticle(j); - track1= trackEvent1->GetParticle(j); - - if( part1->GetPdgCode() != track1->GetPdgCode() ) - { - Fatal("ProcessTracksAndParticlesNonIdentAnal", - "Event %d: Particle %d: PID of simulated particle (%d) not the same of reconstructed track (%d)", - i,j, part1->GetPdgCode(),track1->GetPdgCode() ); - return; - } - - for(ii = 0; iiProcess(part1); - for(ii = 0; iiProcess(track1); - for(ii = 0; iiProcess(track1,part1); - - if (nocorrfctns) 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( (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;iiProcessSameEventParticles(partpair); - - for(ii = 0;iiProcessSameEventParticles(trackpair); - - for(ii = 0;iiProcessSameEventParticles(trackpair,partpair); - } - } - - if ( fBufferSize == 0) continue;//do not mix diff histograms - /***************************************/ - /***** Filling denominators *********/ - /***************************************/ - partbuffer.ResetIter(); - trackbuffer.ResetIter(); - - Int_t nmonitor = 0; - - while ( (partEvent3 = partbuffer.Next() ) != 0x0) - { - trackEvent3 = trackbuffer.Next(); - - if ( (j%fDisplayMixingInfo) == 0) - Info("ProcessTracksAndParticlesNonIdentAnal", - "Mixing particle %d from event %d with particles from event %d",j,i,i-(++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( (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;iiProcessDiffEventParticles(partpair); - - for(ii = 0;iiProcessDiffEventParticles(trackpair); - - for(ii = 0;iiProcessDiffEventParticles(trackpair,partpair); - } - }// for particles event2 - }//while event2 - }//for over particles in event1 - - partEvent2 = partbuffer.Push(partEvent2); - trackEvent2 = trackbuffer.Push(trackEvent2); - + ProcessRecAndSimNonId(rawtrackEvent,rawpartEvent); }//end of loop over events (1) - - partbuffer.SetOwner(kTRUE); - trackbuffer.SetOwner(kTRUE); - - delete partEvent1; - delete trackEvent1; - delete partEvent2; - delete trackEvent2; - delete partpair; - delete trackpair; - } /*************************************************************************************/ void AliHBTAnalysis::ProcessTracksNonIdentAnal() { //Process Tracks only with non identical mode - AliHBTParticle * track1, * track2; - - AliHBTEvent * trackEvent1=0x0; - AliHBTEvent * trackEvent2=0x0; - AliHBTEvent * trackEvent3=0x0; - - AliHBTEvent * rawtrackEvent; - - AliHBTPair * trackpair = new AliHBTPair(); - AliHBTEventBuffer trackbuffer(fBufferSize); - - register UInt_t ii; - - trackEvent1 = new AliHBTEvent(); - trackEvent1->SetOwner(kFALSE); - + AliAOD * rawtrackEvent; fReader->Rewind(); Info("ProcessTracksNonIdentAnal","**************************************"); Info("ProcessTracksNonIdentAnal","***** NON IDENT MODE ****************"); Info("ProcessTracksNonIdentAnal","**************************************"); - for (Int_t i = 0;;i++)//infinite loop { if (fReader->Next()) break; //end when no more events available - rawtrackEvent = fReader->GetTrackEvent(); - - if (rawtrackEvent == 0x0) continue;//in case of any error - - /********************************/ - /* Filtering out */ - /********************************/ - if ( trackEvent2==0x0 )//in case fBufferSize == 0 and pointers are created do not eneter - { - trackEvent2 = new AliHBTEvent(); - trackEvent2->SetOwner(kTRUE); - } - - FilterOut(trackEvent1, trackEvent2, rawtrackEvent); - - for (Int_t j = 0; jGetNumberOfParticles() ; j++) - { - if ( (j%fDisplayMixingInfo) == 0) - Info("ProcessTracksNonIdentAnal", - "Mixing particle %d from event %d with particles from event %d",j,i,i); - - track1= trackEvent1->GetParticle(j); - - for(ii = 0; iiProcess(track1); - - if (fNTrackFunctions == 0x0) continue; - - /***************************************/ - /****** filling numerators ********/ - /****** (particles from event2) ********/ - /***************************************/ - for (Int_t k = 0; k < trackEvent2->GetNumberOfParticles() ; k++) - { - 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;iiProcessSameEventParticles(trackpair); - } - } - if ( fBufferSize == 0) continue;//do not mix diff histograms - /***************************************/ - /***** Filling denominators *********/ - /***************************************/ - trackbuffer.ResetIter(); - Int_t nmonitor = 0; - - while ( (trackEvent3 = trackbuffer.Next() ) != 0x0) - { - - if ( (j%fDisplayMixingInfo) == 0) - Info("ProcessTracksNonIdentAnal", - "Mixing particle %d from event %d with particles from event %d",j,i,i-(++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;iiProcessDiffEventParticles(trackpair); - } - }// for particles event2 - }//while event2 - }//for over particles in event1 - - trackEvent2 = trackbuffer.Push(trackEvent2); - + rawtrackEvent = fReader->GetEventRec(); + ProcessRecNonId(rawtrackEvent,0x0); }//end of loop over events (1) - - trackbuffer.SetOwner(kTRUE); - delete trackpair; - delete trackEvent1; - delete trackEvent2; } /*************************************************************************************/ void AliHBTAnalysis::ProcessParticlesNonIdentAnal() { //process paricles only with non identical mode - AliHBTParticle * part1 = 0x0, * part2 = 0x0; - - AliHBTEvent * partEvent1 = 0x0; - AliHBTEvent * partEvent2 = 0x0; - AliHBTEvent * partEvent3 = 0x0; - - AliHBTEvent * rawpartEvent = 0x0; - - AliHBTPair * partpair = new AliHBTPair(); - AliHBTEventBuffer partbuffer(fBufferSize); - - register UInt_t ii; - - partEvent1 = new AliHBTEvent(); - partEvent1->SetOwner(kFALSE); - + AliAOD * rawpartEvent = 0x0; fReader->Rewind(); Info("ProcessParticlesNonIdentAnal","**************************************"); @@ -1217,111 +1605,25 @@ void AliHBTAnalysis::ProcessParticlesNonIdentAnal() { if (fReader->Next()) break; //end when no more events available - rawpartEvent = fReader->GetParticleEvent(); - if ( rawpartEvent == 0x0 ) continue;//in case of any error - - /********************************/ - /* Filtering out */ - /********************************/ - if (partEvent2==0x0)//in case fBufferSize == 0 and pointers are created do not eneter - { - partEvent2 = new AliHBTEvent(); - partEvent2->SetOwner(kTRUE); - } - - FilterOut(partEvent1, partEvent2, rawpartEvent); - - for (Int_t j = 0; jGetNumberOfParticles() ; j++) - { - if ( (j%fDisplayMixingInfo) == 0) - Info("ProcessParticlesNonIdentAnal", - "Mixing particle %d from event %d with particles from event %d",j,i,i); - - part1= partEvent1->GetParticle(j); - - for(ii = 0; iiProcess(part1); - - if (fNParticleFunctions == 0) continue; - - /***************************************/ - /****** filling numerators ********/ - /****** (particles from event2) ********/ - /***************************************/ - for (Int_t k = 0; k < partEvent2->GetNumberOfParticles() ; k++) - { - 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;iiProcessSameEventParticles(partpair); - } - } - if ( fBufferSize == 0) continue;//do not mix diff histograms - /***************************************/ - /***** Filling denominators *********/ - /***************************************/ - partbuffer.ResetIter(); - Int_t nmonitor = 0; - - while ( (partEvent3 = partbuffer.Next() ) != 0x0) - { - if ( (j%fDisplayMixingInfo) == 0) - Info("ProcessParticlesNonIdentAnal", - "Mixing particle %d from event %d with particles from event %d",j,i,i-(++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;iiProcessDiffEventParticles(partpair); - } - } - }// for particles event2 - }//while event2 - }//for over particles in event1 - partEvent2 = partbuffer.Push(partEvent2); + rawpartEvent = fReader->GetEventSim(); + ProcessSimNonId(0x0,rawpartEvent); }//end of loop over events (1) - - partbuffer.SetOwner(kTRUE); - delete partpair; - delete partEvent1; - delete partEvent2; } /*************************************************************************************/ -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++) @@ -1330,8 +1632,8 @@ void AliHBTAnalysis::FilterOut(AliHBTEvent* outpart1, AliHBTEvent* outpart2, Ali part = inpart->GetParticle(i); track = intrack->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 ( ((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 @@ -1350,25 +1652,25 @@ void AliHBTAnalysis::FilterOut(AliHBTEvent* outpart1, AliHBTEvent* outpart2, Ali if (in2) { - outpart2->AddParticle(new AliHBTParticle(*part)); - outtrack2->AddParticle(new AliHBTParticle(*track)); + outpart2->AddParticle(part); + outtrack2->AddParticle(track); 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; @@ -1377,8 +1679,8 @@ void AliHBTAnalysis::FilterOut(AliHBTEvent* out1, AliHBTEvent* out2, AliHBTEvent 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 @@ -1428,3 +1730,34 @@ Bool_t AliHBTAnalysis::IsNonIdentAnalysis() 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(); + } +} + +/*************************************************************************************/ +