which included commits to RCS files with non-trunk default branches.
--- /dev/null
+
+Unix.*.Root.MacroPath: .:$(ALICE_ROOT)/macros:$(ALICE_ROOT)
+
+Root.Html.OutputDir: $(ALICE_ROOT)/html/roothtml/HBTANALYSIS
+Root.Html.SourceDir: ./
+Root.Html.Author: Piotr Krzysztof Skowronski
+Root.Html.Root: http://root.cern.ch/root/html
--- /dev/null
+
+#include "AliHBTAnalysis.h"
+
+#include <iostream.h>
+
+#include "AliHBTRun.h"
+#include "AliHBTReader.h"
+#include "AliHBTParticle.h"
+#include "AliHBTParticleCut.h"
+#include "AliHBTPair.h"
+#include "AliHBTPairCut.h"
+#include "AliHBTFunction.h"
+
+#include <TList.h>
+
+
+
+ClassImp(AliHBTAnalysis)
+
+const UInt_t AliHBTAnalysis::fgkFctnArraySize = 100;
+const Int_t AliHBTAnalysis::fgkHbtAnalyzeAll = 0;
+
+AliHBTAnalysis::AliHBTAnalysis()
+ {
+ fReader = 0x0;
+
+ fTrackFunctions = new AliHBTTwoPartFctn* [fgkFctnArraySize];
+ fParticleFunctions = new AliHBTTwoPartFctn* [fgkFctnArraySize];
+ fParticleAndTrackFunctions = new AliHBTFourPartFctn* [fgkFctnArraySize];
+
+ fNTrackFunctions = 0;
+ fNParticleFunctions = 0;
+ fNParticleAndTrackFunctions = 0;
+
+ fPairCut = new AliHBTEmptyPairCut();//empty cut - accepts all particles
+
+ }
+
+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
+ delete [] fTrackFunctions;
+ delete [] fParticleFunctions;
+ delete [] fParticleAndTrackFunctions;
+
+ delete fPairCut; // always have an copy of an object - we create we dstroy
+ }
+
+/*************************************************************************************/
+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("AliHBTAnalysis::Process","The reader is not set");
+ return;
+ }
+
+
+ const char *oT = strstr(option,"Tracks");
+ const char *oP = strstr(option,"Particles");
+
+ if(oT && oP)
+ {
+ if (fReader->GetNumberOfPartEvents() <1)
+ {
+ Error("AliHBTAnalysis::Process","There is no Particles. Maybe change the option?");
+ return;
+ }
+ if (fReader->GetNumberOfTrackEvents() <1)
+ {
+ Error("AliHBTAnalysis::Process","There is no Tracks. Maybe change the option?");
+ return;
+ }
+
+ if ( RunCoherencyCheck() )
+ {
+ Error("AliHBTAnalysis::Process",
+ "Coherency check not passed. Maybe change the option?\n");
+ return;
+ }
+ ProcessTracksAndParticles();
+ return;
+ }
+
+ if(oT)
+ {
+ ProcessTracks();
+ return;
+ }
+
+ if(oP)
+ {
+ ProcessParticles();
+ return;
+ }
+
+}
+
+/*************************************************************************************/
+void AliHBTAnalysis::ProcessTracksAndParticles()
+{
+
+//In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
+//the loops are splited
+
+
+ AliHBTParticle * part1, * part2;
+ AliHBTParticle * track1, * track2;
+
+ AliHBTEvent * trackEvent, *partEvent;
+ AliHBTEvent * 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();
+
+ 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%10) == 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++)
+ {
+ part2= partEvent->GetParticle(k);
+ partpair->SetParticles(part1,part2);
+
+ track2= trackEvent->GetParticle(k);
+ 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
+ else
+ { //swaped pair meets all the criteria
+ partpair = partpair->GetSwapedPair();
+ trackpair = trackpair->GetSwapedPair();
+ }
+ }
+ 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 *********/
+ /***************************************/
+ for (Int_t i = 0;i<Nev;i++) //In each event ....
+ {
+
+ 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 ...
+ {
+// if (N>MAXCOMB) break;
+
+ part1= partEvent->GetParticle(j);
+
+ track1= trackEvent->GetParticle(j);
+
+// for (Int_t k = i+1; k<Nev;k++) // ... Loop over all proceeding events ...
+ Int_t NNN;
+
+ if ( (i+2) < Nev) NNN = i+2;
+ else NNN = Nev;
+
+ for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
+ {
+
+ partEvent2= fReader->GetParticleEvent(k);
+ if (!partEvent2) continue;
+
+ trackEvent2 = fReader->GetTrackEvent(k);
+
+ if ( (j%10) == 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
+ {
+ partpair = partpair->GetSwapedPair();
+ trackpair = trackpair->GetSwapedPair();
+ }
+ }
+ 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(Int_t l = 0; l<N2;l++) // ... on all particles
+ }//for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
+ }
+
+ }
+
+ /***************************************/
+
+
+}
+/*************************************************************************************/
+
+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 * trackEvent2;
+
+// 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();
+
+ 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%10) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
+
+ 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
+ { //do not meets crietria of the
+ if( fPairCut->Pass(trackpair->GetSwapedPair()) ) continue;
+ else trackpair = trackpair->GetSwapedPair();
+ }
+
+ UInt_t ii;
+
+ for(ii = 0;ii<fNTrackFunctions;ii++)
+ fTrackFunctions[ii]->ProcessSameEventParticles(trackpair);
+
+
+ }
+ }
+ }
+
+ /***************************************/
+ /***** Filling diff histogram *********/
+ /***************************************/
+ for (Int_t i = 0;i<Nev;i++) //In each event ....
+ {
+ 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);
+
+// for (Int_t k = i+1; k<Nev;k++) // ... Loop over all proceeding events ...
+ Int_t NNN;
+
+ if ( (i+2) < Nev) NNN = i+2;
+ else NNN = Nev;
+
+ for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
+ {
+
+ trackEvent2 = fReader->GetTrackEvent(k);
+ if (!trackEvent2) continue;
+
+ if ( (j%10) == 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 trackpair = trackpair->GetSwapedPair();
+ }
+ UInt_t ii;
+ for(ii = 0;ii<fNTrackFunctions;ii++)
+ fTrackFunctions[ii]->ProcessDiffEventParticles(trackpair);
+
+ }//for(Int_t l = 0; l<N2;l++) // ... on all particles
+ }//for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
+ }
+
+ }
+
+ /***************************************/
+
+
+}
+
+/*************************************************************************************/
+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;
+
+ AliHBTEvent * partEvent;
+ AliHBTEvent * partEvent2;
+
+// Int_t N1, N2, N=0; //number of particles in current event(we prcess two events in one time)
+
+ Int_t Nev = fReader->GetNumberOfPartEvents();
+
+ /***************************************/
+ /****** Looping same events ********/
+ /****** filling numerators ********/
+ /***************************************/
+ AliHBTPair * partpair = new AliHBTPair();
+
+ 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%10) == 0) cout<<"Mixing particle "<<j<<" from event "<<i<<" with particles from event "<<i<<endl;
+
+ 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
+ partpair = partpair->GetSwapedPair();
+ }
+ }
+
+ UInt_t ii;
+
+ for(ii = 0;ii<fNParticleFunctions;ii++)
+ fParticleFunctions[ii]->ProcessSameEventParticles(partpair);
+ }
+ }
+ }
+
+ /***************************************/
+ /***** Filling diff histogram *********/
+ /***************************************/
+ for (Int_t i = 0;i<Nev;i++) //In each event ....
+ {
+ partEvent= fReader->GetParticleEvent(i);
+ if (!partEvent) continue;
+
+// N=0;
+
+ for (Int_t j = 0; j< partEvent->GetNumberOfParticles(); j++) // ... Loop over all particles ...
+ {
+// if (N>MAXCOMB) break;
+
+ part1= partEvent->GetParticle(j);
+
+// for (Int_t k = i+1; k<Nev;k++) // ... Loop over all proceeding events ...
+ Int_t NNN;
+
+ if ( (i+2) < Nev) NNN = i+2; //loop over next event
+ else NNN = Nev;
+
+ for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
+ {
+
+ partEvent2= fReader->GetParticleEvent(k);
+ if (!partEvent2) continue;
+
+ if ( (j%10) == 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 partpair = partpair->GetSwapedPair();
+ }
+ UInt_t ii;
+ for(ii = 0;ii<fNParticleFunctions;ii++)
+ fParticleFunctions[ii]->ProcessDiffEventParticles(partpair);
+
+ }//for(Int_t l = 0; l<N2;l++) // ... on all particles
+ }//for (Int_t k = i+1; k<NNN;k++) // ... Loop over next event
+ }
+
+ }
+
+ /***************************************/
+
+
+}
+
+/*************************************************************************************/
+
+
+void AliHBTAnalysis::Write()
+{
+ UInt_t ii;
+ for(ii = 0;ii<fNParticleFunctions;ii++)
+ fParticleFunctions[ii]->Write();
+
+ for(ii = 0;ii<fNTrackFunctions;ii++)
+ fTrackFunctions[ii]->Write();
+
+ for(ii = 0;ii<fNParticleAndTrackFunctions;ii++)
+ fParticleAndTrackFunctions[ii]->Write();
+}
+/*************************************************************************************/
+void AliHBTAnalysis::SetGlobalPairCut(AliHBTPairCut* cut)
+{
+ if (cut == 0x0)
+ {
+ Error("AliHBTAnalysis::SetGlobalPairCut","Pointer is NULL. Ignoring");
+ }
+ delete fPairCut;
+ fPairCut = (AliHBTPairCut*)cut->Clone();
+}
+
+/*************************************************************************************/
+void AliHBTAnalysis::AddTrackFunction(AliHBTTwoPartFctn* f)
+{
+ if (f == 0x0) return;
+ if (fNTrackFunctions == fgkFctnArraySize)
+ {
+ Error("AliHBTAnalysis::AddTrackFunction","Can not add this function, not enough place in the array.");
+ }
+ fTrackFunctions[fNTrackFunctions] = f;
+ fNTrackFunctions++;
+}
+void AliHBTAnalysis::AddParticleFunction(AliHBTTwoPartFctn* f)
+{
+ if (f == 0x0) return;
+
+ if (fNParticleFunctions == fgkFctnArraySize)
+ {
+ Error("AliHBTAnalysis::AddParticleFunction","Can not add this function, not enough place in the array.");
+ }
+ fParticleFunctions[fNParticleFunctions] = f;
+ fNParticleFunctions++;
+
+
+}
+void AliHBTAnalysis::AddParticleAndTrackFunction(AliHBTFourPartFctn* f)
+{
+ if (f == 0x0) return;
+ if (fNParticleAndTrackFunctions == fgkFctnArraySize)
+ {
+ Error("AliHBTAnalysis::AddParticleAndTrackFunction","Can not add this function, not enough place in the array.");
+ }
+ fParticleAndTrackFunctions[fNParticleAndTrackFunctions] = f;
+ fNParticleAndTrackFunctions++;
+}
+
+
+/*************************************************************************************/
+
+
+/*************************************************************************************/
+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;
+
+ cout<<"Number of events ..."; fflush(stdout);
+
+ 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;
+ for( i = 0; i<fReader->GetNumberOfTrackEvents();i++)
+ {
+ partEvent= fReader->GetParticleEvent(i); //gets the "ith" event
+ trackEvent = fReader->GetTrackEvent(i);
+
+ if ( (partEvent == 0x0) && (partEvent == 0x0) ) continue;
+ if ( (partEvent == 0x0) || (partEvent == 0x0) )
+ {
+ Error("AliHBTAnalysis::RunCoherencyCheck()",
+ "One event is NULL and the other one not. Event Number %d",i);
+ return kTRUE;
+ }
+
+ if ( partEvent->GetNumberOfParticles() != trackEvent->GetNumberOfParticles() )
+ {
+ Error("AliHBTAnalysis::RunCoherencyCheck()",
+ "Event %d: Number of simulated particles (%d) not equal to number of reconstructed tracks (%d)",
+ i,partEvent->GetNumberOfParticles() , trackEvent->GetNumberOfParticles());
+ return kTRUE;
+ }
+ else
+ for (Int_t j = 0; j<partEvent->GetNumberOfParticles(); j++)
+ {
+ if( partEvent->GetParticle(j)->GetPdgCode() != trackEvent->GetParticle(j)->GetPdgCode() )
+ {
+ Error("AliHBTAnalysis::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;
+ return kFALSE;
+}
+
+
+/*************************************************************************************/
+
+
+/*************************************************************************************/
+
+
--- /dev/null
+#ifndef ALIHBTANALYSIS_H
+#define ALIHBTANALYSIS_H
+
+#include <TObject.h>
+
+
+
+class AliHBTParticleCut;
+class AliHBTCut;
+class AliHBTPairCut;
+class AliHBTPair;
+
+class AliHBTRun;
+class AliHBTReader;
+class AliHBTTwoPartFctn;
+class AliHBTFourPartFctn;
+
+
+class TList;
+
+class AliHBTAnalysis: public TObject
+ {
+ public:
+ AliHBTAnalysis();
+
+ virtual ~AliHBTAnalysis();
+
+ virtual void Process(Option_t* option = "TracksAndParticles");
+
+
+ void SetGlobalPairCut(AliHBTPairCut* cut);
+
+ void AddTrackFunction(AliHBTTwoPartFctn*);
+ void AddParticleFunction(AliHBTTwoPartFctn*);
+ void AddParticleAndTrackFunction(AliHBTFourPartFctn*);
+
+ void AddResolutionFunction(AliHBTFourPartFctn* f){AddParticleAndTrackFunction(f);}
+
+ void SetReader(AliHBTReader* r){fReader = r;}
+
+ void Write();
+ protected:
+
+ Bool_t RunCoherencyCheck();
+
+
+ AliHBTReader* fReader;
+
+ virtual void ProcessTracks();
+ virtual void ProcessParticles();
+ virtual void ProcessTracksAndParticles();
+
+
+ AliHBTTwoPartFctn** fTrackFunctions; //array of pointers to functions that analyze rekonstructed tracks
+ AliHBTTwoPartFctn** fParticleFunctions; //array of pointers to functions that analyze generated particles
+ AliHBTFourPartFctn** fParticleAndTrackFunctions; //array of pointers to functions that analyze both
+ //reconstructed tracks and generated particles
+ //i.e. - resolution analyzers
+ UInt_t fNTrackFunctions;
+ UInt_t fNParticleFunctions;
+ UInt_t fNParticleAndTrackFunctions;
+
+ /**********************************************/
+ /* Control parameters */
+
+ AliHBTPairCut *fPairCut;
+
+ // AliHBTCut *fParticleCut;
+ /**********************************************/
+
+
+ private:
+ static const Int_t fgkHbtAnalyzeAll;
+ static const UInt_t fgkFctnArraySize;
+/*********************************************/
+ public:
+ ClassDef(AliHBTAnalysis,0)
+ };
+
+
+
+
+#endif
--- /dev/null
+#include "AliHBTCorrelFctn.h"
+
+
+
+ClassImp(AliHBTQInvCorrelFctn)
+
+//Corroleation function is created from dividing two histograms of QInvariant:
+// of particles from the same evnt
+//by
+// of particles from different events
+
+TH1* AliHBTQInvCorrelFctn::GetResult()
+{
+ return GetRatio(GetDenominator()->GetMaximum()/GetNumerator()->GetMaximum());
+}
+
+
+ClassImp(AliHBTInvMassCorrelFctn)
+
+AliHBTInvMassCorrelFctn::
+AliHBTInvMassCorrelFctn(Int_t nbins, Double_t maxXval, Double_t minXval):
+ AliHBTTwoPartFctn1D(nbins,maxXval,minXval)
+{
+ Rename("InvMass CF","Invariant Mass Correlation Function");
+}
+
+TH1* AliHBTInvMassCorrelFctn::GetResult()
+{
+ return GetRatio(GetDenominator()->GetMaximum()/GetNumerator()->GetMaximum());
+}
--- /dev/null
+#ifndef ALIHBTCORRELFUNCTION_H
+#define ALIHBTCORRELFUNCTION_H
+
+#include "AliHBTFunction.h"
+#include "AliHBTParticle.h"
+//Set of functions:
+// Q Invaraint Correlation Function
+// Invariant Mass Function
+//
+//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
+//Piotr.Skowronski@cern.ch
+
+
+class AliHBTQInvCorrelFctn: public AliHBTTwoPartFctn1D
+{
+//Q Invaraint Correlation Function
+//1D two particle function
+ public:
+ AliHBTQInvCorrelFctn(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0):
+ AliHBTTwoPartFctn1D(nbins,maxXval,minXval){}
+ virtual ~AliHBTQInvCorrelFctn(){};
+ TH1* GetResult();
+ protected:
+ Double_t GetValue(AliHBTPair * pair){return pair->GetQInv();}
+ public:
+ ClassDef(AliHBTQInvCorrelFctn,1)
+
+};
+
+
+class AliHBTInvMassCorrelFctn: public AliHBTTwoPartFctn1D
+{
+// Invariant Mass Function
+ public:
+ AliHBTInvMassCorrelFctn(Int_t nbins = 2000, Double_t maxXval = 2., Double_t minXval = 0.0);
+ virtual ~AliHBTInvMassCorrelFctn(){};
+ TH1* GetResult();
+ protected:
+ Double_t GetValue(AliHBTPair * pair) { return pair->GetInvMass();}
+ public:
+ ClassDef(AliHBTInvMassCorrelFctn,1)
+
+};
+
+
+
+#endif
--- /dev/null
+#include "AliHBTEvent.h"
+#include "AliHBTParticle.h"
+
+ClassImp(AliHBTEvent)
+
+const UInt_t AliHBTEvent::fgkInitEventSize = 10000;
+
+
+
+/**************************************************************************/
+
+AliHBTEvent::AliHBTEvent()
+ {
+ if(fgkInitEventSize<1)
+ {
+ Fatal("AliHBTEvent::AliHBTEvent()",
+ "fgkInitEventSize has a stiupid value (%d). Change it to positive number and recompile",
+ fgkInitEventSize);
+
+ }
+ fSize=fgkInitEventSize;
+ fParticles = new AliHBTParticle* [fSize];
+ fNParticles = 0;
+ }
+/**************************************************************************/
+
+AliHBTEvent::~AliHBTEvent()
+ {
+ this->Reset();//delete all particles
+ if(fParticles)
+ {
+ delete [] fParticles; //and delete array itself
+ }
+ fParticles = 0x0;
+ }
+/**************************************************************************/
+void AliHBTEvent::Reset()
+{
+ //deletes all particles from the event
+ if(fParticles)
+ {
+ for(Int_t i =0; i<fNParticles; i++)
+ delete fParticles[i];
+ }
+ fNParticles = 0;
+}
+
+AliHBTParticle* AliHBTEvent::GetParticleSafely(Int_t n)
+ {
+ if( (n<0) || (fNParticles<=n) ) return 0x0;
+ else return fParticles[n];
+
+ }
+/**************************************************************************/
+
+void AliHBTEvent:: AddParticle(AliHBTParticle* hbtpart)
+ {
+ //Adds new perticle to the event
+ if ( fNParticles+1 >= fSize) Expand(); //if there is no space in array, expand it
+ fParticles[fNParticles++] = hbtpart; //add a pointer
+ }
+
+/**************************************************************************/
+void AliHBTEvent::AddParticle(TParticle* part)
+ {
+ AddParticle( new AliHBTParticle(*part) );
+ }
+/**************************************************************************/
+void AliHBTEvent::
+AddParticle(Int_t pdg, Double_t px, Double_t py, Double_t pz, Double_t etot,
+ Double_t vx, Double_t vy, Double_t vz, Double_t time)
+ {
+ AddParticle(new AliHBTParticle(pdg,px,py,pz,etot,vx,vy,vz,time) );
+ }
+/**************************************************************************/
+
+void AliHBTEvent::Expand()
+ {
+ //expands the array with pointers to particles
+ //about the size defined in fgkInitEventSize
+
+ fSize+=fgkInitEventSize;
+ AliHBTParticle** tmpParticles = new AliHBTParticle* [fSize]; //create new array of pointers
+ //check if we got memory / if not Abort
+ if (!tmpParticles) Fatal("AliHBTEvent::Expand()","No more space in memory");
+
+
+ for(Int_t i = 0; i<fNParticles ;i++)
+ //copy pointers to the new array
+ {
+ tmpParticles[i] = fParticles[i];
+ }
+ delete [] fParticles; //delete old array
+ fParticles = tmpParticles; //copy new pointer to the array of pointers to particles
+
+ }
+
+
+
--- /dev/null
+#ifndef ALIHBTEvent_H
+#define ALIHBTEvent_H
+//This class sters HBT perticles for one event
+//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
+
+#include <TObject.h>
+
+class AliHBTParticle;
+class TParticle;
+class AliHBTEvent: public TObject
+ {
+ public:
+ AliHBTEvent();
+ virtual ~AliHBTEvent();
+ const static UInt_t fgkInitEventSize; //initial number of the array
+ //if expanded, this size is used also
+ AliHBTParticle* GetParticle(Int_t n); //gets particle
+ AliHBTParticle* GetParticleSafely(Int_t n); //gets particle with index check
+
+ void AddParticle(AliHBTParticle*); //adds particle to the event
+ void AddParticle(TParticle*); //adds particle to the event
+ void AddParticle(Int_t pdg, Double_t px, Double_t py, Double_t pz, Double_t etot,
+ Double_t vx, Double_t vy, Double_t vz, Double_t time);
+
+ Int_t GetNumberOfParticles() const;
+ void Reset(); //deletes all entries
+ protected:
+ AliHBTParticle ** fParticles; //!array of pointers to the particles
+ Int_t fNParticles; //!number of particles in Event
+ Int_t fSize; //!current size of the array
+
+ void Expand(); //expands the array if necessary
+ private:
+
+ public:
+ ClassDef(AliHBTEvent,1)
+ };
+/**************************************************************************/
+
+inline
+AliHBTParticle* AliHBTEvent::GetParticle(Int_t n)
+ {
+ //Gets particle without boundary check
+ return fParticles[n];
+ }
+
+/**************************************************************************/
+
+inline
+Int_t AliHBTEvent::GetNumberOfParticles() const
+ {
+ //reurns number of particles in this event
+ return fNParticles;
+ }
+
+#endif
--- /dev/null
+#include "AliHBTFunction.h"
+/******************************************************************/
+/*
+Piotr Krzysztof Skowronski
+Piotr.Skowronski@cern.ch
+Base classes for HBT functions
+
+ function
+ / \
+ / \
+ / \
+ / \
+ / \
+ / \
+ / \
+ two part four part
+ / | \ / | \
+ / | \ / | \
+ 1D 2D 3D 1D 2D 3D
+
+ four particle functions are intendent to be resolution functions:
+ it is mecessary to have simulated particle pair corresponding to given
+ recontructed track pair in order to calculate function simualted value
+ and recontructed value to be further histogrammed
+
+*/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTFunction )
+
+AliHBTFunction::AliHBTFunction()
+{
+
+ fPairCut = new AliHBTEmptyPairCut(); //dummy cut
+}
+
+void AliHBTFunction::
+Write()
+ {
+ if (GetNumerator()) GetNumerator()->Write();
+ if (GetDenominator()) GetDenominator()->Write();
+ TH1* res = GetResult();
+ if (res) GetResult()->Write();
+ }
+/******************************************************************/
+
+TH1* AliHBTFunction::
+GetRatio(Double_t normfactor)
+ {
+ TString str = fName + " ratio";
+ TH1 *result = (TH1*)GetNumerator()->Clone(str.Data());
+
+ result->SetTitle(str.Data());
+ result->Sumw2();
+
+ result->Divide(GetNumerator(),GetDenominator(),normfactor);
+
+ return result;
+
+ }
+/******************************************************************/
+void AliHBTFunction::SetPairCut(AliHBTPairCut* cut)
+{
+//Sets new Pair Cut. Old one is deleted
+//Note that it is created new object instead of simple pointer set
+//I do not want to have pointer
+//to object created somewhere else
+//because in that case I could not believe that
+//it would always exist (sb could delete it)
+//so we have always own copy
+
+ if(!cut)
+ {
+ Error("AliHBTFunction::SetPairCut","argument is NULL");
+ return;
+ }
+ delete fPairCut;
+ fPairCut = (AliHBTPairCut*)cut->Clone();
+
+}
+
+/******************************************************************/
+
+void AliHBTFunction::
+Rename(const Char_t * name)
+ {
+ //renames the function and histograms
+ SetName(name);
+ SetTitle(name);
+
+ TString numstr = fName + " Numerator"; //title and name of the
+ //numerator histogram
+ TString denstr = fName + " Denominator";//title and name of the
+ //denominator histogram
+
+ GetNumerator()->SetName(numstr.Data());
+ GetNumerator()->SetTitle(numstr.Data());
+
+ GetDenominator()->SetName(denstr.Data());
+ GetDenominator()->SetTitle(denstr.Data());
+
+ }
+
+void AliHBTFunction::
+Rename(const Char_t * name, const Char_t * title)
+ {
+ //renames and retitle the function and histograms
+
+ SetName(name);
+ SetTitle(title);
+
+ TString numstrn = fName + " Numerator"; //name of the
+ //numerator histogram
+
+ TString numstrt = fTitle + " Numerator"; //title of the
+ //numerator histogram
+
+ TString denstrn = fName + " Denominator";//name of the
+ //denominator histogram
+
+ TString denstrt = fTitle + " Denominator";//title of the
+ //denominator histogram
+
+
+ GetNumerator()->SetName(numstrn.Data());
+ GetNumerator()->SetTitle(numstrt.Data());
+
+ GetDenominator()->SetName(denstrn.Data());
+ GetDenominator()->SetTitle(denstrt.Data());
+
+
+ }
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTTwoPartFctn )
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTFourPartFctn)
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTTwoPartFctn1D )
+
+AliHBTTwoPartFctn1D::
+AliHBTTwoPartFctn1D(Int_t nbins, Double_t maxXval, Double_t minXval)
+ {
+ //Constructor of Two Part One Dimentional Function
+ // nbins: number of bins in histograms - default 100
+ // maxXval and minXval: range of histgram(s) default 0 - 0.15 (GeV)
+
+
+ TString numstr = fName + " Numerator"; //title and name of the
+ //numerator histogram
+ TString denstr = fName + " Denominator";//title and name of the
+ //denominator histogram
+
+ fNumerator = new TH1D(numstr.Data(),numstr.Data(),nbins,minXval,maxXval);
+ fDenominator = new TH1D(denstr.Data(),denstr.Data(),nbins,minXval,maxXval);
+
+ fNumerator->Sumw2();
+ fDenominator->Sumw2();
+
+ }
+/******************************************************************/
+AliHBTTwoPartFctn1D::~AliHBTTwoPartFctn1D()
+{
+ delete fNumerator;
+ delete fDenominator;
+}
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTTwoPartFctn2D )
+
+AliHBTTwoPartFctn2D::
+AliHBTTwoPartFctn2D(Int_t nXbins, Double_t maxXval, Double_t minXval ,
+ Int_t nYbins, Double_t maxYval, Double_t minYval)
+
+{
+ TString numstr = fName + " Numerator"; //title and name of the
+ //numerator histogram
+ TString denstr = fName + " Denominator";//title and name of the
+ //denominator histogram
+
+ fNumerator = new TH2D(numstr.Data(),numstr.Data(),
+ nXbins,minXval,maxXval,
+ nYbins,minYval,maxYval);
+
+ fDenominator = new TH2D(denstr.Data(),denstr.Data(),
+ nXbins,minXval,maxXval,
+ nYbins,minYval,maxYval);
+
+ fNumerator->Sumw2();
+ fDenominator->Sumw2();
+
+}
+AliHBTTwoPartFctn2D::~AliHBTTwoPartFctn2D()
+{
+ delete fNumerator;
+ delete fDenominator;
+}
+void AliHBTTwoPartFctn2D::ProcessSameEventParticles(AliHBTPair* pair)
+{
+ pair = CheckPair(pair);
+ if(pair)
+ {
+ Double_t x,y;
+ GetValues(pair,x,y);
+ fNumerator->Fill(y,x);
+ }
+}
+
+void AliHBTTwoPartFctn2D::ProcessDiffEventParticles(AliHBTPair* pair)
+{
+ pair = CheckPair(pair);
+ if(pair)
+ {
+ Double_t x,y;
+ GetValues(pair,x,y);
+ fDenominator->Fill(y,x);
+ }
+
+}
+
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTTwoPartFctn3D)
+
+AliHBTTwoPartFctn3D::
+AliHBTTwoPartFctn3D(Int_t nXbins, Double_t maxXval, Double_t minXval,
+ Int_t nYbins, Double_t maxYval, Double_t minYval,
+ Int_t nZbins, Double_t maxZval, Double_t minZval)
+
+{
+ TString numstr = fName + " Numerator"; //title and name of the
+ //numerator histogram
+ TString denstr = fName + " Denominator";//title and name of the
+ //denominator histogram
+
+ fNumerator = new TH3D(numstr.Data(),numstr.Data(),
+ nXbins,minXval,maxXval,
+ nYbins,minYval,maxYval,
+ nZbins,minZval,maxZval);
+
+ fDenominator = new TH3D(denstr.Data(),denstr.Data(),
+ nXbins,minXval,maxXval,
+ nYbins,minYval,maxYval,
+ nZbins,minZval,maxZval);
+
+ fNumerator->Sumw2();
+ fDenominator->Sumw2();
+
+}
+
+
+AliHBTTwoPartFctn3D::~AliHBTTwoPartFctn3D()
+{
+ delete fNumerator;
+ delete fDenominator;
+}
+
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+ClassImp( AliHBTFourPartFctn2D)
+
+
+AliHBTFourPartFctn2D::
+AliHBTFourPartFctn2D(Int_t nXbins, Double_t maxXval, Double_t minXval ,
+ Int_t nYbins, Double_t maxYval, Double_t minYval)
+
+{
+ TString numstr = fName + " Numerator"; //title and name of the
+ //numerator histogram
+ TString denstr = fName + " Denominator";//title and name of the
+ //denominator histogram
+
+ fNumerator = new TH2D(numstr.Data(),numstr.Data(),
+ nXbins,minXval,maxXval,
+ nYbins,minYval,maxYval);
+
+ fDenominator = new TH2D(denstr.Data(),denstr.Data(),
+ nXbins,minXval,maxXval,
+ nYbins,minYval,maxYval);
+
+ fNumerator->Sumw2();
+ fDenominator->Sumw2();
+
+}
+AliHBTFourPartFctn2D::~AliHBTFourPartFctn2D()
+{
+ delete fNumerator;
+ delete fDenominator;
+}
+void AliHBTFourPartFctn2D::
+ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
+{
+ partpair = CheckPair(partpair);
+ trackpair = CheckPair(trackpair);
+ if( partpair && trackpair)
+ {
+ Double_t x,y;
+ GetValues(trackpair,partpair,x,y);
+ fNumerator->Fill(y,x);
+ }
+}
+
+void AliHBTFourPartFctn2D::
+ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair)
+{
+ partpair = CheckPair(partpair);
+ trackpair = CheckPair(trackpair);
+ if( partpair && trackpair)
+ {
+ Double_t x,y;
+ GetValues(trackpair,partpair,x,y);
+ fDenominator->Fill(y,x);
+ }
+
+}
+
--- /dev/null
+//Piotr Skowronski@cern.ch
+
+#ifndef ALIHBTFUNCTION_H
+#define ALIHBTFUNCTION_H
+
+#include "AliHBTParticleCut.h"
+#include "AliHBTPairCut.h"
+#include "AliHBTPair.h"
+
+#include <TH2.h>
+#include <TH3.h>
+
+
+class AliHBTAnalysis;
+
+class AliHBTFunction: public TNamed
+//Abstract base class for HBT functions
+{
+ public:
+ AliHBTFunction();
+ virtual ~AliHBTFunction(){}
+
+ virtual TH1* GetNumerator() =0;
+ virtual TH1* GetDenominator() =0;
+ virtual TH1* GetResult() = 0;
+
+ virtual void Write();
+
+ TH1* GetRatio(Double_t normfactor = 1.0);
+ void Rename(const Char_t * name); //renames the function and histograms ==title is the same that name
+ void Rename(const Char_t * name, const Char_t * title); //renames and retitle the function and histograms
+
+ void SetPairCut(AliHBTPairCut*);
+
+ virtual AliHBTPair* CheckPair(AliHBTPair* pair);
+
+ protected:
+
+ AliHBTPairCut* fPairCut;
+
+ public:
+ ClassDef(AliHBTFunction,1)
+};
+inline AliHBTPair* AliHBTFunction::CheckPair(AliHBTPair* pair)
+{
+ //check if pair and both particles meets the cut criteria
+ if(fPairCut->Pass(pair)) //if the pair is BAD
+ {//it is BAD
+ pair = pair->GetSwapedPair();
+ if(fPairCut->Pass(pair)) //so try reverse combination
+ {
+ return 0x0;//it is BAD as well - so return
+ }
+ }
+ return pair;
+}
+
+
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+class AliHBTTwoPartFctn: public AliHBTFunction
+{
+ public:
+ AliHBTTwoPartFctn(){}
+ virtual ~AliHBTTwoPartFctn(){}
+
+ virtual void ProcessSameEventParticles(AliHBTPair* pair) = 0;
+ virtual void ProcessDiffEventParticles(AliHBTPair* pair) = 0;
+
+
+ protected:
+ public:
+ ClassDef(AliHBTTwoPartFctn,1)
+
+};
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+class AliHBTFourPartFctn: public AliHBTFunction
+{
+ public:
+ AliHBTFourPartFctn(){};
+ virtual ~AliHBTFourPartFctn(){};
+
+ virtual void
+ ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) = 0;
+ virtual void
+ ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair) = 0;
+
+ protected:
+ public:
+ ClassDef(AliHBTFourPartFctn,1)
+
+};
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+
+class AliHBTTwoPartFctn1D: public AliHBTTwoPartFctn
+{
+ public:
+ AliHBTTwoPartFctn1D(Int_t nbins = 100, Double_t maxXval = 0.15, Double_t minXval = 0.0);
+ virtual ~AliHBTTwoPartFctn1D();
+
+
+ TH1* GetNumerator(){return fNumerator;}
+ TH1* GetDenominator(){return fDenominator;}
+
+ void ProcessSameEventParticles(AliHBTPair* pair);
+ void ProcessDiffEventParticles(AliHBTPair* pair);
+
+ protected:
+ //retruns velue to be histogrammed
+ virtual Double_t GetValue(AliHBTPair* pair) = 0;
+
+ TH1D* fNumerator;
+ TH1D* fDenominator;
+
+ public:
+ ClassDef(AliHBTTwoPartFctn1D,1)
+};
+
+inline void
+AliHBTTwoPartFctn1D::ProcessSameEventParticles(AliHBTPair* pair)
+{
+ //Fills the numerator
+ pair = CheckPair(pair);
+ if(pair) fNumerator->Fill(GetValue(pair));
+}
+
+
+inline void
+AliHBTTwoPartFctn1D::ProcessDiffEventParticles(AliHBTPair* pair)
+ {
+ //fills denumerator
+ pair = CheckPair(pair);
+ if(pair) fDenominator->Fill(GetValue(pair));
+
+ }
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliHBTTwoPartFctn2D: public AliHBTTwoPartFctn
+{
+ public:
+ AliHBTTwoPartFctn2D(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0,
+ Int_t nYbins = 200, Double_t maxYval = .15, Double_t minYval =-0.15);
+ ~AliHBTTwoPartFctn2D();
+
+ TH1* GetNumerator(){return fNumerator;}
+ TH1* GetDenominator(){return fDenominator;}
+
+ void ProcessSameEventParticles(AliHBTPair* pair);
+ void ProcessDiffEventParticles(AliHBTPair* pair);
+
+ virtual void GetValues(AliHBTPair* pair, Double_t&, Double_t&) = 0;
+
+ protected:
+ TH2D* fNumerator;
+ TH2D* fDenominator;
+
+ public:
+ ClassDef(AliHBTTwoPartFctn2D,1)
+};
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliHBTTwoPartFctn3D: public AliHBTTwoPartFctn
+{
+ public:
+ AliHBTTwoPartFctn3D(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0,
+ Int_t nYbins = 200, Double_t maxYval = .15, Double_t minYval =-0.15,
+ Int_t nZbins = 200, Double_t maxZval = .15, Double_t minZval =-0.15);
+
+ virtual ~AliHBTTwoPartFctn3D();
+
+ TH1* GetNumerator(){return fNumerator;}
+ TH1* GetDenominator(){return fDenominator;}
+
+ protected:
+ TH3D* fNumerator;
+ TH3D* fDenominator;
+ public:
+ ClassDef(AliHBTTwoPartFctn3D,1)
+};
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+class AliHBTFourPartFctn2D: public AliHBTFourPartFctn
+{
+ public:
+ AliHBTFourPartFctn2D(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0,
+ Int_t nYbins = 200, Double_t maxYval = .15, Double_t minYval =-0.15);
+ ~AliHBTFourPartFctn2D();
+
+ TH1* GetNumerator(){return fNumerator;}
+ TH1* GetDenominator(){return fDenominator;}
+
+ void ProcessSameEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair);
+ void ProcessDiffEventParticles(AliHBTPair* trackpair, AliHBTPair* partpair);
+
+ virtual void GetValues(AliHBTPair*,AliHBTPair*, Double_t&, Double_t&) = 0;
+
+ protected:
+ TH2D* fNumerator;
+ TH2D* fDenominator;
+
+ public:
+ ClassDef(AliHBTFourPartFctn2D,1)
+};
+
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+
+
+#endif
--- /dev/null
+#include "AliHBTPair.h"
+#include "AliHBTParticle.h"
+
+ClassImp(AliHBTPair)
+
+//value of rev
+/************************************************************************/
+AliHBTPair::AliHBTPair(Bool_t rev)
+ {
+//value of rev defines if it is Swaped
+//if you pass kTRUE swpaped pair will NOT be created
+//though you wont be able to get the swaped pair from this pair
+
+ if(!rev) fSwapedPair = new AliHBTPair(kTRUE); //if false create swaped pair
+ else fSwapedPair = 0x0; //if true set the pointer to NULL
+
+ }
+/************************************************************************/
+Double_t AliHBTPair::GetInvMass()
+{
+ if(fInvMassNotCalc)
+ {
+ CalculateInvMassSqr(); //method is inline so we not waste th time for jumping into method
+
+ if(fInvMassSqr<0) fInvMass = TMath::Sqrt(-fInvMassSqr);
+ else fInvMass = TMath::Sqrt(fInvMassSqr);
+
+ fInvMassNotCalc = kFALSE;
+ }
+ return fInvMass;
+}
+/************************************************************************/
+Double_t AliHBTPair::GetQSideCMSLC()
+{
+ //return Q Side in Central Of Mass System in Longitudialy Comoving Frame
+
+ if (fQSideCMSLCNotCalc)
+ {
+ fQSideCMSLC = (fPart1->Px()*fPart2->Py()-fPart2->Px()*fPart1->Py())/GetKt();
+ fQSideCMSLCNotCalc = kFALSE;
+ }
+ return fQSideCMSLC;
+}
+/************************************************************************/
+Double_t AliHBTPair::GetQOutCMSLC()
+{
+ if(fQOutCMSLCNotCalc)
+ {
+ CalculateSums();
+ CalculateDiffs();
+ Double_t k2 = fPxSum*fPxDiff+fPySum*fPyDiff;
+ fQOutCMSLC = 0.5*k2/GetKt();
+ fQOutCMSLCNotCalc = kFALSE;
+ }
+ return fQOutCMSLC;
+}
+/************************************************************************/
+Double_t AliHBTPair::GetQLongCMSLC()
+{
+ if (fQLongCMSLCNotCalc)
+ {
+ CalculateSums();
+ CalculateDiffs();
+ Double_t beta = fPzSum/fESum;
+ Double_t gamma = 1.0/TMath::Sqrt(1.0 - beta*beta);
+ fQLongCMSLC = gamma*( fPzDiff - beta*fEDiff );
+ fQLongCMSLCNotCalc = kFALSE;
+ }
+ return fQLongCMSLC;
+}
+/************************************************************************/
+Double_t AliHBTPair::GetKt()
+{
+ if(fKtNotCalc)
+ {
+ CalculateSums();
+ fKt = 0.5*TMath::Hypot(fPxSum,fPySum);
+ fKtNotCalc = kFALSE;
+ }
+ return fKt;
+}
+/************************************************************************/
+
+Double_t AliHBTPair::GetKStar()
+{
+ if (fKStarNotCalc)
+ {
+
+ CalculateSums();
+
+ Double_t Ptrans = fPxSum*fPxSum + fPySum*fPySum;
+ Double_t Mtrans = fESum*fESum - fPzSum*fPzSum;
+ Double_t Pinv = TMath::Sqrt(Mtrans - Ptrans);
+
+ Double_t Q = ( fPart1->GetMass()*fPart1->GetMass() - fPart2->GetMass()*fPart2->GetMass())/Pinv;
+
+ CalculateQInvL();
+
+ Q = TMath::Sqrt( Q*Q - fQInvL);
+ fKStar = Q/2;
+ fKStarNotCalc = kFALSE;
+ }
+ return fKStar;
+}
+/************************************************************************/
+
+Double_t AliHBTPair::GetQInv()
+{
+ if(fQInvNotCalc)
+ {
+ CalculateQInvL();
+
+ if (fQInvL < 0.) fQInv = TMath::Sqrt(-fQInvL);
+ else fQInv = -TMath::Sqrt( fQInvL);
+
+ fQInvNotCalc = kFALSE;
+ }
+
+ return fQInv;
+}
+
+
+
+
+
+
+
+
--- /dev/null
+#ifndef ALIHBTPAIR_H
+#define ALIHBTPAIR_H
+
+#include <TObject.h>
+
+
+#include "AliHBTParticle.h"
+//class implements pair of particles and taking care of caluclation (almost)
+//all of pair properties (Qinv, InvMass,...)
+//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
+
+class AliHBTPair: public TObject
+{
+ public:
+ AliHBTPair(Bool_t rev = kFALSE); //contructor
+ ~AliHBTPair(){}
+ void SetParticles(AliHBTParticle*,AliHBTParticle*); //sets particles in the pair
+ AliHBTPair* GetSwapedPair() {return fSwapedPair;} //returns pair with swapped particles
+
+ AliHBTParticle* Particle1() const {return fPart1;} //returns pointer to first particle
+ AliHBTParticle* Particle2() const {return fPart2;} //returns pointer to decond particle
+
+ //Center Mass System - Longitudinally Comoving
+
+ Double_t GetInvMass(); //returns invariant mass of the pair
+
+ Double_t GetQSideCMSLC(); //returns Q Side CMS longitudionally co-moving
+ Double_t GetQOutCMSLC(); //returns Q out CMS longitudionally co-moving
+ Double_t GetQLongCMSLC(); //returns Q Long CMS longitudionally co-moving
+
+
+ Double_t GetKt(); //returns K transverse
+ Double_t GetKStar();
+
+ Double_t GetQInv(); //returns Q invariant
+ Double_t GetQSide(); //returns Q side
+ Double_t GetQLong(); //returns Q long
+ Double_t GetQOut(); //returns Q out
+
+
+ protected:
+ AliHBTParticle* fPart1; //pointer to first particle
+ AliHBTParticle* fPart2; //pointer to second particle
+
+ AliHBTPair* fSwapedPair; //pointer to swapped pair
+
+/************************************************************/
+/************CMS (LC) Q's *********************************/
+/************************************************************/
+ //Center Mass System - Longitudinally Comoving
+
+ Double_t fQSideCMSLC; //value of Q side CMS longitudially co-moving
+ Bool_t fQSideCMSLCNotCalc; //flag indicating if fQSideCMSLC is already calculated for this pair
+
+ Double_t fQOutCMSLC; //value of Q out CMS longitudially co-moving
+ Bool_t fQOutCMSLCNotCalc;//flag indicating if fQOutCMSLC is already calculated for this pair
+
+ Double_t fQLongCMSLC; //value of Q long CMS longitudially co-moving
+ Bool_t fQLongCMSLCNotCalc;//flag indicating if fQLongCMSLC is already calculated for this pair
+/************************************************************/
+/************************************************************/
+ Double_t fQInv; //half of differnece of 4-momenta
+ Bool_t fQInvNotCalc;//flag indicating if fQInv is already calculated for this pair
+
+ Double_t fInvMass; //invariant mass
+ Bool_t fInvMassNotCalc;//flag indicating if fInvMass is already calculated for this pair
+
+ Double_t fKt; //K == sum vector of particle's momenta. Kt transverse component
+ Bool_t fKtNotCalc;//flag indicating if fKt is already calculated for this pair
+
+ Double_t fKStar;
+ Bool_t fKStarNotCalc;
+
+ Double_t fPInv; //invariant momentum
+ Double_t fQSide; //Q Side
+ Double_t fOut;//Q Out
+ Double_t fQLong;//Q Long
+
+
+ Double_t fInvMassSqr;//squre of invariant mass
+ Bool_t fMassSqrNotCalc; //
+ void CalculateInvMassSqr();
+
+ Double_t fQInvL;
+ Bool_t fQInvLNotCalc;
+ void CalculateQInvL();
+
+ Double_t fPxSum;
+ Double_t fPySum;
+ Double_t fPzSum;
+ Double_t fESum;
+ Bool_t fSumsNotCalc;
+ void CalculateSums();
+
+ Double_t fPxDiff;
+ Double_t fPyDiff;
+ Double_t fPzDiff;
+ Double_t fEDiff;
+ Bool_t fDiffsNotCalc;
+ void CalculateDiffs();
+
+
+ /***************************************************/
+ void CalculateBase();
+ Bool_t fChanged;
+
+
+ private:
+ public:
+ ClassDef(AliHBTPair,1)
+};
+/****************************************************************/
+inline
+void AliHBTPair::SetParticles(AliHBTParticle* p1,AliHBTParticle* p2)
+{
+ //sets the particle to the pair
+
+ fPart1 = p1;
+ fPart2 = p2;
+ if (fSwapedPair) //if we have Swaped (so we are not)
+ fSwapedPair->SetParticles(p2,p1); //set particles for him too
+
+ // Resel all calculations (flags)
+
+ fChanged = kTRUE;
+ fSumsNotCalc = kTRUE;
+ fDiffsNotCalc = kTRUE;
+ fMassSqrNotCalc = kTRUE;
+ fInvMassNotCalc = kTRUE;
+ fQInvNotCalc = kTRUE;
+ fQSideCMSLCNotCalc = kTRUE;
+ fQOutCMSLCNotCalc = kTRUE;
+ fQLongCMSLCNotCalc = kTRUE;
+ fKtNotCalc = kTRUE;
+ fKStarNotCalc = kTRUE;
+ fQInvLNotCalc = kTRUE;
+
+ //and do nothing until will be asked for
+}
+/****************************************************************/
+inline
+Double_t AliHBTPair::GetQSide()
+{
+ return 0.0;
+}
+/****************************************************************/
+inline
+void AliHBTPair::CalculateInvMassSqr()
+ {
+ if (fMassSqrNotCalc)
+ {
+ CalculateSums();
+
+ Double_t fPart12s= (fPxSum*fPxSum) + (fPySum*fPySum) + (fPzSum*fPzSum);
+
+ fInvMassSqr=fESum*fESum-fPart12s;
+
+ fMassSqrNotCalc = kFALSE;
+ }
+ }
+/****************************************************************/
+inline
+void AliHBTPair::CalculateQInvL()
+ {
+ if (fQInvLNotCalc)
+ {
+ CalculateDiffs();
+ fQInvL = fEDiff*fEDiff - ( fPxDiff*fPxDiff + fPyDiff*fPyDiff + fPzDiff*fPzDiff );
+ fQInvLNotCalc = kFALSE;
+ }
+ }
+/****************************************************************/
+inline
+void AliHBTPair::CalculateSums()
+ {
+ if(fSumsNotCalc)
+ {
+ fPxSum = fPart1->Px()+fPart2->Px();
+ fPySum = fPart1->Py()+fPart2->Py();
+ fPzSum = fPart1->Pz()+fPart2->Pz();
+ fESum = fPart1->Energy() + fPart2->Energy();
+ fSumsNotCalc = kFALSE;
+ }
+ }
+/****************************************************************/
+inline
+void AliHBTPair::CalculateDiffs()
+ {
+ if(fDiffsNotCalc)
+ {
+ fPxDiff = fPart1->Px()-fPart2->Px();
+ fPyDiff = fPart1->Py()-fPart2->Py();
+ fPzDiff = fPart1->Pz()-fPart2->Pz();
+ fEDiff = fPart1->Energy() - fPart2->Energy();
+ fDiffsNotCalc = kFALSE;
+ }
+ }
+
+#endif
--- /dev/null
+#include "AliHBTPairCut.h"
+#include "AliHBTPair.h"
+#include <iostream.h>
+
+
+ClassImp(AliHBTPairCut)
+const Int_t AliHBTPairCut::fkgMaxCuts = 50;
+/**********************************************************/
+
+AliHBTPairCut::AliHBTPairCut()
+{
+//constructor
+ fFirstPartCut = new AliHBTEmptyParticleCut(); //empty cuts
+ fSecondPartCut= new AliHBTEmptyParticleCut(); //empty cuts
+
+ fCuts = new AliHbtBasePairCut*[fkgMaxCuts];
+ fNCuts = 0;
+}
+/**********************************************************/
+
+AliHBTPairCut::AliHBTPairCut(const AliHBTPairCut& in)
+{
+ //copy constructor
+ fCuts = new AliHbtBasePairCut*[fkgMaxCuts];
+ fNCuts = in.fNCuts;
+
+ fFirstPartCut = (AliHBTParticleCut*)in.fFirstPartCut->Clone();
+ fSecondPartCut = (AliHBTParticleCut*)in.fSecondPartCut->Clone();
+
+ for (Int_t i = 0;i<fNCuts;i++)
+ {
+ fCuts[i] = (AliHbtBasePairCut*)in.fCuts[i]->Clone();//create new object (clone) and rember pointer to it
+ }
+}
+/**********************************************************/
+
+AliHBTPairCut::~AliHBTPairCut()
+{
+//destructor
+ for (Int_t i = 0;i<fNCuts;i++)
+ {
+ delete fCuts[i];
+ }
+ delete []fCuts;
+}
+/**********************************************************/
+
+/**********************************************************/
+
+void AliHBTPairCut::AddBasePairCut(AliHbtBasePairCut* basecut)
+ {
+ //adds the base pair cut (cut on one value)
+
+ if (!basecut) return;
+ if( fNCuts == (fkgMaxCuts-1) )
+ {
+ Warning("AddBasePairCut","Not enough place for another cut");
+ return;
+ }
+ fCuts[fNCuts++]=basecut;
+ }
+/**********************************************************/
+
+Bool_t AliHBTPairCut::Pass(AliHBTPair* pair)
+{
+//methods which checks if given pair meets all criteria of the cut
+//if it meets returns FALSE
+//if NOT returns TRUE
+ if(!pair)
+ {
+ Warning("Pass()","No Pasaran! We never accept NULL pointers");
+ return kTRUE;
+ }
+
+ //check particle's cuts
+ if( ( fFirstPartCut->Pass( pair->Particle1() ) ) ||
+ ( fSecondPartCut->Pass(pair->Particle2() ) ) )
+ {
+ return kTRUE;
+ }
+
+ //examine all base pair cuts
+ for (Int_t i = 0;i<fNCuts;i++)
+ {
+ if ( (fCuts[i]->Pass(pair)) ) return kTRUE; //if one of the cuts reject, then reject
+ }
+
+// cout<<"passed "<<pair->Particle1()->GetPdgCode()<<" "<<pair->Particle2()->GetPdgCode()<<endl;
+ return kFALSE;
+}
+/**********************************************************/
+
+void AliHBTPairCut::SetFirstPartCut(AliHBTParticleCut* cut)
+{
+ if(!cut)
+ {
+ Error("SetFirstPartCut","argument is NULL");
+ return;
+ }
+ delete fFirstPartCut;
+ fFirstPartCut = (AliHBTParticleCut*)cut->Clone();
+
+}
+/**********************************************************/
+
+void AliHBTPairCut::SetSecondPartCut(AliHBTParticleCut* cut)
+{
+ if(!cut)
+ {
+ Error("SetSecondPartCut","argument is NULL");
+ return;
+ }
+ delete fSecondPartCut;
+ fSecondPartCut = (AliHBTParticleCut*)cut->Clone();
+}
+/**********************************************************/
+
+void AliHBTPairCut::SetPartCut(AliHBTParticleCut* cut)
+{
+//sets the the same cut on both particles
+ if(!cut)
+ {
+ Error("SetFirstPartCut","argument is NULL");
+ return;
+ }
+ delete fFirstPartCut;
+ fFirstPartCut = (AliHBTParticleCut*)cut->Clone();
+
+ delete fSecondPartCut; //even if null should not be harmful
+ fSecondPartCut = fFirstPartCut;
+}
+/**********************************************************/
+
+void AliHBTPairCut::SetQInvRange(Double_t min, Double_t max)
+{
+ AliHBTQInvCut* cut= (AliHBTQInvCut*)FindCut(kHbtPairCutPropQInv);
+ if(cut) cut->SetRange(min,max);
+ else fCuts[fNCuts++] = new AliHBTQInvCut(min,max);
+
+}
+
+AliHbtBasePairCut* AliHBTPairCut::FindCut(AliHBTPairCutProperty property)
+{
+ for (Int_t i = 0;i<fNCuts;i++)
+ {
+ if (fCuts[i]->GetProperty() == property)
+ return fCuts[i]; //we found the cut we were searching for
+ }
+
+ return 0x0; //we did not found this cut
+
+}
+/**********************************************************/
+
+void AliHBTPairCut::Streamer(TBuffer &b)
+{
+ // Stream all objects in the array to or from the I/O buffer.
+
+ UInt_t R__s, R__c;
+ if (b.IsReading())
+ {
+ Version_t v = b.ReadVersion(&R__s, &R__c);
+ TObject::Streamer(b);
+ b >> fFirstPartCut;
+ b >> fSecondPartCut;
+ b >> fNCuts;
+ for (Int_t i = 0;i<fNCuts;i++)
+ {
+ b >> fCuts[i];
+ }
+ b.CheckByteCount(R__s, R__c,AliHBTPairCut::IsA());
+ }
+ else
+ {
+ R__c = b.WriteVersion(AliHBTPairCut::IsA(), kTRUE);
+ TObject::Streamer(b);
+ b << fFirstPartCut;
+ b << fSecondPartCut;
+ b << fNCuts;
+ for (Int_t i = 0;i<fNCuts;i++)
+ {
+ b << fCuts[i];
+ }
+ b.SetByteCount(R__c, kTRUE);
+ }
+}
+
+ClassImp(AliHBTEmptyPairCut)
+
+void AliHBTEmptyPairCut::Streamer(TBuffer &b)
+ {
+ AliHBTPairCut::Streamer(b);
+ }
+
+ClassImp(AliHbtBasePairCut)
+
+ClassImp(AliHBTQInvCut)
--- /dev/null
+//Piotr Skowronski@cern.ch
+//Class implemnts cut on the pair of particles
+//
+//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
+
+#ifndef ALIHBTPAIRCUT_H
+#define ALIHBTPAIRCUT_H
+
+
+#include "AliHBTParticleCut.h"
+#include "AliHBTPair.h"
+
+class AliHbtBasePairCut;
+
+enum AliHBTPairCutProperty
+ {
+ kHbtPairCutPropQInv, //Q invariant
+ kHbtPairCutPropNone
+ };
+
+class AliHBTPairCut: public TObject
+{
+ public:
+ AliHBTPairCut();
+ AliHBTPairCut(const AliHBTPairCut&);
+
+ virtual ~AliHBTPairCut();
+ virtual Bool_t Pass(AliHBTPair*);
+
+ void SetFirstPartCut(AliHBTParticleCut*); //sets the cut on the first particle
+ void SetSecondPartCut(AliHBTParticleCut*); //sets the cut on the first particle
+
+ void SetPartCut(AliHBTParticleCut*);//sets the the same cut on both particles
+
+ void AddBasePairCut(AliHbtBasePairCut*);
+
+ void SetQInvRange(Double_t min, Double_t max);
+
+ protected:
+ AliHBTParticleCut* fFirstPartCut;
+ AliHBTParticleCut* fSecondPartCut;
+
+ AliHbtBasePairCut** fCuts; //!
+ Int_t fNCuts;
+
+
+ AliHbtBasePairCut* FindCut(AliHBTPairCutProperty);
+ private:
+ static const Int_t fkgMaxCuts;
+ public:
+ ClassDef(AliHBTPairCut,1)
+
+};
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliHBTEmptyPairCut: public AliHBTPairCut
+{
+//Empty - it passes possitively all particles - it means returns always False
+//Class describing cut on pairs of particles
+ public:
+ AliHBTEmptyPairCut(){};
+ AliHBTEmptyPairCut(const AliHBTEmptyPairCut&){};
+ virtual ~AliHBTEmptyPairCut(){};
+
+ Bool_t Pass(AliHBTPair*) {return kFALSE;} //accpept everything
+
+ ClassDef(AliHBTEmptyPairCut,1)
+
+};
+
+
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliHbtBasePairCut: public TObject
+ {
+ //This class defines the range of some property - pure virtual
+ //Property is coded by AliHBTCutTypes type
+
+ public:
+
+ AliHbtBasePairCut(Double_t min = 0.0, Double_t max = 0.0, AliHBTPairCutProperty prop= kHbtPairCutPropNone):
+ fMin(min),fMax(max),fProperty(prop){}
+
+ virtual ~AliHbtBasePairCut(){}
+
+ Bool_t Pass(AliHBTPair*);
+
+ void SetRange(Double_t min, Double_t max){fMin = min; fMax = max;}
+
+ void SetMinimum(Double_t min){fMin = min;}
+ void SetMaximum(Double_t max){fMax = max;}
+
+ Double_t GetMinimum() const {return fMin;}
+ Double_t GetMaximum() const {return fMax;}
+
+ AliHBTPairCutProperty GetProperty() const {return fProperty;}
+
+ protected:
+ virtual Double_t GetValue(AliHBTPair*) = 0;
+
+ Double_t fMin;
+ Double_t fMax;
+
+ AliHBTPairCutProperty fProperty;
+
+ private:
+ public:
+ ClassDef(AliHbtBasePairCut,1)
+
+ };
+
+inline Bool_t
+AliHbtBasePairCut::Pass(AliHBTPair* pair)
+ {
+ Double_t value = GetValue(pair);
+ if ( (value > fMin) && (value <fMax ) ) return kFALSE; //accepted
+ else return kTRUE; //rejected
+ }
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+class AliHBTQInvCut: public AliHbtBasePairCut
+ {
+ public:
+ AliHBTQInvCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBasePairCut(min,max,kHbtPairCutPropQInv){}
+ virtual ~AliHBTQInvCut(){}
+ protected:
+ virtual Double_t GetValue(AliHBTPair* pair){return pair->GetQInv();}
+ private:
+ public:
+ ClassDef(AliHBTQInvCut,1)
+ };
+
+#endif
--- /dev/null
+//Simplified TParticle class
+
+
+#include "AliHBTParticle.h"
+
+#include <TParticle.h>
+
+ClassImp(AliHBTParticle)
+
+//______________________________________________________________________________
+AliHBTParticle::AliHBTParticle():
+ fPdgCode(0), fPx(0), fPy(0),fPz(0),fE(0), fVx(0), fVy(0),fVz(0),fVt(0)
+{//empty particle
+}
+
+
+//______________________________________________________________________________
+AliHBTParticle::AliHBTParticle(Int_t pdg, Double_t px, Double_t py, Double_t pz, Double_t etot,
+ Double_t vx, Double_t vy, Double_t vz, Double_t time):
+ fPdgCode(pdg), fPx(px), fPy(py),fPz(pz),fE(etot),
+ fVx(vx), fVy(vy),fVz(vz),fVt(time)
+{
+//mormal constructor
+
+ if (GetPDG()) {
+ fCalcMass = GetPDG()->Mass();
+ } else {
+ Double_t a2 = fE*fE -fPx*fPx -fPy*fPy -fPz*fPz;
+ if (a2 >= 0) fCalcMass = TMath::Sqrt(a2);
+ else fCalcMass = -TMath::Sqrt(-a2);
+ }
+}
+
+//______________________________________________________________________________
+AliHBTParticle::AliHBTParticle(const TParticle &p):
+ fPdgCode(p.GetPdgCode()),
+ fPx(p.Px()),fPy(p.Py()),fPz(p.Pz()),fE(p.Energy()),
+ fVx(p.Vx()),fVy(p.Vy()),fVz(p.Vz()),fVt(p.T())
+{
+ //all copied in the initialization
+}
+
+//______________________________________________________________________________
+const Char_t* AliHBTParticle::GetName() const
+{
+ static char def[4] = "XXX";
+ const TParticlePDG *ap = TDatabasePDG::Instance()->GetParticle(fPdgCode);
+ if (ap) return ap->GetName();
+ else return def;
+}
+
+
+//______________________________________________________________________________
+
+
--- /dev/null
+#ifndef ALIHBTPARTICLE
+#define ALIHBTPARTICLE
+
+//Ali HBT Particle: simplified class TParticle
+//Simplified in order to minimize the size of object
+// - we want to keep a lot of such a objects in memory
+
+#include <TObject.h>
+#include <TLorentzVector.h>
+#include <TMath.h>
+#include <TDatabasePDG.h>
+
+class TParticle;
+
+class AliHBTParticle : public TObject
+{
+
+public:
+ // ****** constructors and destructor
+ AliHBTParticle();
+
+ AliHBTParticle(Int_t pdg, Double_t px, Double_t py, Double_t pz, Double_t etot,
+ Double_t vx, Double_t vy, Double_t vz, Double_t time);
+
+ AliHBTParticle(const TParticle& p);
+
+ virtual ~AliHBTParticle(){};
+
+
+ Int_t GetPdgCode () const { return fPdgCode; }
+ Double_t GetCalcMass () const { return fCalcMass; }
+ Double_t GetMass () { return GetPDG()->Mass();}
+
+ TParticlePDG* GetPDG (){return TDatabasePDG::Instance()->GetParticle(fPdgCode);}
+
+ Int_t Beauty () { return GetPDG()->Beauty(); }
+ Int_t Charm () { return GetPDG()->Charm(); }
+ Int_t Strangeness () { return GetPDG()->Strangeness();}
+ void ProductionVertex(TLorentzVector &v) { v.SetXYZT(fVx,fVy,fVz,fVt);}
+
+
+ Double_t Vx () const { return fVx;}
+ Double_t Vy () const { return fVy;}
+ Double_t Vz () const { return fVz;}
+ Double_t T () const { return fVt;}
+
+ Double_t Px () const { return fPx; } //X coordinate of the momentum
+ Double_t Py () const { return fPy; } //Y coordinate of the momentum
+ Double_t Pz () const { return fPz; } //Z coordinate of the momentum
+ Double_t P () const //momentum
+ { return TMath::Sqrt(fPx*fPx+fPy*fPy+fPz*fPz); }
+
+ void Momentum(TLorentzVector &v) { v.SetPxPyPzE(fPx,fPy,fPz,fE);}
+
+ Double_t Pt () const //transverse momentum
+ { return TMath::Sqrt(fPx*fPx+fPy*fPy); }
+ Double_t Energy() const { return fE; }
+
+ //Pseudo Rapidity
+ Double_t Eta () const { if (P() != fPz) return 0.5*TMath::Log((P()+fPz)/(P()-fPz));
+ else return 1.e30;}
+
+ //Rapidity
+ Double_t Y () const { if (fE != fPz) return 0.5*TMath::Log((fE+fPz)/(fE-fPz));
+ else return 1.e30;}
+
+ Double_t Phi () const { return kPI+TMath::ATan2(-fPy,-fPx); }
+
+ Double_t Theta () const { return (fPz==0)?kPI/2:TMath::ACos(fPz/P()); }
+
+ // setters
+
+
+ void SetMomentum(Double_t px, Double_t py, Double_t pz, Double_t e)
+ {fPx=px; fPy=py; fPz=pz; fE=e;}
+ void SetMomentum(const TLorentzVector& p)
+ {SetMomentum(p.Px(),p.Py(),p.Pz(),p.Energy());}
+
+ void SetProductionVertex(Double_t vx, Double_t vy, Double_t vz, Double_t t)
+ {fVx=vx; fVy=vy; fVz=vz; fVt=t;}
+ void SetProductionVertex(const TLorentzVector& v)
+ {SetProductionVertex(v.X(),v.Y(),v.Z(),v.T());}
+ const Char_t* GetName() const;
+
+protected:
+
+ Int_t fPdgCode; // PDG code of the particle
+ Double_t fCalcMass; // Calculated mass
+
+ Double_t fPx; // x component of momentum
+ Double_t fPy; // y component of momentum
+ Double_t fPz; // z component of momentum
+ Double_t fE; // Energy
+
+ Double_t fVx; // x of production vertex
+ Double_t fVy; // y of production vertex
+ Double_t fVz; // z of production vertex
+ Double_t fVt; // t of production vertex
+
+public:
+ ClassDef(AliHBTParticle,1) // TParticle vertex particle information
+};
+
+#endif
--- /dev/null
+#include "AliHBTParticleCut.h"
+
+#include <iostream.h>
+
+ClassImp(AliHBTParticleCut)
+const Int_t AliHBTParticleCut::fkgMaxCuts = 50;
+/******************************************************************/
+
+AliHBTParticleCut::AliHBTParticleCut()
+ {
+ fCuts = new AliHbtBaseCut* [fkgMaxCuts];//last property in the property
+ //property enum => defines number of properties
+ fNCuts = 0;
+ fPID = 0;
+ }
+/******************************************************************/
+
+AliHBTParticleCut::AliHBTParticleCut(const AliHBTParticleCut& in)
+{
+ fCuts = new AliHbtBaseCut* [fkgMaxCuts];//last property in the property
+ //property enum => defines number of properties
+ fNCuts = in.fNCuts;
+ fPID = in.fPID;
+ for (Int_t i = 0;i<fNCuts;i++)
+ {
+ fCuts[i] = (AliHbtBaseCut*)in.fCuts[i]->Clone();//create new object (clone) and rember pointer to it
+ }
+}
+/******************************************************************/
+
+AliHBTParticleCut::~AliHBTParticleCut()
+{
+ for (Int_t i = 0;i<fNCuts;i++)
+ {
+ delete fCuts[i];
+ }
+ delete []fCuts;
+}
+/******************************************************************/
+
+Bool_t AliHBTParticleCut::Pass(AliHBTParticle* p)
+{
+//method checks all the cuts that are set (in the list)
+//If any of the baseCuts rejects particle False(rejection) is returned
+ if(!p)
+ {
+ Warning("Pass()","No Pasaran! We never accept NULL pointers");
+ return kTRUE;
+ }
+ if( (p->GetPdgCode() != fPID) && ( fPID != 0)) return kTRUE;
+
+ for (Int_t i = 0;i<fNCuts;i++)
+ {
+ if ( (fCuts[i]->Pass(p)) ) return kTRUE; //if one of the cuts rejects, then reject
+ }
+ return kFALSE;
+}
+/******************************************************************/
+
+void AliHBTParticleCut::AddBasePartCut(AliHbtBaseCut* basecut)
+{
+ //adds the base pair cut (cut on one value)
+
+ if (!basecut) return;
+ if( fNCuts == (fkgMaxCuts-1) )
+ {
+ Warning("AddBasePartCut","Not enough place for another cut");
+ return;
+ }
+ fCuts[fNCuts++]=basecut;
+
+}
+
+/******************************************************************/
+AliHbtBaseCut* AliHBTParticleCut::FindCut(AliHBTCutProperty property)
+{
+
+ for (Int_t i = 0;i<fNCuts;i++)
+ {
+ if (fCuts[i]->GetProperty() == property)
+ return fCuts[i]; //we found the cut we were searching for
+ }
+
+ return 0x0; //we did not found this cut
+
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetMomentumRange(Double_t min, Double_t max)
+{
+ AliHBTMomentumCut* cut= (AliHBTMomentumCut*)FindCut(kHbtP);
+ if(cut) cut->SetRange(min,max);
+ else fCuts[fNCuts++] = new AliHBTMomentumCut(min,max);
+}
+/******************************************************************/
+
+
+void AliHBTParticleCut::SetPtRange(Double_t min, Double_t max)
+{
+ AliHBTPtCut* cut= (AliHBTPtCut*)FindCut(kHbtPt);
+ if(cut) cut->SetRange(min,max);
+ else fCuts[fNCuts++] = new AliHBTPtCut(min,max);
+
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetEnergyRange(Double_t min, Double_t max)
+{
+ AliHBTEnergyCut* cut= (AliHBTEnergyCut*)FindCut(kHbtE);
+ if(cut) cut->SetRange(min,max);
+ else fCuts[fNCuts++] = new AliHBTEnergyCut(min,max);
+
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetRapidityRange(Double_t min, Double_t max)
+{
+ AliHbtBaseCut* cut = FindCut(kHbtRapidity);
+ if(cut) cut->SetRange(min,max);
+ else fCuts[fNCuts++] = new AliHBTRapidityCut(min,max);
+
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetPseudoRapidityRange(Double_t min, Double_t max)
+{
+ AliHbtBaseCut* cut = FindCut(kHbtPseudoRapidity);
+ if(cut) cut->SetRange(min,max);
+ else fCuts[fNCuts++] = new AliHBTPseudoRapidityCut(min,max);
+
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetPxRange(Double_t min, Double_t max)
+{
+ AliHbtBaseCut* cut = FindCut(kHbtPx);
+ if(cut) cut->SetRange(min,max);
+ else fCuts[fNCuts++] = new AliHBTPxCut(min,max);
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetPyRange(Double_t min, Double_t max)
+{
+ AliHbtBaseCut* cut = FindCut(kHbtPy);
+ if(cut) cut->SetRange(min,max);
+ else fCuts[fNCuts++] = new AliHBTPyCut(min,max);
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetPzRange(Double_t min, Double_t max)
+{
+ AliHbtBaseCut* cut = FindCut(kHbtPz);
+ if(cut) cut->SetRange(min,max);
+ else fCuts[fNCuts++] = new AliHBTPzCut(min,max);
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetPhiRange(Double_t min, Double_t max)
+{
+ AliHbtBaseCut* cut = FindCut(kHbtPhi);
+ if(cut) cut->SetRange(min,max);
+ else fCuts[fNCuts++] = new AliHBTPhiCut(min,max);
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetThetaRange(Double_t min, Double_t max)
+{
+ AliHbtBaseCut* cut = FindCut(kHbtTheta);
+ if(cut) cut->SetRange(min,max);
+ else fCuts[fNCuts++] = new AliHBTThetaCut(min,max);
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetVxRange(Double_t min, Double_t max)
+{
+ AliHbtBaseCut* cut = FindCut(kHbtVx);
+ if(cut) cut->SetRange(min,max);
+ else fCuts[fNCuts++] = new AliHBTVxCut(min,max);
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetVyRange(Double_t min, Double_t max)
+{
+ AliHbtBaseCut* cut = FindCut(kHbtVy);
+ if(cut) cut->SetRange(min,max);
+ else fCuts[fNCuts++] = new AliHBTVyCut(min,max);
+}
+/******************************************************************/
+
+void AliHBTParticleCut::SetVzRange(Double_t min, Double_t max)
+{
+ AliHbtBaseCut* cut = FindCut(kHbtVz);
+ if(cut) cut->SetRange(min,max);
+ else fCuts[fNCuts++] = new AliHBTVzCut(min,max);
+}
+
+/******************************************************************/
+void AliHBTParticleCut::Streamer(TBuffer &b)
+{
+ // Stream all objects in the array to or from the I/O buffer.
+
+ UInt_t R__s, R__c;
+ if (b.IsReading())
+ {
+ Version_t v = b.ReadVersion(&R__s, &R__c);
+ TObject::Streamer(b);
+ b >> fPID;
+ b >> fNCuts;
+ for (Int_t i = 0;i<fNCuts;i++)
+ {
+ b >> fCuts[i];
+ }
+ b.CheckByteCount(R__s, R__c,AliHBTParticleCut::IsA());
+ }
+ else
+ {
+ R__c = b.WriteVersion(AliHBTParticleCut::IsA(), kTRUE);
+ TObject::Streamer(b);
+ b << fPID;
+ b << fNCuts;
+ for (Int_t i = 0;i<fNCuts;i++)
+ {
+ b << fCuts[i];
+ }
+ b.SetByteCount(R__c, kTRUE);
+ }
+}
+
+void AliHBTParticleCut::Print(void)
+{
+ cout<<"Printing AliHBTParticleCut, this = "<<this<<endl;
+ cout<<"fPID "<<fPID<<endl;
+ cout<<"fNCuts "<<fNCuts <<endl;
+ for (Int_t i = 0;i<fNCuts;i++)
+ {
+ cout<<" fCuts["<<i<<"] "<<fCuts[i]<<endl<<" ";
+ fCuts[i]->Print();
+ }
+}
+
+/******************************************************************/
+/******************************************************************/
+
+ClassImp(AliHBTEmptyParticleCut)
+void AliHBTEmptyParticleCut::Streamer(TBuffer &b)
+ {
+ AliHBTParticleCut::Streamer(b);
+ }
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp(AliHbtBaseCut)
+void AliHbtBaseCut::Print(void)
+{
+ cout<<"fMin="<<fMin <<", fMax=" <<fMax<<" ";
+ PrintProperty();
+}
+void AliHbtBaseCut::PrintProperty(void)
+{
+ switch (fProperty)
+ {
+ case kHbtP:
+ cout<<"kHbtP"; break;
+ case kHbtPt:
+ cout<<"kHbtPt"; break;
+ case kHbtE:
+ cout<<"kHbtE"; break;
+ case kHbtRapidity:
+ cout<<"kHbtRapidity"; break;
+ case kHbtPseudoRapidity:
+ cout<<"kHbtPseudoRapidity"; break;
+ case kHbtPx:
+ cout<<"kHbtPx"; break;
+ case kHbtPy:
+ cout<<"kHbtPy"; break;
+ case kHbtPz:
+ cout<<"kHbtPz"; break;
+ case kHbtPhi:
+ cout<<"kHbtPhi"; break;
+ case kHbtTheta:
+ cout<<"kHbtTheta"; break;
+ case kHbtVx:
+ cout<<"kHbtVx"; break;
+ case kHbtVy:
+ cout<<"kHbtVy"; break;
+ case kHbtVz:
+ cout<<"kHbtVz"; break;
+ case kHbtNone:
+ cout<<"kHbtNone"; break;
+ default:
+ cout<<"Property Not Found";
+ }
+ cout<<endl;
+}
+ClassImp( AliHBTMomentumCut )
+
+ClassImp( AliHBTPtCut )
+ClassImp( AliHBTEnergyCut )
+ClassImp( AliHBTRapidityCut )
+ClassImp( AliHBTPseudoRapidityCut )
+ClassImp( AliHBTPxCut )
+ClassImp( AliHBTPyCut )
+ClassImp( AliHBTPzCut )
+ClassImp( AliHBTPhiCut )
+ClassImp( AliHBTThetaCut )
+ClassImp( AliHBTVxCut )
+ClassImp( AliHBTVyCut )
+ClassImp( AliHBTVzCut )
+
+
--- /dev/null
+//Piotr Skowronski@cern.ch
+//Classes for single particle cuts
+//User should use only AliHBTParticleCut, eventually EmptyCut which passes all particles
+//There is all interface for setting cuts on all particle properties
+//The main method is Pass - which returns
+// True in order to reject particle
+// False in case it meets all the criteria of the given cut
+//
+//User should create (and also destroy) cuts himself
+//and then pass them to the Analysis or Function by a proper method
+//
+//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
+
+
+#ifndef ALIHBTPARTICLECUT_H
+#define ALIHBTPARTICLECUT_H
+
+#include <TObject.h>
+#include "AliHBTParticle.h"
+
+
+class AliHBTEmptyParticleCut;
+class AliHBTParticleCut;
+class AliHBTPairCut;
+class AliHBTPair;
+class AliHbtBaseCut;
+
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+enum AliHBTCutProperty
+ {
+//codes particle property
+ kHbtP, //Momentum
+ kHbtPt, //Transverse momentum
+ kHbtE, //Energy
+ kHbtRapidity, //
+ kHbtPseudoRapidity,
+ kHbtPx, //X coordinate of the momentum
+ kHbtPy, //Y coordinate of the momentum
+ kHbtPz, //Z coordinate of the momentum
+ kHbtPhi,//angle
+ kHbtTheta,//angle
+ kHbtVx, // vertex X coordinate
+ kHbtVy, // vertex Y coordinate
+ kHbtVz, // vertex Z coordinate
+//_____________________________
+ kHbtNone
+ };
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliHBTParticleCut: public TObject
+{
+//Class describing cut on pairs of particles
+ public:
+ AliHBTParticleCut();
+ AliHBTParticleCut(const AliHBTParticleCut&);
+ virtual ~AliHBTParticleCut();
+
+ virtual Bool_t Pass(AliHBTParticle*);
+
+ void AddBasePartCut(AliHbtBaseCut*);
+
+ Int_t GetPID() const { return fPID;}
+ void SetPID(Int_t pid){fPID=pid;}
+ void SetMomentumRange(Double_t min, Double_t max);
+ void SetPRange(Double_t min, Double_t max){SetMomentumRange(min,max);}
+ void SetPtRange(Double_t min, Double_t max);
+ void SetEnergyRange(Double_t min, Double_t max);
+ void SetRapidityRange(Double_t min, Double_t max);
+ void SetYRange(Double_t min, Double_t max){SetRapidityRange(min,max);}
+ void SetPseudoRapidityRange(Double_t min, Double_t max);
+ void SetPxRange(Double_t min, Double_t max);
+ void SetPyRange(Double_t min, Double_t max);
+ void SetPzRange(Double_t min, Double_t max);
+ void SetPhiRange(Double_t min, Double_t max);
+ void SetThetaRange(Double_t min, Double_t max);
+ void SetVxRange(Double_t min, Double_t max);
+ void SetVyRange(Double_t min, Double_t max);
+ void SetVzRange(Double_t min, Double_t max);
+
+ void Print(void);
+ protected:
+
+ AliHbtBaseCut* FindCut(AliHBTCutProperty);
+
+ AliHbtBaseCut ** fCuts;//! Array with cuts
+ Int_t fNCuts;
+
+ Int_t fPID; //particle PID - if=0 (rootino) all pids are accepted
+
+ private:
+ static const Int_t fkgMaxCuts;
+ public:
+ ClassDef(AliHBTParticleCut,1)
+
+};
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliHBTEmptyParticleCut: public AliHBTParticleCut
+{
+//Empty - it passes possitively all particles - it means returns always False
+//Class describing cut on pairs of particles
+ public:
+ AliHBTEmptyParticleCut(){};
+ virtual ~AliHBTEmptyParticleCut(){};
+
+ Bool_t Pass(AliHBTParticle*) {return kFALSE;} //accpept everything
+
+ ClassDef(AliHBTEmptyParticleCut,1)
+
+};
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+class AliHbtBaseCut: public TObject
+ {
+ //This class defines the range of some property - pure virtual
+ //Property is coded by AliHBTCutTypes type
+
+ public:
+
+ AliHbtBaseCut(Double_t min = 0.0, Double_t max = 0.0,AliHBTCutProperty prop = kHbtNone):
+ fProperty(prop),fMin(min),fMax(max){}
+
+ virtual ~AliHbtBaseCut(){}
+
+ Bool_t Pass(AliHBTParticle *);
+
+ void SetRange(Double_t min, Double_t max){fMin = min; fMax = max;}
+
+ void SetMinimum(Double_t min){fMin = min;}
+ void SetMaximum(Double_t max){fMax = max;}
+
+ Double_t GetMinimum() const {return fMin;}
+ Double_t GetMaximum() const {return fMax;}
+
+ AliHBTCutProperty GetProperty() const {return fProperty;}
+ void Print(void);
+ protected:
+ virtual Double_t GetValue(AliHBTParticle *) = 0;
+
+ AliHBTCutProperty fProperty;
+ Double_t fMin;
+ Double_t fMax;
+ private:
+ void PrintProperty(void);
+ public:
+ ClassDef(AliHbtBaseCut,1)
+
+ };
+
+inline Bool_t
+AliHbtBaseCut::Pass(AliHBTParticle *p)
+ {
+ if ( (GetValue(p) > fMin) && (GetValue(p) <fMax ) ) return kFALSE; //accepted
+ else return kTRUE; //rejected
+ }
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+
+class AliHBTMomentumCut: public AliHbtBaseCut
+ {
+ public:
+ AliHBTMomentumCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtP){}
+ virtual ~AliHBTMomentumCut(){}
+ protected:
+ Double_t GetValue(AliHBTParticle * p){return p->P();}
+ public:
+ ClassDef(AliHBTMomentumCut,1)
+
+ };
+
+
+class AliHBTPtCut: public AliHbtBaseCut
+ {
+ public:
+ AliHBTPtCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtPt){}
+ virtual ~AliHBTPtCut(){}
+ protected:
+ Double_t GetValue(AliHBTParticle * p){return p->Pt();}
+ public:
+ ClassDef(AliHBTPtCut,1)
+
+ };
+
+
+class AliHBTEnergyCut: public AliHbtBaseCut
+ {
+ public:
+ AliHBTEnergyCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtE){}
+ virtual ~AliHBTEnergyCut(){}
+ protected:
+ Double_t GetValue(AliHBTParticle * p){return p->Energy();}
+ public:
+ ClassDef(AliHBTEnergyCut,1)
+
+ };
+
+class AliHBTRapidityCut: public AliHbtBaseCut
+ {
+ public:
+ AliHBTRapidityCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtRapidity){}
+ virtual ~AliHBTRapidityCut(){}
+ protected:
+ Double_t GetValue(AliHBTParticle * p){return p->Y();}
+ public:
+ ClassDef(AliHBTRapidityCut,1)
+
+ };
+
+class AliHBTPseudoRapidityCut: public AliHbtBaseCut
+ {
+ public:
+ AliHBTPseudoRapidityCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtPseudoRapidity){}
+ virtual ~AliHBTPseudoRapidityCut(){}
+ protected:
+ Double_t GetValue(AliHBTParticle * p){return p->Eta();}
+ public:
+ ClassDef(AliHBTPseudoRapidityCut,1)
+
+ };
+
+class AliHBTPxCut: public AliHbtBaseCut
+ {
+ public:
+ AliHBTPxCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtPx){}
+ virtual ~AliHBTPxCut(){}
+ protected:
+ Double_t GetValue(AliHBTParticle * p){return p->Px();}
+ public:
+ ClassDef(AliHBTPxCut,1)
+
+ };
+
+class AliHBTPyCut: public AliHbtBaseCut
+ {
+ public:
+ AliHBTPyCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtPy){}
+ virtual ~AliHBTPyCut(){}
+ protected:
+ Double_t GetValue(AliHBTParticle * p){return p->Py();}
+ public:
+ ClassDef(AliHBTPyCut,1)
+
+ };
+
+
+class AliHBTPzCut: public AliHbtBaseCut
+ {
+ public:
+ AliHBTPzCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtPz){}
+ virtual ~AliHBTPzCut(){}
+ protected:
+ Double_t GetValue(AliHBTParticle * p){return p->Pz();}
+ public:
+ ClassDef(AliHBTPzCut,1)
+
+ };
+
+class AliHBTPhiCut: public AliHbtBaseCut
+ {
+ public:
+ AliHBTPhiCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtPhi){}
+ virtual ~AliHBTPhiCut(){}
+ protected:
+ Double_t GetValue(AliHBTParticle * p){return p->Phi();}
+ public:
+ ClassDef(AliHBTPhiCut,1)
+
+ };
+
+class AliHBTThetaCut: public AliHbtBaseCut
+ {
+ public:
+ AliHBTThetaCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtTheta){}
+ virtual ~AliHBTThetaCut(){}
+ protected:
+ Double_t GetValue(AliHBTParticle * p){return p->Theta();}
+ public:
+ ClassDef(AliHBTThetaCut,1)
+
+ };
+
+class AliHBTVxCut: public AliHbtBaseCut
+ {
+ //Cut of the X coordinate of the vertex position
+ public:
+ AliHBTVxCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtVx){}
+ virtual ~AliHBTVxCut(){}
+ protected:
+ Double_t GetValue(AliHBTParticle * p){return p->Vx();} //retruns value of the vertex
+ public:
+ ClassDef(AliHBTVxCut,1)
+
+ };
+
+
+class AliHBTVyCut: public AliHbtBaseCut
+ {
+ //Cut of the X coordinate of the vertex position
+ public:
+ AliHBTVyCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtVy){}
+ virtual ~AliHBTVyCut(){}
+ protected:
+ Double_t GetValue(AliHBTParticle * p){return p->Vy();} //retruns value of the vertex
+ public:
+ ClassDef(AliHBTVyCut,1)
+
+ };
+
+class AliHBTVzCut: public AliHbtBaseCut
+ {
+ //Cut of the X coordinate of the vertex position
+ public:
+ AliHBTVzCut(Double_t min = 0.0, Double_t max = 0.0):AliHbtBaseCut(min,max,kHbtVz){}
+ virtual ~AliHBTVzCut(){}
+ protected:
+ Double_t GetValue(AliHBTParticle * p){return p->Vz();} //retruns value of the vertex
+ public:
+ ClassDef(AliHBTVzCut,1)
+
+ };
+
+
+
+
+
+#endif
--- /dev/null
+#include "AliHBTQResolutionFctns.h"
+
+AliHBTQOutResolVSQInvFctn* xxx = new AliHBTQOutResolVSQInvFctn();
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTQOutResolVSQInvFctn )
+AliHBTQOutResolVSQInvFctn::
+AliHBTQOutResolVSQInvFctn(Int_t nXbins, Double_t maxXval, Double_t minXval,
+ Int_t nYbins, Double_t maxYval, Double_t minYval):
+ AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QOutResolVSQInv","Q_{Out} Resolution vs. Q_{Inv}");
+}
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTQSideResolVSQInvFctn )
+
+AliHBTQSideResolVSQInvFctn::
+AliHBTQSideResolVSQInvFctn(Int_t nXbins, Double_t maxXval, Double_t minXval,
+ Int_t nYbins, Double_t maxYval, Double_t minYval):
+ AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QSideResolVSQInv","Q_{Side} Resolution vs. Q_{Inv}");
+}
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTQLongResolVSQInvFctn )
+
+AliHBTQLongResolVSQInvFctn::
+AliHBTQLongResolVSQInvFctn(Int_t nXbins, Double_t maxXval, Double_t minXval,
+ Int_t nYbins, Double_t maxYval, Double_t minYval):
+ AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QLongResolVSQInv","Q_{Long} Resolution vs. Q_{Inv}");
+}
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTQInvResolVSKtFctn )
+
+AliHBTQInvResolVSKtFctn::
+AliHBTQInvResolVSKtFctn(Int_t nXbins, Double_t maxXval, Double_t minXval,
+ Int_t nYbins, Double_t maxYval, Double_t minYval):
+ AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QInvResolVSKt","Q_{Inv} Resolution vs. K_{t}");
+}
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+ClassImp( AliHBTQOutResolVSKtFctn )
+
+AliHBTQOutResolVSKtFctn::
+AliHBTQOutResolVSKtFctn(Int_t nXbins, Double_t maxXval, Double_t minXval,
+ Int_t nYbins, Double_t maxYval, Double_t minYval):
+ AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QOutResolVSKt","Q_{Out} Resolution vs. K_{t} ");
+}
+
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+ClassImp( AliHBTQSideResolVSKtFctn )
+
+AliHBTQSideResolVSKtFctn::
+AliHBTQSideResolVSKtFctn(Int_t nXbins, Double_t maxXval, Double_t minXval,
+ Int_t nYbins, Double_t maxYval, Double_t minYval):
+ AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QSideResolVSKt","Q_{Side} Resolution vs. K_{t} ");
+}
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+ClassImp( AliHBTQLongResolVSKtFctn )
+
+AliHBTQLongResolVSKtFctn::
+AliHBTQLongResolVSKtFctn(Int_t nXbins, Double_t maxXval, Double_t minXval,
+ Int_t nYbins, Double_t maxYval, Double_t minYval):
+ AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QLongResolVSKt","Q_{Long} Resolution vs. K_{t} ");
+}
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+ClassImp( AliHBTQOutResolVSQOutFctn)
+
+AliHBTQOutResolVSQOutFctn::
+AliHBTQOutResolVSQOutFctn(Int_t nXbins, Double_t maxXval, Double_t minXval,
+ Int_t nYbins, Double_t maxYval, Double_t minYval):
+ AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QOutResolVSQOut","Q_{Out} Resolution vs. Q_{Out} ");
+}
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+ClassImp( AliHBTQSideResolVSQSideFctn )
+
+AliHBTQSideResolVSQSideFctn::
+AliHBTQSideResolVSQSideFctn(Int_t nXbins, Double_t maxXval, Double_t minXval,
+ Int_t nYbins, Double_t maxYval, Double_t minYval):
+ AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QSideResolVSQSide","Q_{Side} Resolution vs. Q_{Side} ");
+}
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+ClassImp( AliHBTQLongResolVSQLongFctn )
+
+AliHBTQLongResolVSQLongFctn::
+AliHBTQLongResolVSQLongFctn(Int_t nXbins, Double_t maxXval, Double_t minXval,
+ Int_t nYbins, Double_t maxYval, Double_t minYval):
+ AliHBTFourPartFctn2D(nXbins,maxXval,minXval,nYbins,maxYval,minYval)
+{
+ Rename("QLongResolVSQLong","Q_{Long} Resolution vs. Q_{Long} ");
+}
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+
+
+
--- /dev/null
+#ifndef ALIHBTQOUTVSQINVRESOLFCTN_H
+#define ALIHBTQOUTVSQINVRESOLFCTN_H
+//General Remark:
+//CMSLC means
+//Center of Mass System Longitudially Co-moving
+
+
+//this class creates resolution function of Qout
+//(difference of simulated pair Qout and recontructed pair)
+//in function of QInv
+//it inherits from AliHBTFourPartFctn2D
+// it needs two pairs to compare
+// and is two dimentional: numerator and denominator are TH2D
+
+class AliHBTQOutResolVSQInvFctn; //QOutCMSLC Res VS QInvCMSLC
+class AliHBTQSideResolVSQInvFctn; //QSideCMSLC Res VS QInvCMSLC
+class AliHBTQLongResolVSQInvFctn; //QLongCMSLC Res VS QInvCMSLC
+
+class AliHBTQInvResolVSKtFctn; //QInvCMSLC Res VS Kt
+class AliHBTQOutResolVSKtFctn; //QOutCMSLC Res VS Kt
+class AliHBTQSideResolVSKtFctn; //QSideCMSLC Res VS Kt
+class AliHBTQLongResolVSKtFctn; //QLongCMSLC Res VS Kt
+
+class AliHBTQOutResolVSQOutFctn; //QOutCMSLC Res VS QOut
+class AliHBTQSideResolVSQSideFctn;//QSideCMSLC Res VS QSide
+class AliHBTQLongResolVSQLongFctn;//QLongCMSLC Res VS QLong
+
+
+#include "AliHBTFunction.h"
+/***********************************************************************/
+/***********************************************************************/
+class AliHBTQOutResolVSQInvFctn: public AliHBTFourPartFctn2D
+ {
+ public:
+ AliHBTQOutResolVSQInvFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0,
+ Int_t nYbins = 500, Double_t maxYval = .15, Double_t minYval =-0.15);
+
+ virtual ~AliHBTQOutResolVSQInvFctn(){}
+
+ TH1* GetResult(){return fNumerator;}
+ void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+ {
+ x = partpair->GetQOutCMSLC() - trackpair->GetQOutCMSLC();
+ y = partpair->GetQInv();
+ }
+ protected:
+ private:
+ public:
+ ClassDef(AliHBTQOutResolVSQInvFctn,1)
+ };
+
+/***********************************************************************/
+/***********************************************************************/
+class AliHBTQSideResolVSQInvFctn: public AliHBTFourPartFctn2D
+ {
+ public:
+ AliHBTQSideResolVSQInvFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0,
+ Int_t nYbins = 500, Double_t maxYval = .05, Double_t minYval =-0.05);
+ virtual ~AliHBTQSideResolVSQInvFctn(){}
+
+ void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+ {
+ x = partpair->GetQSideCMSLC() - trackpair->GetQSideCMSLC();
+ y = partpair->GetQInv();
+ }
+ TH1* GetResult(){return fNumerator;}
+ protected:
+ private:
+ public:
+ ClassDef(AliHBTQSideResolVSQInvFctn,1)
+ };
+
+/***********************************************************************/
+/***********************************************************************/
+class AliHBTQLongResolVSQInvFctn: public AliHBTFourPartFctn2D
+ {
+ public:
+ AliHBTQLongResolVSQInvFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0,
+ Int_t nYbins = 500, Double_t maxYval = .05, Double_t minYval =-0.05);
+ virtual ~AliHBTQLongResolVSQInvFctn(){}
+
+ void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+ {
+ x = partpair->GetQLongCMSLC() - trackpair->GetQLongCMSLC();
+ y = partpair->GetQInv();
+ }
+ TH1* GetResult(){return fNumerator;}
+ protected:
+ private:
+ public:
+ ClassDef(AliHBTQLongResolVSQInvFctn,1)
+ };
+
+/***********************************************************************/
+/***********************************************************************/
+class AliHBTQInvResolVSKtFctn: public AliHBTFourPartFctn2D
+ {
+ public:
+ AliHBTQInvResolVSKtFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0,
+ Int_t nYbins = 500, Double_t maxYval = .05, Double_t minYval =-0.05);
+ virtual ~AliHBTQInvResolVSKtFctn(){};
+
+ void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+ {
+ x = partpair->GetQInv() - trackpair->GetQInv();
+ y = partpair->GetKt();
+ }
+ TH1* GetResult(){return fNumerator;}
+ protected:
+ private:
+ public:
+ ClassDef(AliHBTQInvResolVSKtFctn,1)
+ };
+/***********************************************************************/
+/***********************************************************************/
+class AliHBTQOutResolVSKtFctn: public AliHBTFourPartFctn2D
+ {
+ public:
+ AliHBTQOutResolVSKtFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0,
+ Int_t nYbins = 500, Double_t maxYval = .15, Double_t minYval =-0.15);
+ virtual ~AliHBTQOutResolVSKtFctn(){}
+ TH1* GetResult(){return GetNumerator();}
+ void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+ {
+ x = partpair->GetQOutCMSLC() - trackpair->GetQOutCMSLC();
+ y = partpair->GetKt();
+ }
+ protected:
+ private:
+ public:
+ ClassDef(AliHBTQOutResolVSKtFctn,1)
+ };
+/***********************************************************************/
+/***********************************************************************/
+class AliHBTQSideResolVSKtFctn: public AliHBTFourPartFctn2D
+ {
+ public:
+ AliHBTQSideResolVSKtFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0,
+ Int_t nYbins = 500, Double_t maxYval = .05, Double_t minYval =-0.05);
+ virtual ~AliHBTQSideResolVSKtFctn(){}
+ TH1* GetResult(){return GetNumerator();}
+ void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+ {
+ x = partpair->GetQSideCMSLC() - trackpair->GetQSideCMSLC();
+ y = partpair->GetKt();
+ }
+ protected:
+ private:
+ public:
+ ClassDef(AliHBTQSideResolVSKtFctn,1)
+ };
+/***********************************************************************/
+/***********************************************************************/
+class AliHBTQLongResolVSKtFctn: public AliHBTFourPartFctn2D
+ {
+ public:
+ AliHBTQLongResolVSKtFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0,
+ Int_t nYbins = 500, Double_t maxYval = .05, Double_t minYval =-0.05);
+ virtual ~AliHBTQLongResolVSKtFctn(){}
+
+ void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+ {
+ x = partpair->GetQLongCMSLC() - trackpair->GetQLongCMSLC();
+ y = partpair->GetKt();
+ }
+ TH1* GetResult(){return fNumerator;}
+ protected:
+ private:
+ public:
+ ClassDef(AliHBTQLongResolVSKtFctn,1)
+ };
+/***********************************************************************/
+/***********************************************************************/
+class AliHBTQOutResolVSQOutFctn: public AliHBTFourPartFctn2D
+ {
+ public:
+ AliHBTQOutResolVSQOutFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0,
+ Int_t nYbins = 500, Double_t maxYval = .15, Double_t minYval =-0.15);
+ virtual ~AliHBTQOutResolVSQOutFctn(){}
+
+ void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+ {
+ y = partpair->GetQOutCMSLC();
+ x = y - trackpair->GetQOutCMSLC();
+ }
+ TH1* GetResult(){return fNumerator;}
+ protected:
+ private:
+ public:
+ ClassDef(AliHBTQOutResolVSQOutFctn,1)
+ };
+
+/***********************************************************************/
+/***********************************************************************/
+
+class AliHBTQSideResolVSQSideFctn: public AliHBTFourPartFctn2D
+ {
+ public:
+ AliHBTQSideResolVSQSideFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0,
+ Int_t nYbins = 500, Double_t maxYval = .15, Double_t minYval =-0.15);
+ virtual ~AliHBTQSideResolVSQSideFctn(){}
+
+ void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+ {
+ y = partpair->GetQSideCMSLC();
+ x = y - trackpair->GetQSideCMSLC();
+ }
+ TH1* GetResult(){return fNumerator;}
+ protected:
+ private:
+ public:
+ ClassDef(AliHBTQSideResolVSQSideFctn,1)
+ };
+
+
+/***********************************************************************/
+/***********************************************************************/
+
+class AliHBTQLongResolVSQLongFctn: public AliHBTFourPartFctn2D
+ {
+ public:
+ AliHBTQLongResolVSQLongFctn(Int_t nXbins = 200, Double_t maxXval = 1.5, Double_t minXval = 0.0,
+ Int_t nYbins = 500, Double_t maxYval = .05, Double_t minYval =-0.05);
+ virtual ~AliHBTQLongResolVSQLongFctn(){}
+
+ void GetValues(AliHBTPair* trackpair, AliHBTPair* partpair, Double_t& x, Double_t& y)
+ {
+ y = partpair->GetQLongCMSLC();
+ x = y - trackpair->GetQLongCMSLC();
+ }
+ TH1* GetResult(){return fNumerator;}
+ protected:
+ private:
+ public:
+ ClassDef(AliHBTQLongResolVSQLongFctn,1)
+ };
+
+
+
+#endif
--- /dev/null
+#include "AliHBTReader.h"
+
+#include "AliHBTParticleCut.h"
+
+
+ClassImp(AliHBTReader)
+//pure virtual
+
+/*************************************************************************************/
+
+AliHBTReader::AliHBTReader()
+{
+//constructor
+ fCuts = new TObjArray();
+}
+
+/*************************************************************************************/
+
+AliHBTReader::~AliHBTReader()
+{
+//destructor
+ fCuts->SetOwner();
+ delete fCuts;
+}
+
+/*************************************************************************************/
+
+void AliHBTReader::AddParticleCut(AliHBTParticleCut* cut)
+{
+ //sets the new cut
+
+ if (!cut) //if cut is NULL return with error
+ {
+ Error("AddParticleType","NULL pointers are not accepted any more.\nIf You want to accept all particles of this type, set an empty cut ");
+ return;
+ }
+ AliHBTParticleCut *c = (AliHBTParticleCut*)cut->Clone();
+ fCuts->Add(c);
+}
+
+/*************************************************************************************/
+
+Bool_t AliHBTReader::Pass(AliHBTParticle* p)
+ {
+ //Method examines whether particle meets all cut and particle type criteria
+
+ if(p==0x0)//of corse we not pass NULL pointers
+ {
+ Warning("Pass()","No Pasaran! We never accept NULL pointers");
+ return kTRUE;
+ }
+ //if no particle is specified, we pass all particles
+ //excluding NULL pointers, of course
+
+ for(Int_t i=0; i<fCuts->GetEntriesFast(); i++)
+ {
+ AliHBTParticleCut &cut = *((AliHBTParticleCut*)fCuts->At(i));
+ if(!cut.Pass(p)) return kFALSE; //accepted
+ }
+
+ return kTRUE;//not accepted
+
+ }
+/*************************************************************************************/
+
+Bool_t AliHBTReader::Pass(Int_t pid)
+{
+//this method checks if any of existing cuts accepts this pid particles
+//or any cuts accepts all particles
+
+ if(pid == 0)
+ return kTRUE;
+
+ for(Int_t i=0; i<fCuts->GetEntriesFast(); i++)
+ {
+ AliHBTParticleCut &cut = *((AliHBTParticleCut*)fCuts->At(i));
+ //if some of cuts accepts all particles or some accepts particles of this type, accept
+ if ( (cut.GetPID() == 0) || (cut.GetPID() == pid) ) return kFALSE;
+ }
+ return kTRUE;
+}
+/*************************************************************************************/
--- /dev/null
+#ifndef ALIHBTREADER_H
+#define ALIHBTREADER_H
+
+#include <TObject.h>
+
+//Reader Base class (reads particles and tracks and
+//puts it to the AliHBTRun object
+//Piotr.Skowronski@cern.ch
+
+class AliHBTRun;
+class AliHBTEvent;
+class AliHBTParticleCut;
+class TObjArray;
+class AliHBTParticle;
+
+class AliHBTReader: public TObject
+
+{
+ public:
+ AliHBTReader();
+ virtual ~AliHBTReader();
+ //in the future this class is will read global tracking
+ virtual Int_t Read(AliHBTRun* particles, AliHBTRun *tracks) = 0;
+
+ virtual AliHBTEvent* GetParticleEvent(Int_t) = 0;
+ virtual AliHBTEvent* GetTrackEvent(Int_t) = 0;
+ virtual Int_t GetNumberOfPartEvents() = 0;
+ virtual Int_t GetNumberOfTrackEvents() = 0;
+
+
+ void AddParticleCut(AliHBTParticleCut* cut);
+
+ protected:
+
+ TObjArray *fCuts;//
+
+ Bool_t Pass(AliHBTParticle*);
+ Bool_t Pass(Int_t pid);
+
+ private:
+
+ public:
+ ClassDef(AliHBTReader,1)
+};
+
+#endif
--- /dev/null
+
+#include "AliHBTReader.h"
+
+ClassImp(AliHBTReaderITSv1)
+
+AliHBTReaderITSv1::AliHBTReaderITSv1(const Char_t* goodtracksfilename):
+ fGoodITSTracksFileName(goodtracksfilename)
+
+ {
+ fParticles = new AliHBTRun();
+ fTracks = new AliHBTRun();
+ fIsRead = kFALSE;
+ }
+AliHBTReaderITSv1::AliHBTReaderITSv1()
+{
+ delete fParticles;
+ delete fTracks;
+}
+AliHBTEvent* AliHBTReaderITSv1::GetParticleEvent(Int_t n)
+ {
+ //returns Nth event with simulated particles
+ if (!fIsRead) Read(fParticles,fTracks);
+ return fParticles->GetEvent(n);
+ }
+/********************************************************************/
+AliHBTEvent* AliHBTReaderITSv1::GetTrackEvent(Int_t n)
+ {
+ //returns Nth event with reconstructed tracks
+ if (!fIsRead) Read(fParticles,fTracks);
+ return fTracks->GetEvent(n);
+ }
+/********************************************************************/
+
+Int_t AliHBTReaderITSv1::GetNumberOfPartEvents()
+ {
+ //returns number of events of particles
+ if (!fIsRead) Read(fParticles,fTracks);
+ return fParticles->GetNumberOfEvents();
+ }
+
+/********************************************************************/
+Int_t AliHBTReaderITSv1::GetNumberOfTrackEvents()
+ {
+ //returns number of events of tracks
+ if (!fIsRead) Read(fParticles,fTracks);
+ return fTracks->GetNumberOfEvents();
+ }
+/********************************************************************/
+
+Int_t AliHBTReaderITSv1::Read(AliHBTRun* particles, AliHBTRun *tracks)
+ {
+ AliGoodTracksITSv1 *goodITStracks = new AliGoodTracksITSv1(fGoodITSTracksFileName);
+
+ if (!goodITStracks)
+ {
+ Error("AliHBTReaderITSv1::Read","Exiting due to problems with opening files. Errorcode %d",i);
+ return 1;
+ }
+ for (Int_t currentEvent = 0; currentEvent < goodITStracks->fNevents; currentEvent++)
+ {
+ for(Int_t i =0; i<goodITStracks->fGoodInEvent[currentEvent]; i++)
+ {
+ const struct GoodTrackITSv1 & gt = goodTPCTracks->GetTrack(currentEvent,i);
+
+ if(Pass(gt.code)) continue;
+
+ Double_t mass = TDatabasePDG::Instance()->GetParticle(gt.code)->Mass();
+ Double_t pEtot = TMath::Sqrt(gt.px*gt.px + gt.py*gt.py + gt.pz*gt.pz + mass*mass);
+
+ AliHBTParticle* part = new AliHBTParticle(gt.code, gt.px, gt.py, gt.pz, pEtot, gt.x, gt.y, gt.z, 0.0);
+ if(Pass(part)) { delete part; continue;}
+
+ Double_t tEtot = TMath::Sqrt( gt.pxg*gt.pxg + gt.pyg*gt.pyg + gt.pzg*gt.pzg + mass*mass);
+
+ AliHBTParticle* track = new AliHBTParticle(gt.code, gt.pxg, gt.pyg , gt.pzg, tEtot, 0., 0., 0., 0.);
+ if(Pass(track)) { delete track;continue;}
+
+ particles->AddParticle(currentEvent,part);
+ tracks->AddParticle(currentEvent,track);
+ }
+ }
+ delete goodITStracks;
+ fIsRead = kTRUE;
+ return 0;
+
+ }
+/********************************************************************/
+
+/********************************************************************/
+/********************************************************************/
+/********************************************************************/
+
+
+AliGoodTracksITSv1::~AliGoodTracksITSv1()
+{
+ delete [] fGoodInEvent;
+ for (Int_t i = 0;i<fNevents;i++)
+ delete [] fData[i];
+ delete [] fData;
+}
+/********************************************************************/
+AliGoodTracksITSv1::AliGoodTracksITSv1(const TString& infilename)
+{
+
+ cout<<"AliGoodTracksITSv1::AliGoodTracksITSv1() ....\n";
+
+
+ fNevents = 0;
+ Int_t maxevents = 500;
+
+ ifstream in(infilename.Data());
+
+ if(!in)
+ {
+ cerr<<"Can not open file with Good ITSv1 Tracks named:"<<infilename.Data()<<endl;
+ delete this;
+ return;
+ }
+
+
+ fGoodInEvent = new Int_t[maxevents];
+ fData = new struct GoodTrack* [maxevents];
+
+ Int_t i;
+ Int_t lastevent = 0;
+ fGoodInEvent[0] =0;
+ fData[0] = new struct GoodTrack[50000];
+
+ Int_t evno;
+ while(in>>evno)
+ {
+ if(lastevent>evno)
+ {
+ for(i = lastevent;i<=evno;i++)
+ {
+ fGoodInEvent[i] =0;
+ fData[i] = new struct GoodTrack[50000];
+ }
+ lastevent = evno;
+ }
+ if(fGoodInEvent[evno]>=50000)
+ {
+ cerr<<"AliGoodTracksITSv1::AliGoodTracksITSv1() : Not enough place in the array\n";
+ continue;
+ }
+ in>>fData[evno][fGoodInEvent[evno]].lab;
+ in>>fData[evno][fGoodInEvent[evno]].code;
+ in>>fData[evno][fGoodInEvent[evno]].px;
+ in>>fData[evno][fGoodInEvent[evno]].py;
+ in>>fData[evno][fGoodInEvent[evno]].pz;
+ in>>fData[evno][fGoodInEvent[evno]].x;
+ in>>fData[evno][fGoodInEvent[evno]].y;
+ in>>fData[evno][fGoodInEvent[evno]].z;
+
+ /* cout<<evno<<" ";
+ cout<<fData[evno][fGoodInEvent[evno]].lab;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].code;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].px;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].py;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].pz;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].x;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].y;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].z;
+ cout<<"\n";
+ */
+ fGoodInEvent[evno]++;
+ }
+ fNevents = evno+1;
+ in.close();
+ cout<<"AliGoodTracksITSv1::AliGoodTracksITSv1() .... Done\n";
+}
+
+
+
+const GoodTrack& AliGoodTracksITSv1::GetTrack(Int_t event, Int_t n) const
+ {
+
+ if( (event>fNevents) || (event<0))
+ {
+ gROOT->Fatal("AliGoodTracksITSv1::GetTrack","No such Event %d",event);
+ }
+ if( (n>fGoodInEvent[event]) || (n<0))
+ {
+ gROOT->Fatal("AliGoodTracksITSv1::GetTrack","No such Good TPC Track %d",n);
+ }
+ return fData[event][n];
+
+ }
--- /dev/null
+#ifndef ALIHBTREADERITSV1_H
+#define ALIHBTREADERITSV1_H
+
+#include "AliHBTReader.h"
+
+#include <TString.h>
+
+
+class AliHBTReaderITSv1: public AliHBTReader
+{
+ public:
+ AliHBTReaderITS(const Char_t* goodtracksfilename = "itsgood_tracks");
+ virtual ~AliHBTReaderITS();
+
+ Int_t Read(AliHBTRun* particles, AliHBTRun *tracks);//reads tracks and particles and puts them in runs
+
+ AliHBTEvent* GetParticleEvent(Int_t);//returns pointer to event with particles
+ AliHBTEvent* GetTrackEvent(Int_t);//returns pointer to event with particles
+ Int_t GetNumberOfPartEvents();//returns number of particle events
+ Int_t GetNumberOfTrackEvents();//returns number of track events
+
+ protected:
+ TString fGoodITSTracksFileName;
+
+ AliHBTRun* fParticles; //!simulated particles
+ AliHBTRun* fTracks; //!reconstructed tracks (particles)
+
+ Bool_t fIsRead;//flag indicating if the data are already read
+
+ private:
+ public:
+ ClassDef(AliHBTReaderITS,1)
+};
+
+struct GoodTrackITSv1 //good tracks produced by ITSComparison V1
+{
+ Int_t fEventN; //event number
+ Int_t lab;
+ Int_t code;
+ Float_t px,py,pz;
+ Float_t x,y,z;
+ Float_t pxg,pyg,pzg,ptg;
+ Bool_t flag;
+};
+
+class AliGoodTracksITSv1
+ {
+ //container for good tracks ITS tracking V1
+ //this class is for internal use only
+ friend class AliHBTReaderITSv1;
+
+ private:
+ AliGoodTracksITSv1(const TString& infilename = TString("itsgood_tracks"));
+ ~AliGoodTracksITSv1();
+
+ const GoodTrackITSv1& GetTrack(Int_t event, Int_t n) const;
+
+ Int_t fNevents; //Number of events
+ Int_t* fGoodInEvent; //Numbers of good track in event
+ struct GoodTrack **fData;
+ };
+
--- /dev/null
+
+#include "AliHBTReaderITSv2.h"
+
+ClassImp(AliHBTReaderITSv2)
+
--- /dev/null
+#ifndef ALIHBTREADERITSV2_H
+#define ALIHBTREADERITSV2_H
+
+#include "AliHBTReader.h"
+
+#include <TString.h>
+
+
+class AliHBTReaderITSv2: public AliHBTReaderTPC
+{
+ public:
+ AliHBTReaderITSv2(const Char_t* trackfilename = "AliITStracksV2.root",
+ const Char_t* clusterfilename = "AliITSclustersV2.root",
+ const Char_t* goodtracksfilename = "good_tracks_its",
+ const Char_t* galicefilename = "");
+
+ virtual ~AliHBTReaderITSv2();
+
+ Int_t Read(AliHBTRun* particles, AliHBTRun *tracks);//reads tracks and particles and puts them in runs
+
+ protected:
+ private:
+ public:
+ ClassDef(AliHBTReaderITSv2,1)
+};
--- /dev/null
+#include "AliHBTReaderPPprod.h"
+
+#include <iostream.h>
+#include <fstream.h>
+#include <TString.h>
+#include <TTree.h>
+#include <TFile.h>
+
+#include <AliTPCtrack.h>
+#include <AliTPCParam.h>
+#include <AliTPCtracker.h>
+
+#include "AliRun.h"
+#include "AliHBTRun.h"
+#include "AliHBTEvent.h"
+#include "AliHBTParticle.h"
+#include "AliHBTParticleCut.h"
+
+
+ClassImp(AliHBTReaderPPprod)
+
+AliHBTReaderPPprod::
+ AliHBTReaderPPprod(const Char_t* trackfilename,const Char_t* clusterfilename,
+ const Char_t* goodtracksfilename,const Char_t* galicefilename):
+ fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
+ fGAliceFileName(galicefilename),
+ fGoodTPCTracksFileName(goodtracksfilename)
+{
+ //constructor, only file names are set
+ //Defaults:
+ // trackfilename = "AliTPCtracks.root"
+ // clusterfilename = "AliTPCclusters.root"
+ // goodtracksfilename = "good_tracks_tpc"
+ // galicefilename = "" - this means: Do not open gAlice file -
+ // just leave the global pointer untached
+
+ fParticles = new AliHBTRun();
+ fTracks = new AliHBTRun();
+
+ fTracksFile = 0x0; //files are opened during reading only
+ fClustersFile = 0x0;
+
+ fIsRead = kFALSE;
+}
+/********************************************************************/
+
+AliHBTReaderPPprod::~AliHBTReaderPPprod()
+ {
+ delete fParticles;
+ delete fTracks;
+ }
+/********************************************************************/
+
+AliHBTEvent* AliHBTReaderPPprod::GetParticleEvent(Int_t n)
+ {
+ //returns Nth event with simulated particles
+ if (!fIsRead) Read(fParticles,fTracks);
+ return fParticles->GetEvent(n);
+ }
+/********************************************************************/
+AliHBTEvent* AliHBTReaderPPprod::GetTrackEvent(Int_t n)
+ {
+ //returns Nth event with reconstructed tracks
+ if (!fIsRead) Read(fParticles,fTracks);
+ return fTracks->GetEvent(n);
+ }
+/********************************************************************/
+
+Int_t AliHBTReaderPPprod::GetNumberOfPartEvents()
+ {
+ //returns number of events of particles
+ if (!fIsRead) Read(fParticles,fTracks);
+ return fParticles->GetNumberOfEvents();
+ }
+
+/********************************************************************/
+Int_t AliHBTReaderPPprod::GetNumberOfTrackEvents()
+ {
+ //returns number of events of tracks
+ if (!fIsRead) Read(fParticles,fTracks);
+ return fTracks->GetNumberOfEvents();
+ }
+/********************************************************************/
+
+
+Int_t AliHBTReaderPPprod::Read(AliHBTRun* particles, AliHBTRun *tracks)
+ {
+ //reads data and puts put to the particles and tracks objects
+ //reurns 0 if everything is OK
+ //
+ Int_t i; //iterator and some temprary values
+ Int_t Nevents;
+ if (!particles) //check if an object is instatiated
+ {
+ Error("AliHBTReaderPPprod::Read"," particles object must instatiated before passing it to the reader");
+ }
+ if (!tracks) //check if an object is instatiated
+ {
+ Error("AliHBTReaderPPprod::Read"," tracks object must instatiated before passing it to the reader");
+ }
+ particles->Reset();//clear runs == delete all old events
+ tracks->Reset();
+
+ if( (i=OpenFiles()) )
+ {
+ Error("AliHBTReaderPPprod::Read","Exiting due to problems with opening files. Errorcode %d",i);
+ return i;
+ }
+
+ AliGoodTracksPP *goodTPCTracks = new AliGoodTracksPP(fGoodTPCTracksFileName);
+ if (!goodTPCTracks)
+ {
+ Error("AliHBTReaderPPprod::Read","Exiting due to problems with opening files. Errorcode %d",i);
+ return 1;
+ }
+
+ Nevents = 100;
+
+ fClustersFile->cd();
+ AliTPCParam *TPCParam= (AliTPCParam*)fClustersFile->Get("75x40_100x60");
+ if (!TPCParam)
+ {
+ Error("AliHBTReaderPPprod::Read","TPC parameters have not been found !\n");
+ return 1;
+ }
+
+ TObjArray *tarray = new TObjArray(5000);
+ tarray->SetOwner();// this causes memory leak, but in some cases deleting is infinite loop
+
+
+ for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)
+ {
+ cout<<"Reading Event "<<currentEvent<<endl;
+ /**************************************/
+ /**************************************/
+ /**************************************/
+ fTracksFile->cd();
+ Char_t treename[100];
+ sprintf(treename,"TreeT_TPC_%d",currentEvent);
+
+ TTree *tracktree=0;
+
+ tracktree=(TTree*)fTracksFile->Get(treename);
+ if (!tracktree)
+ {
+ Error("AliHBTReaderPPprod::Read","Can't get a tree with TPC tracks !\n");
+ return 1;
+ }
+
+ TBranch *trackbranch=tracktree->GetBranch("tracks");
+ if (!trackbranch)
+ {
+ Error("AliHBTReaderPPprod::Read","Can't get a branch with TPC tracks !\n");
+ return 2;
+ }
+ Int_t NTPCtracks=(Int_t)tracktree->GetEntries();
+ cout<<"Found "<<NTPCtracks<<" TPC tracks.\n";
+ //Copy tracks to array
+
+ AliTPCtrack *iotrack=0;
+
+ fClustersFile->cd();
+ AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent);
+ if (!tracker)
+ {
+ Error("AliHBTReaderPPprod::Read","Can't get a tracker !\n");
+ return 3;
+ }
+ tracker->LoadInnerSectors();
+ tracker->LoadOuterSectors();
+
+ for (i=0; i<NTPCtracks; i++)
+ {
+ iotrack=new AliTPCtrack;
+ trackbranch->SetAddress(&iotrack);
+ tracktree->GetEvent(i);
+ tracker->CookLabel(iotrack,0.1);
+ tarray->AddLast(iotrack);
+ }
+
+
+ fTracksFile->Delete(treename);//delete tree from memmory (and leave untached on disk)- we do not need it any more
+ fTracksFile->Delete("tracks");
+
+ delete tracker;
+
+ tracker = 0x0;
+ trackbranch = 0x0;
+ tracktree = 0x0;
+
+ Int_t & ngood = goodTPCTracks->fGoodInEvent[currentEvent];
+
+ Double_t xk;
+ Double_t par[5];
+ Float_t phi, lam, pt;
+ Int_t label;
+ Bool_t found;
+
+ for (i=0; i<ngood; i++)
+ {
+ const struct GoodTrack & gt = goodTPCTracks->GetTrack(currentEvent,i);
+
+ if(Pass(gt.code)) continue;
+
+ label = gt.lab;
+ found = kFALSE; //guard in case we don't find track with such a label
+ for (Int_t j=0;j<NTPCtracks;j++)
+ {
+ iotrack = (AliTPCtrack*)tarray->At(j);
+ if (iotrack->GetLabel() == label)
+ {
+ found = kTRUE;
+ break;
+ }
+ }
+ if(!found)
+ {
+ Warning("Read",
+ "Sth is going wrong with tracks - there is no TPC track corresponding to goodtrack.\nGood tack label %d",label);
+ continue; //put comunicate on the screen and continue loop
+ }
+
+ Double_t mass = TDatabasePDG::Instance()->GetParticle(gt.code)->Mass();
+ Double_t pEtot = TMath::Sqrt(gt.px*gt.px + gt.py*gt.py + gt.pz*gt.pz + mass*mass);
+
+ AliHBTParticle* part = new AliHBTParticle(gt.code, gt.px, gt.py, gt.pz, pEtot, gt.x, gt.y, gt.z, 0.0);
+ if(Pass(part)) continue;
+
+
+ iotrack->PropagateTo(gt.x);
+ iotrack->GetExternalParameters(xk,par);
+ phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
+ if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
+ if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
+ lam=par[3];
+ pt=1.0/TMath::Abs(par[4]);
+
+ Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
+ Double_t tpy = pt * TMath::Sin(phi); //track x coordinate of momentum
+ Double_t tpz = pt * lam;
+
+ Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);
+
+ AliHBTParticle* track = new AliHBTParticle(gt.code, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
+ if(Pass(track)) continue;
+
+ particles->AddParticle(currentEvent,part);
+ tracks->AddParticle(currentEvent,track);
+
+ }
+ tarray->Clear();
+ /**************************************/
+ /**************************************/
+ /**************************************/
+ }
+
+ CloseFiles();
+ delete tarray;
+ delete goodTPCTracks;
+ fIsRead = kTRUE;
+ return 0;
+ }
+
+/********************************************************************/
+Int_t AliHBTReaderPPprod::OpenFiles()
+{
+
+ fTracksFile = 0;
+ fTracksFile=TFile::Open("AliTPCtracks.root");
+ if (!fTracksFile->IsOpen())
+ {
+ Error("AliHBTReaderPPprod::OpenFiles","Can't open AliTPCtracks.root");
+ return 1;
+ }
+
+ fClustersFile = 0;
+
+ fClustersFile=TFile::Open("AliTPCclusters.root");
+ if (!fClustersFile->IsOpen())
+ {
+ Error("AliHBTReaderPPprod::OpenFiles","Can't open AliTPCclusters.root");
+ return 2;
+ }
+
+
+
+ return 0;
+}
+
+
+
+/********************************************************************/
+
+void AliHBTReaderPPprod::CloseFiles()
+{
+ fTracksFile->Close();
+ fClustersFile->Close();
+}
+
+/********************************************************************/
+
+/********************************************************************/
+/********************************************************************/
+/********************************************************************/
+
+
+AliGoodTracksPP::~AliGoodTracksPP()
+{
+ delete [] fGoodInEvent;
+ for (Int_t i = 0;i<fNevents;i++)
+ delete [] fData[i];
+ delete [] fData;
+}
+/********************************************************************/
+AliGoodTracksPP::AliGoodTracksPP(const TString& infilename)
+{
+
+ fNevents = 100;
+ ifstream in(infilename.Data());
+
+ if(!in)
+ {
+ cerr<<"Can not open file with Good TPC Tracks named:"<<infilename.Data()<<endl;
+ delete this;
+ return;
+ }
+
+
+ fGoodInEvent = new Int_t[fNevents];
+ fData = new struct GoodTrack* [fNevents];
+
+ Int_t i;
+ for( i = 0;i<fNevents;i++)
+ {
+ fGoodInEvent[i] =0;
+ fData[i] = new struct GoodTrack[500];
+ }
+ Float_t tmp;
+ Int_t evno;
+ while(in>>evno)
+ {
+ if(fGoodInEvent[evno]>=500)
+ {
+ cerr<<"AliGoodTracksPP::AliGoodTracksPP() : Not enough place in the array\n";
+ continue;
+ }
+ in>>fData[evno][fGoodInEvent[evno]].lab;
+ in>>fData[evno][fGoodInEvent[evno]].code;
+ in>>fData[evno][fGoodInEvent[evno]].px;
+ in>>fData[evno][fGoodInEvent[evno]].py;
+ in>>fData[evno][fGoodInEvent[evno]].pz;
+ in>>fData[evno][fGoodInEvent[evno]].x;
+ in>>fData[evno][fGoodInEvent[evno]].y;
+ in>>fData[evno][fGoodInEvent[evno]].z;
+
+ in>>tmp;
+ in>>tmp;
+ in>>tmp;
+
+ /* cout<<evno<<" ";
+ cout<<fData[evno][fGoodInEvent[evno]].lab;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].code;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].px;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].py;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].pz;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].x;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].y;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].z;
+ cout<<"\n";
+ */
+ fGoodInEvent[evno]++;
+ }
+ in.close();
+ cout<<"AliGoodTracksPP::AliGoodTracksPP() .... Done\n";
+}
+
+
+
+const GoodTrack& AliGoodTracksPP::GetTrack(Int_t event, Int_t n) const
+ {
+
+ if( (event>fNevents) || (event<0))
+ {
+ gAlice->Fatal("AliGoodTracksPP::GetTrack","No such Event %d",event);
+ }
+ if( (n>fGoodInEvent[event]) || (n<0))
+ {
+ gAlice->Fatal("AliGoodTracksPP::GetTrack","No such Good TPC Track %d",n);
+ }
+ return fData[event][n];
+
+ }
--- /dev/null
+#ifndef ALIHBTREADERPPPROD_H
+#define ALIHBTREADERPPPROD_H
+
+#include "AliHBTReader.h"
+#include "AliHBTReaderTPC.h"
+
+//This reader reads tracks AliTPCtracks.root
+// particles form tpc_good_tracks
+//I am aware that this file is temporary however we do not have any other PID yet
+//Piotr.Skowronski@cern.ch
+
+#include <TString.h>
+class TFile;
+
+class AliHBTReaderPPprod: public AliHBTReader
+{
+ public:
+ AliHBTReaderPPprod(const Char_t* trackfilename = "AliTPCtracks.root",
+ const Char_t* clusterfilename = "AliTPCclusters.root",
+ const Char_t* goodtracksfilename = "good_tracks_tpc",
+ const Char_t* galicefilename = "");
+
+ virtual ~AliHBTReaderPPprod();
+
+ Int_t Read(AliHBTRun* particles, AliHBTRun *tracks); //reads tracks and particles and puts them in runs
+
+ AliHBTEvent* GetParticleEvent(Int_t); //returns pointer to event with particles
+ AliHBTEvent* GetTrackEvent(Int_t);//returns pointer to event with particles
+ Int_t GetNumberOfPartEvents(); //returns number of particle events
+ Int_t GetNumberOfTrackEvents();//returns number of track events
+
+ protected:
+ //in the future this class is will read global tracking
+
+
+ Int_t OpenFiles(); //opens files to be read
+ void CloseFiles(); //close files
+
+ AliHBTRun* fParticles; //!simulated particles
+ AliHBTRun* fTracks; //!reconstructed tracks (particles)
+
+ TString fTrackFileName; //name of the file with tracks
+ TString fClusterFileName;//name of the file with clusters
+ TString fGAliceFileName;//name of the file with galice.root
+ TString fGoodTPCTracksFileName; //name of text file with good tracks
+
+ TFile *fTracksFile; //file with tracks
+ TFile *fClustersFile;//file with clusters
+
+ Bool_t fIsRead; //flag indicating if the data are already read
+
+
+ private:
+ public:
+ ClassDef(AliHBTReaderPPprod,1)
+};
+
+
+class AliGoodTracksPP
+ {
+ //container for good tracks
+ //this class is for internal use only
+
+ friend class AliHBTReaderPPprod;
+
+ private:
+ AliGoodTracksPP(const TString& infilename = TString("good_tracks_tpc")); //constructor
+ ~AliGoodTracksPP(); //dctor
+
+ const GoodTrack& GetTrack(Int_t event, Int_t n) const; //returns reference to the nth good track in event "event"
+
+ Int_t fNevents; //Number of events
+ Int_t* fGoodInEvent; //Numbers of good track in event
+ struct GoodTrack **fData;
+ };
+
+
+#endif
--- /dev/null
+#include "AliHBTReaderTPC.h"
+
+#include <iostream.h>
+#include <fstream.h>
+#include <TString.h>
+#include <TTree.h>
+#include <TFile.h>
+
+#include <AliTPCtrack.h>
+#include <AliTPCParam.h>
+#include <AliTPCtracker.h>
+
+#include "AliRun.h"
+#include "AliHBTRun.h"
+#include "AliHBTEvent.h"
+#include "AliHBTParticle.h"
+#include "AliHBTParticleCut.h"
+
+
+ClassImp(AliHBTReaderTPC)
+//reader for TPC tracking
+//needs galice.root, AliTPCtracks.root, AliTPCclusters.root, good_tracks_tpc
+//
+//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
+//Piotr.Skowronski@cern.ch
+
+AliHBTReaderTPC::
+ AliHBTReaderTPC(const Char_t* trackfilename,const Char_t* clusterfilename,
+ const Char_t* goodtracksfilename,const Char_t* galicefilename):
+ fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
+ fGAliceFileName(galicefilename),
+ fGoodTPCTracksFileName(goodtracksfilename)
+{
+ //constructor, only file names are set
+ //Defaults:
+ // trackfilename = "AliTPCtracks.root"
+ // clusterfilename = "AliTPCclusters.root"
+ // goodtracksfilename = "good_tracks_tpc"
+ // galicefilename = "" - this means: Do not open gAlice file -
+ // just leave the global pointer untached
+
+ fParticles = new AliHBTRun();
+ fTracks = new AliHBTRun();
+
+ fTracksFile = 0x0; //files are opened during reading only
+ fClustersFile = 0x0;
+
+ fIsRead = kFALSE;
+}
+/********************************************************************/
+
+AliHBTReaderTPC::~AliHBTReaderTPC()
+ {
+ //desctructor
+ delete fParticles;
+ delete fTracks;
+ }
+/********************************************************************/
+
+AliHBTEvent* AliHBTReaderTPC::GetParticleEvent(Int_t n)
+ {
+ //returns Nth event with simulated particles
+ if (!fIsRead) Read(fParticles,fTracks);
+ return fParticles->GetEvent(n);
+ }
+/********************************************************************/
+AliHBTEvent* AliHBTReaderTPC::GetTrackEvent(Int_t n)
+ {
+ //returns Nth event with reconstructed tracks
+ if (!fIsRead) Read(fParticles,fTracks);
+ return fTracks->GetEvent(n);
+ }
+/********************************************************************/
+
+Int_t AliHBTReaderTPC::GetNumberOfPartEvents()
+ {
+ //returns number of events of particles
+ if (!fIsRead) Read(fParticles,fTracks);
+ return fParticles->GetNumberOfEvents();
+ }
+
+/********************************************************************/
+Int_t AliHBTReaderTPC::GetNumberOfTrackEvents()
+ {
+ //returns number of events of tracks
+ if (!fIsRead) Read(fParticles,fTracks);
+ return fTracks->GetNumberOfEvents();
+ }
+/********************************************************************/
+
+
+Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks)
+ {
+ //reads data and puts put to the particles and tracks objects
+ //reurns 0 if everything is OK
+ //
+ Int_t i; //iterator and some temprary values
+ Int_t Nevents;
+ if (!particles) //check if an object is instatiated
+ {
+ Error("AliHBTReaderTPC::Read"," particles object must instatiated before passing it to the reader");
+ }
+ if (!tracks) //check if an object is instatiated
+ {
+ Error("AliHBTReaderTPC::Read"," tracks object must instatiated before passing it to the reader");
+ }
+ particles->Reset();//clear runs == delete all old events
+ tracks->Reset();
+
+ if( (i=OpenFiles()) )
+ {
+ Error("AliHBTReaderTPC::Read","Exiting due to problems with opening files. Errorcode %d",i);
+ return i;
+ }
+
+ AliGoodTracks *goodTPCTracks = new AliGoodTracks(fGoodTPCTracksFileName);
+ if (!goodTPCTracks)
+ {
+ Error("AliHBTReaderTPC::Read","Exiting due to problems with opening files. Errorcode %d",i);
+ return 1;
+ }
+
+
+ if (gAlice->TreeE())//check if tree E exists
+ {
+ Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice
+ cout<<"Found "<<Nevents<<endl;
+ }
+ else
+ {//if not return an error
+ Error("AliHBTReaderPPprod::Read","Can not find Header tree (TreeE) in gAlice");
+ return 1;
+ }
+
+ fClustersFile->cd();//set cluster file active
+ AliTPCParam *TPCParam= (AliTPCParam*)fClustersFile->Get("75x40_100x60");
+ if (!TPCParam)
+ {
+ Error("AliHBTReaderTPC::Read","TPC parameters have not been found !\n");
+ return 1;
+ }
+
+ TObjArray *tarray = new TObjArray(5000); //cotainer for tpc tracks
+ tarray->SetOwner(); //set the ownership of the objects it contains
+ //when array is is deleted or cleared all objects
+ //that it contains are deleted
+
+ for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
+ {
+ cout<<"Reading Event "<<currentEvent<<endl;
+ /**************************************/
+ /**************************************/
+ /**************************************/
+ fTracksFile->cd();//set track file active
+
+ Char_t treename[100];
+ sprintf(treename,"TreeT_TPC_%d",currentEvent);//prepare name of the tree
+
+ TTree *tracktree=0;
+
+ tracktree=(TTree*)fTracksFile->Get(treename);//get the tree
+ if (!tracktree) //check if we got the tree
+ {//if not return with error
+ Error("AliHBTReaderTPC::Read","Can't get a tree with TPC tracks !\n");
+ return 1;
+ }
+
+ TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks
+ if (!trackbranch) ////check if we got the branch
+ {//if not return with error
+ Error("AliHBTReaderTPC::Read","Can't get a branch with TPC tracks !\n");
+ return 2;
+ }
+ Int_t NTPCtracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks
+ cout<<"Found "<<NTPCtracks<<" TPC tracks.\n";
+ //Copy tracks to array
+
+ AliTPCtrack *iotrack=0;
+
+ fClustersFile->cd();//set cluster file active
+ AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent);//create the tacker for this event
+ if (!tracker) //check if it has created succeffuly
+ {//if not return with error
+ Error("AliHBTReaderTPC::Read","Can't get a tracker !\n");
+ return 3;
+ }
+ tracker->LoadInnerSectors();
+ tracker->LoadOuterSectors();
+
+ for (i=0; i<NTPCtracks; i++) //loop over all tpc tracks
+ {
+ iotrack=new AliTPCtrack; //create new tracks
+ trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file)
+ tracktree->GetEvent(i); //stream track i to the iotrack
+ tracker->CookLabel(iotrack,0.1); //calculate (cook) the label of the tpc track
+ //which is the label of corresponding simulated particle
+ tarray->AddLast(iotrack); //put the track in the array
+ }
+
+ fTracksFile->Delete(treename);//delete tree from memmory (and leave untached on disk)- we do not need it any more
+ fTracksFile->Delete("tracks");//delete branch from memmory
+ delete tracker; //delete tracker
+
+ tracker = 0x0;
+ trackbranch = 0x0;
+ tracktree = 0x0;
+
+ Int_t & ngood = goodTPCTracks->fGoodInEvent[currentEvent]; //number of good tracks in the current event
+
+ Double_t xk;
+ Double_t par[5];
+ Float_t phi, lam, pt;//angles and transverse momentum
+ Int_t label; //label of the current track
+ Bool_t found; //flag indicated wether we managed to match good_tpc_track with track
+
+ for (i=0; i<ngood; i++) //loop over all good tracks
+ {
+ const struct GoodTrack & gt = goodTPCTracks->GetTrack(currentEvent,i); //get ith goog track
+
+ if(Pass(gt.code)) continue; //check if we are intersted with particles of this type
+ //if not take next partilce
+
+ label = gt.lab;
+ found = kFALSE; //guard in case we don't find track with such a label
+ for (Int_t j=0;j<NTPCtracks;j++)//lopp over all tpc tracks
+ {
+ iotrack = (AliTPCtrack*)tarray->At(j);
+ if (iotrack->GetLabel() == label) //if the label is the same
+ {
+ found = kTRUE; //we found the track
+ break;
+ }
+ }
+ if(!found) //check if we found the track
+ {
+ Warning("Read",
+ "Sth is going wrong with tracks - there is no TPC track corresponding to goodtrack.\nGood tack label %d",label);
+ continue; //put comunicate on the screen and continue loop
+ }
+
+ Double_t mass = TDatabasePDG::Instance()->GetParticle(gt.code)->Mass();//CMS mass of this particle
+ Double_t pEtot = TMath::Sqrt(gt.px*gt.px + gt.py*gt.py + gt.pz*gt.pz + mass*mass); //particle total energy
+
+ AliHBTParticle* part = new AliHBTParticle(gt.code, gt.px, gt.py, gt.pz, pEtot, gt.x, gt.y, gt.z, 0.0);
+ if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
+ //if it does not delete it and take next good track
+
+ iotrack->PropagateTo(gt.x);
+ iotrack->GetExternalParameters(xk,par); //get properties of the track
+ phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
+ if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
+ if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
+ lam=par[3];
+ pt=1.0/TMath::Abs(par[4]);
+
+ Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
+ Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
+ Double_t tpz = pt * lam; //track z coordinate of momentum
+
+ Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
+
+ AliHBTParticle* track = new AliHBTParticle(gt.code, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
+ if(Pass(track)) { delete track;continue;}//check if meets all criteria of any of our cuts
+ //if it does not delete it and take next good track
+
+ particles->AddParticle(currentEvent,part);//put track and particle on the run
+ tracks->AddParticle(currentEvent,track);
+
+ }
+ tarray->Clear(); //clear the array
+
+ /**************************************/
+ /**************************************/
+ /**************************************/
+ }
+
+ //save environment (resouces) --
+ //clean your place after the work
+ CloseFiles();
+ delete tarray;
+ delete goodTPCTracks;
+ fIsRead = kTRUE;
+ return 0;
+ }
+
+/********************************************************************/
+Int_t AliHBTReaderTPC::OpenFiles()
+{
+ //opens all the files
+ fTracksFile = 0;
+ fTracksFile=TFile::Open(fTrackFileName.Data());
+ if (!fTracksFile->IsOpen())
+ {
+ Error("AliHBTReaderTPC::OpenFiles","Can't open file with tacks named ",fTrackFileName.Data());
+ return 1;
+ }
+
+ fClustersFile = 0;
+
+ fClustersFile=TFile::Open(fClusterFileName.Data());
+ if (!fClustersFile->IsOpen())
+ {
+ Error("AliHBTReaderTPC::OpenFiles","Can't open file with TPC clusters named ",fClusterFileName.Data());
+ return 2;
+ }
+
+ return 0;
+}
+
+
+
+/********************************************************************/
+
+void AliHBTReaderTPC::CloseFiles()
+{
+ //closes the files
+ fTracksFile->Close();
+ fClustersFile->Close();
+}
+
+/********************************************************************/
+
+/********************************************************************/
+/********************************************************************/
+/********************************************************************/
+
+
+AliGoodTracks::~AliGoodTracks()
+{
+//destructor
+ delete [] fGoodInEvent;
+ for (Int_t i = 0;i<fNevents;i++)
+ delete [] fData[i];
+ delete [] fData;
+}
+/********************************************************************/
+AliGoodTracks::AliGoodTracks(const TString& infilename)
+{
+
+ cout<<"AliGoodTracks::AliGoodTracks() ....\n";
+ if(!gAlice)
+ {
+ cerr<<"There is no gAlice"<<endl;
+ delete this;
+ return;
+ }
+
+ if (!gAlice->TreeE())
+ {
+ cerr<<"Can not find Header tree (TreeE) in gAlice"<<endl;
+ delete this;
+ return;
+ }
+
+ fNevents = (Int_t)gAlice->TreeE()->GetEntries();
+ //fNevents = 100;
+ cout<<fNevents<<" FOUND"<<endl;
+ ifstream in(infilename.Data());
+
+ if(!in)
+ {
+ cerr<<"Can not open file with Good TPC Tracks named:"<<infilename.Data()<<endl;
+ delete this;
+ return;
+ }
+
+
+ fGoodInEvent = new Int_t[fNevents];
+ fData = new struct GoodTrack* [fNevents];
+
+ Int_t i;
+ for( i = 0;i<fNevents;i++)
+ {
+ fGoodInEvent[i] =0;
+ fData[i] = new struct GoodTrack[50000];
+ }
+
+ Int_t evno;
+ while(in>>evno)
+ {
+ if(fGoodInEvent[evno]>=50000)
+ {
+ cerr<<"AliGoodTracks::AliGoodTracks() : Not enough place in the array\n";
+ continue;
+ }
+ in>>fData[evno][fGoodInEvent[evno]].lab;
+ in>>fData[evno][fGoodInEvent[evno]].code;
+ in>>fData[evno][fGoodInEvent[evno]].px;
+ in>>fData[evno][fGoodInEvent[evno]].py;
+ in>>fData[evno][fGoodInEvent[evno]].pz;
+ in>>fData[evno][fGoodInEvent[evno]].x;
+ in>>fData[evno][fGoodInEvent[evno]].y;
+ in>>fData[evno][fGoodInEvent[evno]].z;
+
+ /* cout<<evno<<" ";
+ cout<<fData[evno][fGoodInEvent[evno]].lab;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].code;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].px;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].py;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].pz;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].x;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].y;
+ cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].z;
+ cout<<"\n";
+ */
+ fGoodInEvent[evno]++;
+ }
+ in.close();
+ cout<<"AliGoodTracks::AliGoodTracks() .... Done\n";
+}
+
+
+
+const GoodTrack& AliGoodTracks::GetTrack(Int_t event, Int_t n) const
+ {
+
+ if( (event>fNevents) || (event<0))
+ {
+ gAlice->Fatal("AliGoodTracks::GetTrack","No such Event %d",event);
+ }
+ if( (n>fGoodInEvent[event]) || (n<0))
+ {
+ gAlice->Fatal("AliGoodTracks::GetTrack","No such Good TPC Track %d",n);
+ }
+ return fData[event][n];
+
+ }
--- /dev/null
+#ifndef ALIHBTREADERTPC_H
+#define ALIHBTREADERTPC_H
+
+#include "AliHBTReader.h"
+
+
+//This reader reads tracks AliTPCtracks.root
+// particles form tpc_good_tracks
+//I am aware that this file is temporary however we do not have any other PID
+//Piotr.Skowronski@cern.ch
+//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
+
+#include <TString.h>
+class TFile;
+
+class AliHBTReaderTPC: public AliHBTReader
+{
+ public:
+ AliHBTReaderTPC(const Char_t* trackfilename = "AliTPCtracks.root",
+ const Char_t* clusterfilename = "AliTPCclusters.root",
+ const Char_t* goodtracksfilename = "good_tracks_tpc",
+ const Char_t* galicefilename = "");
+
+ virtual ~AliHBTReaderTPC();
+
+ Int_t Read(AliHBTRun* particles, AliHBTRun *tracks);//reads tracks and particles and puts them in runs
+
+ AliHBTEvent* GetParticleEvent(Int_t);//returns pointer to event with particles
+ AliHBTEvent* GetTrackEvent(Int_t);//returns pointer to event with particles
+ Int_t GetNumberOfPartEvents();//returns number of particle events
+ Int_t GetNumberOfTrackEvents();//returns number of track events
+
+ protected:
+ //in the future this class is will read global tracking
+
+
+ Int_t OpenFiles();//opens files to be read
+ void CloseFiles();//close files
+
+ AliHBTRun* fParticles; //!simulated particles
+ AliHBTRun* fTracks; //!reconstructed tracks (particles)
+
+ TString fTrackFileName;//name of the file with tracks
+ TString fClusterFileName;//name of the file with clusters
+ TString fGAliceFileName;//name of the file with galice.root
+ TString fGoodTPCTracksFileName;//name of text file with good tracks
+
+ TFile *fTracksFile;//file with tracks
+ TFile *fClustersFile;//file with clusters
+
+ Bool_t fIsRead;//flag indicating if the data are already read
+
+
+ private:
+ public:
+ ClassDef(AliHBTReaderTPC,1)
+};
+
+
+struct GoodTrack //data of good tracks produced by AliTPCComparison.C
+ {
+ Int_t lab;
+ Int_t code;
+ Float_t px,py,pz;
+ Float_t x,y,z;
+ };
+
+class AliGoodTracks
+ {
+ //this class is for internal use only
+ friend class AliHBTReaderTPC;
+
+ private:
+ AliGoodTracks(const TString& infilename = TString("good_tracks_tpc"));
+ ~AliGoodTracks();
+
+ const GoodTrack& GetTrack(Int_t event, Int_t n) const;
+
+ Int_t fNevents; //Number of events
+ Int_t* fGoodInEvent; //Numbers of good track in event
+ struct GoodTrack **fData;
+ };
+
+
+#endif
--- /dev/null
+#include "AliHBTRun.h"
+
+#include <TObjArray.h>
+
+
+ClassImp(AliHBTRun)
+/**************************************************************************/
+
+AliHBTRun::AliHBTRun()
+ {
+ //contructor
+
+ fEvents = new TObjArray();//create array for AliHBTEvents
+ if(!fEvents) Fatal("AliHBTRun::AliHBTRun","Can not allocate memory");
+ fEvents->SetOwner(); //array is an owner: when is deleted or cleared it deletes objects that it contains
+ }
+/**************************************************************************/
+AliHBTRun::~AliHBTRun()
+ {
+ //destructor
+
+ delete fEvents;//delete array with events
+ }
+
+
+/**************************************************************************/
+void AliHBTRun::Reset()
+ {
+ fEvents->Clear();//clear an array with events.
+ //All events are deleted because array is an owner (set in constructor)
+ }
+/**************************************************************************/
+
+
+void AliHBTRun::AddParticle(Int_t event, AliHBTParticle* part)
+{
+ //Adds particle to event
+
+ //if there is no event of this number, crate it and add to the collection
+ if(!GetEvent(event)) fEvents->AddAtAndExpand(new AliHBTEvent, event);
+
+ GetEvent(event)->AddParticle(part);
+}
+/**************************************************************************/
+
+
+void AliHBTRun::AddParticle(Int_t event, TParticle* part)
+{
+ //if there is no event of this number, crate it and add to the collection
+ if(!GetEvent(event)) fEvents->AddAtAndExpand(new AliHBTEvent, event);
+ GetEvent(event)->AddParticle(part);
+}
+/**************************************************************************/
+
+
+void AliHBTRun::AddParticle(Int_t event, Int_t pdg,
+ Double_t px, Double_t py, Double_t pz, Double_t etot,
+ Double_t vx, Double_t vy, Double_t vz, Double_t time)
+{
+ //if there is no event of this number, crate it and add to the collection
+ if(!GetEvent(event)) fEvents->AddAtAndExpand(new AliHBTEvent, event);
+ GetEvent(event)->AddParticle(pdg,px,py,pz,etot,vx,vy,vz,time);
+}
+
+/**************************************************************************/
--- /dev/null
+#ifndef ALIHBTRUN_H
+#define ALIHBTRUN_H
+
+#include "AliHBTEvent.h"
+#include <TObjArray.h>
+
+//class describing set of events (the run)
+//designed for fast acces
+//Piotr.Skowronski@cern.ch
+
+class AliHBTParticle;
+
+class AliHBTEvent;
+
+class AliHBTRun: public TObject
+ {
+ public:
+ AliHBTRun();
+ virtual ~AliHBTRun();
+
+ void AddParticle(Int_t event, AliHBTParticle* part); //inerface to AliHBTEvent::AddParticle(AliHBTParticle*)
+ void AddParticle(Int_t event, TParticle* part);//inerface to AliHBTEvent::AddParticle(TParticle*)
+
+ //inerface to AliHBTEvent::AddParticle(Int_t.Double_t,Double_t,Double_t,Double_t,Double_t,Double_t,Double_t,Double_t,Double_t)
+ void AddParticle(Int_t event, Int_t pdg,
+ Double_t px, Double_t py, Double_t pz, Double_t etot,
+ Double_t vx, Double_t vy, Double_t vz, Double_t time);
+
+ AliHBTParticle* GetParticle(Int_t event, Int_t n); //returns nth particle from event
+ AliHBTEvent* GetEvent(Int_t event) const; //returns AliHBTEvent number "event"
+
+ Int_t GetNumberOfEvents() const; //returns number of events
+ Int_t GetNumberOfParticlesInEvent(Int_t event) const; //returns number of particles in event number "event"
+ void Reset();//clears all events in the array (deletes)
+ protected:
+ TObjArray* fEvents;//!Array containig AliHBTEvents
+ private:
+
+ public:
+ ClassDef(AliHBTRun,1)
+ };
+
+
+/**************************************************************************/
+
+inline
+AliHBTEvent* AliHBTRun::GetEvent(Int_t event) const
+{
+//returns pointer to AliHBTEvent number "event"
+
+ return (AliHBTEvent*)fEvents->At(event);
+}
+/**************************************************************************/
+inline
+AliHBTParticle* AliHBTRun::GetParticle(Int_t event, Int_t n)
+{
+ //returns nth particle from event number event
+ return GetEvent(event)->GetParticle(n);
+}
+
+/**************************************************************************/
+
+inline
+Int_t AliHBTRun::GetNumberOfEvents() const
+ {
+//returns number of events in collection
+
+ return fEvents->GetEntriesFast(); //there may be empty slots but we do not care
+ //Analysis checks it if return is not NULL
+ }
+/**************************************************************************/
+
+inline
+Int_t AliHBTRun::GetNumberOfParticlesInEvent(Int_t event) const
+{
+//returns number of Particles in event
+ return GetEvent(event)->GetNumberOfParticles();
+}
+
+#endif
--- /dev/null
+#ifdef __CINT__
+
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+
+#pragma link C++ class AliHBTAnalysis+;
+#pragma link C++ class AliHBTParticle+;
+#pragma link C++ class AliHBTEvent+;
+#pragma link C++ class AliHBTRun+;
+#pragma link C++ class AliHBTFunction+;
+#pragma link C++ class AliHBTTwoPartFctn+;
+#pragma link C++ class AliHBTFourPartFctn+;
+#pragma link C++ class AliHBTTwoPartFctn1D+;
+#pragma link C++ class AliHBTTwoPartFctn2D+;
+#pragma link C++ class AliHBTTwoPartFctn3D+;
+
+#pragma link C++ class AliHBTFourPartFctn2D+;
+
+#pragma link C++ class AliHBTPair+;
+#pragma link C++ class AliHBTParticleCut-;
+#pragma link C++ class AliHBTEmptyParticleCut-;
+#pragma link C++ class AliHBTPairCut-;
+#pragma link C++ class AliHBTEmptyPairCut-;
+
+#pragma link C++ class AliHbtBaseCut+;
+#pragma link C++ class AliHbtBasePairCut+;
+#pragma link C++ class AliHBTQInvCut+;
+
+#pragma link C++ class AliHBTMomentumCut+;
+#pragma link C++ class AliHBTPtCut+;
+#pragma link C++ class AliHBTEnergyCut+;
+#pragma link C++ class AliHBTRapidityCut+;
+#pragma link C++ class AliHBTPseudoRapidityCut+;
+#pragma link C++ class AliHBTPxCut+;
+#pragma link C++ class AliHBTPyCut+;
+#pragma link C++ class AliHBTPzCut+;
+#pragma link C++ class AliHBTPhiCut+;
+#pragma link C++ class AliHBTThetaCut+;
+#pragma link C++ class AliHBTVxCut+;
+#pragma link C++ class AliHBTVyCut+;
+#pragma link C++ class AliHBTVzCut+;
+
+
+#pragma link C++ class AliHBTReader+;
+#pragma link C++ class AliHBTReaderTPC+;
+#pragma link C++ class AliHBTReaderPPprod+;
+
+#pragma link C++ class AliHBTQInvCorrelFctn+;
+#pragma link C++ class AliHBTInvMassCorrelFctn+;
+
+#pragma link C++ class AliHBTQOutResolVSQInvFctn+;
+#pragma link C++ class AliHBTQSideResolVSQInvFctn+;
+#pragma link C++ class AliHBTQLongResolVSQInvFctn+;
+#pragma link C++ class AliHBTQInvResolVSKtFctn+;
+#pragma link C++ class AliHBTQOutResolVSKtFctn+;
+#pragma link C++ class AliHBTQSideResolVSKtFctn+;
+#pragma link C++ class AliHBTQLongResolVSKtFctn+;
+#pragma link C++ class AliHBTQOutResolVSQOutFctn+;
+#pragma link C++ class AliHBTQSideResolVSQSideFctn+;
+#pragma link C++ class AliHBTQLongResolVSQLongFctn+;
+
+
+
+#endif
--- /dev/null
+############################### HBTAnalysis Makefile ##################################
+
+# Include machine specific definitions
+
+include $(ALICE_ROOT)/conf/GeneralDef
+include $(ALICE_ROOT)/conf/MachineDef.$(ALICE_TARGET)
+
+
+OPT = -g
+PACKAGE = HBTAnalysis
+
+# C++ sources
+
+SRCS = AliHBTAnalysis.cxx AliHBTFunction.cxx \
+ AliHBTEvent.cxx AliHBTRun.cxx \
+ AliHBTParticle.cxx AliHBTParticleCut.cxx \
+ AliHBTPair.cxx AliHBTPairCut.cxx\
+ AliHBTCorrelFctn.cxx \
+ AliHBTReader.cxx AliHBTReaderTPC.cxx\
+ AliHBTQResolutionFctns.cxx AliHBTReaderPPprod.cxx
+
+# Fortran sources
+
+FSRCS =
+
+# C++ Headers
+
+HDRS = $(SRCS:.cxx=.h) HBTAnalysisLinkDef.h
+
+# Library dictionary
+
+DICT = HBTAnalysisCint.cxx
+DICTH = $(DICT:.cxx=.h)
+DICTO = $(patsubst %.cxx,tgt_$(ALICE_TARGET)/%.o,$(DICT))
+
+# FORTRAN Objectrs
+
+FOBJS = $(patsubst %.F,tgt_$(ALICE_TARGET)/%.o,$(FSRCS))
+
+# C Objects
+
+COBJS = $(patsubst %.c,tgt_$(ALICE_TARGET)/%.o,$(CSRCS))
+
+# C++ Objects
+
+OBJS = $(patsubst %.cxx,tgt_$(ALICE_TARGET)/%.o,$(SRCS)) $(DICTO)
+
+# C++ compilation flags
+
+INCLUDES = -I$(ALICE_ROOT)/TPC/ -I$(ALICE_ROOT)/CONTAINERS/
+
+CXXFLAGS = $(CXXOPTS) -I$(ROOTSYS)/include -I. -I$(ALICE_ROOT)/include/ $(INCLUDES)
+
+# FORTRAN compilation flags
+
+FFLAGS = $(FOPT) -I$(ALICE_ROOT)/GEANT321
+
+##### TARGETS #####
+
+# Target
+
+SLIBRARY = $(LIBDIR)/libHBTAnalysis.$(SL)
+ALIBRARY = $(LIBDIR)/libHBTAnalysis.a
+
+default: $(SLIBRARY)
+
+
+$(LIBDIR)/libHBTAnalysis.$(SL): $(OBJS)
+
+$(DICT): $(HDRS)
+
+DEPINC += -I$(ALICE_ROOT)/GEANT321
+
+depend: $(SRCS) $(FSRCS)
+
+TOCLEAN = $(OBJS) $(FOBJS) *Cint.cxx *Cint.h
+
+CHECKS = $(patsubst %.cxx,check/%.viol,$(SRCS))
+
+############################### General Macros ################################
+
+include $(ALICE_ROOT)/conf/GeneralMacros
+
+############################ Dependencies #####################################
+
+-include tgt_$(ALICE_TARGET)/Make-depend
+
+
--- /dev/null
+SRCS = \
+AliHBTAnalysis.cxx AliHBTPair.cxx AliHBTQResolutionFctns.cxx AliHBTReaderPPprod.cxx\
+AliHBTCorrelFctn.cxx AliHBTPairCut.cxx AliHBTReader.cxx AliHBTReaderTPC.cxx\
+AliHBTEvent.cxx AliHBTParticle.cxx AliHBTRun.cxx\
+AliHBTFunction.cxx AliHBTParticleCut.cxx
+
+
+HDRS:= $(SRCS:.cxx=.h)
+
+DHDR= HBTAnalysisLinkDef.h
+
+EINCLUDE:= TPC CONTAINERS
+
+