From dc2c3f36e49f1d42b10f9ccce90bf724267ab918 Mon Sep 17 00:00:00 2001 From: skowron Date: Fri, 6 Sep 2002 15:27:48 +0000 Subject: [PATCH] Faster mixing implemented for non-identical particles --- HBTAN/AliHBTAnalysis.cxx | 511 ++++++++++++++++++++++++++++++++++++++- HBTAN/AliHBTAnalysis.h | 12 +- 2 files changed, 517 insertions(+), 6 deletions(-) diff --git a/HBTAN/AliHBTAnalysis.cxx b/HBTAN/AliHBTAnalysis.cxx index 7d95ce8ec8d..ace1f2058d7 100644 --- a/HBTAN/AliHBTAnalysis.cxx +++ b/HBTAN/AliHBTAnalysis.cxx @@ -4,6 +4,7 @@ #include #include "AliHBTRun.h" +#include "AliHBTEvent.h" #include "AliHBTReader.h" #include "AliHBTParticle.h" #include "AliHBTParticleCut.h" @@ -11,6 +12,7 @@ #include "AliHBTPairCut.h" #include "AliHBTFunction.h" +#include #include @@ -89,6 +91,8 @@ void AliHBTAnalysis::Process(Option_t* option) const char *oT = strstr(option,"Tracks"); const char *oP = strstr(option,"Particles"); + Bool_t nonid = IsNonIdentAnalysis(); + if(oT && oP) { if (fReader->GetNumberOfPartEvents() <1) @@ -108,7 +112,8 @@ void AliHBTAnalysis::Process(Option_t* option) "Coherency check not passed. Maybe change the option?\n"); return; } - ProcessTracksAndParticles(); + if (nonid) ProcessTracksAndParticlesNonIdentAnal(); + else ProcessTracksAndParticles(); return; } @@ -119,7 +124,8 @@ void AliHBTAnalysis::Process(Option_t* option) Error("Process","There is no data to analyze."); return; } - ProcessTracks(); + if (nonid) ProcessTracksNonIdentAnal(); + else ProcessTracks(); return; } @@ -130,7 +136,14 @@ void AliHBTAnalysis::Process(Option_t* option) Error("Process","There is no data to analyze."); return; } - ProcessParticles(); + if (nonid) ProcessParticlesNonIdentAnal(); + else ProcessParticles(); + +// cout<<"NON ID"<Start("id"); for (Int_t i = 0;iStop("id"); +// gBenchmark->Show("id"); } @@ -701,10 +717,497 @@ Bool_t AliHBTAnalysis::RunCoherencyCheck() 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; + + trackEvent1 = new AliHBTEvent(); + partEvent1 = new AliHBTEvent(); + trackEvent1->SetOwner(kFALSE); + partEvent1->SetOwner(kFALSE);; + + cout<<"**************************************"<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); + + for (Int_t j = 0; jGetNumberOfParticles() ; j++) + { + if ( (j%100) == 0) cout<<"Mixing particle "<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;iiProcessSameEventParticles(partpair); + + for(ii = 0;iiProcessSameEventParticles(trackpair); + + for(ii = 0;iiProcessSameEventParticles(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 "<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;iiProcessDiffEventParticles(partpair); + + for(ii = 0;iiProcessDiffEventParticles(trackpair); + + for(ii = 0;iiProcessDiffEventParticles(trackpair,partpair); + } + }// for particles event2 + }//while event2 + }//for over particles in event1 + + pbuffer.AddFirst(partEvent2); + tbuffer.AddFirst(trackEvent2); + ninbuffer++; + + }//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; + + + Int_t Nev = fReader->GetNumberOfTrackEvents(); + + AliHBTPair * trackpair = new AliHBTPair(); + + TList tbuffer; + Int_t ninbuffer = 0; + + trackEvent1 = new AliHBTEvent(); + trackEvent1->SetOwner(kFALSE); + + cout<<"**************************************"<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; jGetNumberOfParticles() ; j++) + { + if ( (j%100) == 0) cout<<"Mixing particle "<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;iiProcessSameEventParticles(trackpair); + } + } + /***************************************/ + /***** Filling denominators *********/ + /***************************************/ + TIter titer(&tbuffer); + Int_t nmonitor = 0; + + while ( (trackEvent3 = (AliHBTEvent*)titer.Next()) != 0x0) + { + + if ( (j%100) == 0) cout<<"Mixing particle "<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;iiProcessDiffEventParticles(trackpair); + + } + }// for particles event2 + }//while event2 + }//for over particles in event1 + + tbuffer.AddFirst(trackEvent2); + ninbuffer++; + + }//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);; + + cout<<"**************************************"<Start("non_id"); + for (Int_t i = 0;iGetParticleEvent(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); + + for (Int_t j = 0; jGetNumberOfParticles() ; j++) + { + if ( (j%100) == 0) cout<<"Mixing particle "<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;iiProcessSameEventParticles(partpair); + + } + } + /***************************************/ + /***** Filling denominators *********/ + /***************************************/ + TIter piter(&pbuffer); + Int_t nmonitor = 0; + + while ( (partEvent3 = (AliHBTEvent*)piter.Next()) != 0x0) + { + if ( (j%100) == 0) cout<<"Mixing particle "<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;iiProcessDiffEventParticles(partpair); + + } + }// for particles event2 + }//while event2 + }//for over particles in event1 + + pbuffer.AddFirst(partEvent2); + ninbuffer++; + + }//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) +{ + //Puts particles accepted as a first particle by global cut in out1 + //and as a second particle in out2 + + AliHBTParticle* 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++) + { + in1 = in2 = kTRUE; + 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 (gDebug)//to be removed in real analysis + if ( in1 && in2 ) //both cuts accepted, should never happen, just in case + { + //Particle accpted by both cuts + Error("FilterOut","Particle accepted by both cuts"); + continue; + } + + if (in1) + { + outpart1->AddParticle(part); + outtrack1->AddParticle(track); + continue; + } + + if (in2) + { + outpart2->AddParticle(part); + outtrack2->AddParticle(track); + continue; + } + } + +} + /*************************************************************************************/ +void AliHBTAnalysis::FilterOut(AliHBTEvent* out1, AliHBTEvent* out2, AliHBTEvent* in) +{ + //Puts particles accepted as a first particle by global cut in out1 + //and as a second particle in out2 + AliHBTParticle* part; + + out1->Reset(); + out2->Reset(); + + AliHBTParticleCut *cut1 = fPairCut->GetFirstPartCut(); + AliHBTParticleCut *cut2 = fPairCut->GetSecondPartCut(); + + Bool_t in1, in2; + + for (Int_t i = 0; i < in->GetNumberOfParticles(); i++) + { + 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 (gDebug)//to be removed in real analysis + if ( in1 && in2 ) //both cuts accepted, should never happen, just in case + { + //Particle accpted by both cuts + Error("FilterOut","Particle accepted by both cuts"); + continue; + } + if (in1) + { + out1->AddParticle(part); + continue; + } + + if (in2) + { + out2->AddParticle(part); + continue; + } + } +} +/*************************************************************************************/ +Bool_t AliHBTAnalysis::IsNonIdentAnalysis() +{ + //checks if it is possible to use special analysis for non identical particles + //it means - in global pair cut first particle id is different than second one + //and both are different from 0 + //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; + + Int_t id1 = fPairCut->GetFirstPartCut()->GetPID(); + Int_t id2 = fPairCut->GetSecondPartCut()->GetPID(); + + if ( (id1==0) || (id2==0) || (id1==id2) ) return kFALSE; + + return kTRUE; +} diff --git a/HBTAN/AliHBTAnalysis.h b/HBTAN/AliHBTAnalysis.h index 45e7e31d896..8dbddb57f2c 100644 --- a/HBTAN/AliHBTAnalysis.h +++ b/HBTAN/AliHBTAnalysis.h @@ -11,6 +11,7 @@ class AliHBTPairCut; class AliHBTPair; class AliHBTRun; +class AliHBTEvent; class AliHBTReader; class AliHBTOnePairFctn; class AliHBTTwoPairFctn; @@ -41,11 +42,15 @@ class AliHBTAnalysis: public TObject void WriteFunctions(); void SetBufferSize(Int_t buffsize){fBufferSize=buffsize;} - + + Bool_t IsNonIdentAnalysis(); protected: Bool_t RunCoherencyCheck(); + void FilterOut(AliHBTEvent* outpart1, AliHBTEvent* outpart2, AliHBTEvent* inpart, + AliHBTEvent* outtrack1, AliHBTEvent* outtrack2, AliHBTEvent* intrack); + void FilterOut(AliHBTEvent* out1, AliHBTEvent* out2, AliHBTEvent* in); AliHBTReader* fReader;//! @@ -53,10 +58,13 @@ class AliHBTAnalysis: public TObject virtual void ProcessParticles(); virtual void ProcessTracksAndParticles(); + virtual void ProcessTracksAndParticlesNonIdentAnal(); + virtual void ProcessParticlesNonIdentAnal(); + virtual void ProcessTracksNonIdentAnal(); AliHBTOnePairFctn** fTrackFunctions; //!array of pointers to functions that analyze rekonstructed tracks AliHBTOnePairFctn** fParticleFunctions; //!array of pointers to functions that analyze generated particles - AliHBTTwoPairFctn** fParticleAndTrackFunctions; //!array of pointers to functions that analyze both + AliHBTTwoPairFctn** fParticleAndTrackFunctions; //!array of pointers to functions that analyze both //reconstructed tracks and generated particles //i.e. - resolution analyzers UInt_t fNTrackFunctions; //! -- 2.43.0