From cc5faabc22455ac82a15d545f3b97a33647ff27f Mon Sep 17 00:00:00 2001 From: akisiel Date: Fri, 2 Nov 2007 09:37:15 +0000 Subject: [PATCH] Port of changes from v4-07-Release and additional rule conformance --- .../AliFemtoUser/AliFemtoESDTrackCut.cxx | 164 +++++++++- .../AliFemtoUser/AliFemtoESDTrackCut.h | 22 +- .../AliFemtoEventReaderESDKine.cxx | 305 +++++------------- .../AliFemtoUser/AliFemtoEventReaderESDKine.h | 32 +- .../AliFemtoModelBPLCMSCorrFctn.cxx | 167 ++++------ .../AliFemtoModelBPLCMSCorrFctn.h | 25 +- .../AliFemtoModelGausRinvFreezeOutGenerator.h | 2 +- .../AliFemtoUser/AliFemtoQPairCut.cxx | 25 +- .../AliFemtoUser/AliFemtoQPairCut.h | 31 +- .../AliFemtoShareQualityCorrFctn.cxx | 99 ++++-- .../AliFemtoShareQualityPairCut.cxx | 61 +++- .../AliFemtoShareQualityPairCut.h | 39 ++- 12 files changed, 533 insertions(+), 439 deletions(-) diff --git a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoESDTrackCut.cxx b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoESDTrackCut.cxx index 438976d9fa2..741f0083b62 100644 --- a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoESDTrackCut.cxx +++ b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoESDTrackCut.cxx @@ -11,6 +11,9 @@ *************************************************************************** * * $Log$ + * Revision 1.3 2007/05/22 09:01:42 akisiel + * Add the possibiloity to save cut settings in the ROOT file + * * Revision 1.2 2007/05/21 10:38:25 akisiel * More coding rule conformance * @@ -47,6 +50,38 @@ ClassImp(AliFemtoESDTrackCut) #endif + +// electron +// 0.13 - 1.8 +// 0 7.594129e-02 8.256141e-03 +// 1 -5.535827e-01 8.170825e-02 +// 2 1.728591e+00 3.104210e-01 +// 3 -2.827893e+00 5.827802e-01 +// 4 2.503553e+00 5.736207e-01 +// 5 -1.125965e+00 2.821170e-01 +// 6 2.009036e-01 5.438876e-02 + +// pion +// 0.13 - 2.0 +// 0 1.063457e+00 8.872043e-03 +// 1 -4.222208e-01 2.534402e-02 +// 2 1.042004e-01 1.503945e-02 + +// kaon +// 0.18 - 2.0 +// 0 -7.289406e-02 1.686074e-03 +// 1 4.415666e-01 1.143939e-02 +// 2 -2.996790e-01 1.840964e-02 +// 3 6.704652e-02 7.783990e-03 + +// proton +// 0.26 - 2.0 +// 0 -3.730200e-02 2.347311e-03 +// 1 1.163684e-01 1.319316e-02 +// 2 8.354116e-02 1.997948e-02 +// 3 -4.608098e-02 8.336400e-03 + + AliFemtoESDTrackCut::AliFemtoESDTrackCut() : fCharge(0), fLabel(0), @@ -54,7 +89,9 @@ AliFemtoESDTrackCut::AliFemtoESDTrackCut() : fminTPCclsF(0), fminITScls(0), fNTracksPassed(0), - fNTracksFailed(0) + fNTracksFailed(0), + fRemoveKinks(kFALSE), + fMostProbable(0) { // Default constructor fNTracksPassed = fNTracksFailed = 0; @@ -70,18 +107,18 @@ AliFemtoESDTrackCut::AliFemtoESDTrackCut() : fStatus=0; fminTPCclsF=0; fminITScls=0; - } //------------------------------ -//AliFemtoESDTrackCut::~AliFemtoESDTrackCut(){ -// /* noop */ -//} +AliFemtoESDTrackCut::~AliFemtoESDTrackCut(){ + /* noop */ +} //------------------------------ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track) { // test the particle and return // true if it meets all the criteria // false if it doesn't meet at least one of the criteria + float tMost[5]; //cout<<"AliFemtoESD cut"<Label()<<" "<Flags()<<" "<TPCnclsF()<<" "<ITSncls()<Flags()&fStatus)!=fStatus) { - // cout<Flags()<<" "<Flags()<<" "<PidProbElectron() << " " +// << track->PidProbMuon() << " " +// << track->PidProbPion() << " " +// << track->PidProbKaon() << " " +// << track->PidProbProton() << " " +// << track->PidProbElectron()+track->PidProbMuon()+track->PidProbPion()+track->PidProbKaon()+track->PidProbProton() << endl; + + if ((track->PidProbElectron()PidProbElectron()>fPidProbElectron[1])) { fNTracksFailed++; @@ -181,6 +227,23 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track) //cout<KinkIndex(0)) || (track->KinkIndex(1)) || (track->KinkIndex(2))) + return false; + } + + if (fMostProbable) { + tMost[0] = track->PidProbElectron()*PidFractionElectron(track->P().mag()); + tMost[1] = 0.0; + tMost[2] = track->PidProbPion()*PidFractionPion(track->P().mag()); + tMost[3] = track->PidProbKaon()*PidFractionKaon(track->P().mag()); + tMost[4] = track->PidProbProton()*PidFractionProton(track->P().mag()); + int imost=0; + float ipidmax = 0.0; + for (int ip=0; ip<5; ip++) + if (tMost[ip] > ipidmax) { ipidmax = tMost[ip]; imost = ip; }; + if (imost != fMostProbable) return false; + } // cout<<"Go Through the cut"<Mass()); @@ -257,6 +321,92 @@ TList *AliFemtoESDTrackCut::ListSettings() tListSetttings->AddLast(new TObjString(buf)); snprintf(buf, 200, "AliFemtoESDTrackCut.rapidity.maximum=%lf", fRapidity[1]); tListSetttings->AddLast(new TObjString(buf)); - + snprintf(buf, 200, "AliFemtoESDTrackCut.removekinks=%i", fRemoveKinks); + tListSetttings->AddLast(new TObjString(buf)); + if (fMostProbable) { + if (fMostProbable == 2) + snprintf(buf, 200, "AliFemtoESDTrackCut.mostprobable=%s", "Pion"); + if (fMostProbable == 3) + snprintf(buf, 200, "AliFemtoESDTrackCut.mostprobable=%s", "Kaon"); + if (fMostProbable == 4) + snprintf(buf, 200, "AliFemtoESDTrackCut.mostprobable=%s", "Proton"); + tListSetttings->AddLast(new TObjString(buf)); + } return tListSetttings; } +void AliFemtoESDTrackCut::SetRemoveKinks(const bool& flag) +{ + fRemoveKinks = flag; +} + + // electron +// 0.13 - 1.8 +// 0 7.594129e-02 8.256141e-03 +// 1 -5.535827e-01 8.170825e-02 +// 2 1.728591e+00 3.104210e-01 +// 3 -2.827893e+00 5.827802e-01 +// 4 2.503553e+00 5.736207e-01 +// 5 -1.125965e+00 2.821170e-01 +// 6 2.009036e-01 5.438876e-02 +float AliFemtoESDTrackCut::PidFractionElectron(float mom) const +{ + // Provide a parameterized fraction of electrons dependent on momentum + if (mom<0.13) return 0.0; + if (mom>1.8) return 0.0; + return (7.594129e-02 + -5.535827e-01*mom + +1.728591e+00*mom*mom + -2.827893e+00*mom*mom*mom + +2.503553e+00*mom*mom*mom*mom + -1.125965e+00*mom*mom*mom*mom*mom + +2.009036e-01*mom*mom*mom*mom*mom*mom); +} + +// pion +// 0.13 - 2.0 +// 0 1.063457e+00 8.872043e-03 +// 1 -4.222208e-01 2.534402e-02 +// 2 1.042004e-01 1.503945e-02 +float AliFemtoESDTrackCut::PidFractionPion(float mom) const +{ + // Provide a parameterized fraction of pions dependent on momentum + if (mom<0.13) return 0.0; + if (mom>2.0) return 0.0; + return ( 1.063457e+00 + -4.222208e-01*mom + +1.042004e-01*mom*mom); +} + +// kaon +// 0.18 - 2.0 +// 0 -7.289406e-02 1.686074e-03 +// 1 4.415666e-01 1.143939e-02 +// 2 -2.996790e-01 1.840964e-02 +// 3 6.704652e-02 7.783990e-03 +float AliFemtoESDTrackCut::PidFractionKaon(float mom) const +{ + // Provide a parameterized fraction of kaons dependent on momentum + if (mom<0.18) return 0.0; + if (mom>2.0) return 0.0; + return (-7.289406e-02 + +4.415666e-01*mom + -2.996790e-01*mom*mom + +6.704652e-02*mom*mom*mom); +} + +// proton +// 0.26 - 2.0 +// 0 -3.730200e-02 2.347311e-03 +// 1 1.163684e-01 1.319316e-02 +// 2 8.354116e-02 1.997948e-02 +// 3 -4.608098e-02 8.336400e-03 +float AliFemtoESDTrackCut::PidFractionProton(float mom) const +{ + // Provide a parameterized fraction of protons dependent on momentum + if (mom<0.26) return 0.0; + if (mom>2.0) return 0.0; + return (-3.730200e-02 + +1.163684e-01*mom + +8.354116e-02*mom*mom + -4.608098e-02*mom*mom*mom); +} diff --git a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoESDTrackCut.h b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoESDTrackCut.h index 4322b13c292..a90e5cf9214 100644 --- a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoESDTrackCut.h +++ b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoESDTrackCut.h @@ -22,7 +22,7 @@ class AliFemtoESDTrackCut : public AliFemtoTrackCut public: AliFemtoESDTrackCut(); - //~AliFemtoESDTrackCut(); + virtual ~AliFemtoESDTrackCut(); virtual bool Pass(const AliFemtoTrack* aTrack); @@ -31,7 +31,7 @@ class AliFemtoESDTrackCut : public AliFemtoTrackCut void SetPt(const float& lo, const float& hi); void SetRapidity(const float& lo, const float& hi); - void SetCharge(const int&); + void SetCharge(const int& ch); void SetPidProbElectron(const float& lo, const float& hi); void SetPidProbPion(const float& lo, const float& hi); void SetPidProbKaon(const float& lo, const float& hi); @@ -41,7 +41,12 @@ class AliFemtoESDTrackCut : public AliFemtoTrackCut void SetStatus(const long& w); void SetminTPCclsF(const short& s); void SetminITScls(const int& s); - + void SetRemoveKinks(const bool& flag); + void SetMostProbablePion(); + void SetMostProbableKaon(); + void SetMostProbableProton(); + void SetNoMostProbable(); + private: // here are the quantities I want to cut on... int fCharge; // particle charge @@ -58,6 +63,13 @@ class AliFemtoESDTrackCut : public AliFemtoTrackCut int fminITScls; // min number of clusters assigned in the ITS long fNTracksPassed; // passed tracks count long fNTracksFailed; // failed tracks count + bool fRemoveKinks; // if true particles with any kink label will not pass + int fMostProbable; // this particle type is required to be most probable + + float PidFractionElectron(float mom) const; + float PidFractionPion(float mom) const; + float PidFractionKaon(float mom) const; + float PidFractionProton(float mom) const; #ifdef __ROOT__ ClassDef(AliFemtoESDTrackCut, 1) @@ -77,6 +89,10 @@ inline void AliFemtoESDTrackCut::SetLabel(const bool& flag){fLabel=flag;} inline void AliFemtoESDTrackCut::SetStatus(const long& status){fStatus=status;} inline void AliFemtoESDTrackCut::SetminTPCclsF(const short& minTPCclsF){fminTPCclsF=minTPCclsF;} inline void AliFemtoESDTrackCut::SetminITScls(const int& minITScls){fminITScls=minITScls;} +inline void AliFemtoESDTrackCut::SetMostProbablePion() { fMostProbable = 2; } +inline void AliFemtoESDTrackCut::SetMostProbableKaon() { fMostProbable = 3; } +inline void AliFemtoESDTrackCut::SetMostProbableProton() { fMostProbable = 4; } +inline void AliFemtoESDTrackCut::SetNoMostProbable() { fMostProbable = 0; } #endif diff --git a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoEventReaderESDKine.cxx b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoEventReaderESDKine.cxx index 1714402225b..761795d5d45 100644 --- a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoEventReaderESDKine.cxx +++ b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoEventReaderESDKine.cxx @@ -11,6 +11,9 @@ /* *$Id$ *$Log$ + *Revision 1.1 2007/05/25 12:42:54 akisiel + *Adding a reader for the Kine information + * *Revision 1.2 2007/05/22 09:01:42 akisiel *Add the possibiloity to save cut settings in the ROOT file * @@ -34,9 +37,10 @@ #include "AliFemtoEventReaderESDKine.h" #include "TFile.h" -#include "TTree.h" -#include "AliESD.h" +#include "TChain.h" +#include "AliESDEvent.h" #include "AliESDtrack.h" +#include "AliESDVertex.h" #include "AliStack.h" #include "AliAODParticle.h" #include "TParticle.h" @@ -65,25 +69,11 @@ AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine(): fConstrained(true), fNumberofEvent(0), fCurEvent(0), - fCurFile(0), - fListOfFiles(0x0), + fCurRLEvent(0), fTree(0x0), fEvent(0x0), - fEsdFile(0x0), - fEventFriend(0), - fRunLoader(0x0), - fSharedList(0x0), - fClusterPerPadrow(0x0) + fRunLoader(0x0) { - // default constructor - fClusterPerPadrow = (list **) malloc(sizeof(list *) * AliESDfriendTrack::kMaxTPCcluster); - for (int tPad=0; tPad(); - } - fSharedList = (list **) malloc(sizeof(list *) * AliESDfriendTrack::kMaxTPCcluster); - for (int tPad=0; tPad(); - } } AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine(const AliFemtoEventReaderESDKine &aReader) : @@ -92,15 +82,10 @@ AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine(const AliFemtoEventReader fConstrained(true), fNumberofEvent(0), fCurEvent(0), - fCurFile(0), - fListOfFiles(0x0), + fCurRLEvent(0), fTree(0x0), fEvent(0x0), - fEsdFile(0x0), - fEventFriend(0), - fRunLoader(0x0), - fSharedList(0x0), - fClusterPerPadrow(0x0) + fRunLoader(0x0) { // copy constructor fInputFile = aReader.fInputFile; @@ -108,30 +93,7 @@ AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine(const AliFemtoEventReader fConstrained = aReader.fConstrained; fNumberofEvent = aReader.fNumberofEvent; fCurEvent = aReader.fCurEvent; - fCurFile = aReader.fCurFile; - fTree = aReader.fTree->CloneTree(); - // fEvent = new AliESD(*aReader.fEvent); - fEvent = new AliESD(); - fEsdFile = new TFile(aReader.fEsdFile->GetName()); - fEventFriend = aReader.fEventFriend; - fClusterPerPadrow = (list **) malloc(sizeof(list *) * AliESDfriendTrack::kMaxTPCcluster); - for (int tPad=0; tPad(); - list::iterator iter; - for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) { - fClusterPerPadrow[tPad]->push_back(*iter); - } - } - fSharedList = (list **) malloc(sizeof(list *) * AliESDfriendTrack::kMaxTPCcluster); - for (int tPad=0; tPad(); - list::iterator iter; - for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) { - fSharedList[tPad]->push_back(*iter); - } - } - for (unsigned int veciter = 0; veciterclear(); - delete fClusterPerPadrow[tPad]; - } - delete [] fClusterPerPadrow; - for (int tPad=0; tPadclear(); - delete fSharedList[tPad]; - } - delete [] fSharedList; } //__________________ @@ -168,57 +118,18 @@ AliFemtoEventReaderESDKine& AliFemtoEventReaderESDKine::operator=(const AliFemto fConstrained = aReader.fConstrained; fNumberofEvent = aReader.fNumberofEvent; fCurEvent = aReader.fCurEvent; - fCurFile = aReader.fCurFile; + fCurRLEvent = aReader.fCurRLEvent; if (fTree) delete fTree; - fTree = aReader.fTree->CloneTree(); + // fTree = aReader.fTree->CloneTree(); if (fEvent) delete fEvent; - fEvent = new AliESD(); - if (fEsdFile) delete fEsdFile; - fEsdFile = new TFile(aReader.fEsdFile->GetName()); + fEvent = new AliESDEvent(); if (fRunLoader) delete fRunLoader; fRunLoader = new AliRunLoader(*aReader.fRunLoader); - fEventFriend = aReader.fEventFriend; - - if (fClusterPerPadrow) { - for (int tPad=0; tPadclear(); - delete fClusterPerPadrow[tPad]; - } - delete [] fClusterPerPadrow; - } - - if (fSharedList) { - for (int tPad=0; tPadclear(); - delete fSharedList[tPad]; - } - delete [] fSharedList; - } - - fClusterPerPadrow = (list **) malloc(sizeof(list *) * AliESDfriendTrack::kMaxTPCcluster); - for (int tPad=0; tPad(); - list::iterator iter; - for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) { - fClusterPerPadrow[tPad]->push_back(*iter); - } - } - fSharedList = (list **) malloc(sizeof(list *) * AliESDfriendTrack::kMaxTPCcluster); - for (int tPad=0; tPad(); - list::iterator iter; - for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) { - fSharedList[tPad]->push_back(*iter); - } - } - for (unsigned int veciter = 0; veciterAddFile(buffer); delete tree; } esdFile->Close(); @@ -258,18 +172,6 @@ void AliFemtoEventReaderESDKine::SetInputFile(const char* inputFile) } } -//setting the next file to read -bool AliFemtoEventReaderESDKine::GetNextFile() -{ - // Begin reading the next file - if (fCurFile>=fListOfFiles.size()) - return false; - fFileName=fListOfFiles.at(fCurFile); - cout<<"FileName set on "<Get("esdTree"); - fTree->SetBranchAddress("ESD", &fEvent); - - // Attach the friend tree with additional information - tFriendFileName = fFileName.c_str(); - tFriendFileName.ReplaceAll("s.root","friends.root"); - // tFriendFileName.insert(tFriendFileName.find("s.root"),"friend"); - cout << "Reading friend " << tFriendFileName.Data() << endl;; - fTree->AddFriend("esdFriendTree",tFriendFileName.Data()); - fTree->SetBranchAddress("ESDfriend",&fEventFriend); +// fEsdFile=TFile::Open(fFileName.c_str(),"READ"); +// fTree = (TTree*) fEsdFile->Get("esdTree"); + + fTree->SetBranchStatus("MuonTracks*",0); + fTree->SetBranchStatus("PmdTracks*",0); + fTree->SetBranchStatus("TrdTracks*",0); + fTree->SetBranchStatus("V0s*",0); + fTree->SetBranchStatus("Cascades*",0); + fTree->SetBranchStatus("Kinks*",0); + fTree->SetBranchStatus("CaloClusters*",0); + fTree->SetBranchStatus("AliRawDataErrorLogs*",0); + fTree->SetBranchStatus("ESDfriend*",0); + fEvent->ReadFromTree(fTree); // chain->SetBranchStatus("*",0); // chain->SetBranchStatus("fUniqueID",1); // chain->SetBranchStatus("fTracks",1); // chain->SetBranchStatus("fTracks.*",1); // chain->SetBranchStatus("fTracks.fTPCindex[160]",1); - fTree->SetBranchStatus("fTracks.fCalibContainer",0); - - - fNumberofEvent=fTree->GetEntries(); - cout<<"Number of Entries in file "<LoadHeader()) - { - cout << "Could not read RunLoader header in file " << tGAliceFilename.Data() << endl; - exit(0); - } - fRunLoader->LoadKinematics(); - +// fTree->SetBranchStatus("fTracks.fCalibContainer",0); + + + fNumberofEvent=fTree->GetEntries(); + + if (fNumberofEvent == 0) { + cout<<"no event in input "<GetEvent(fCurEvent);//getting next event - fEvent->SetESDfriend(fEventFriend); // vector tLabelTable;//to check labels - fRunLoader->GetEvent(fCurEvent); + cout << "fFileName is " << fFileName.Data() << endl; + cout << "Current file is " << fTree->GetCurrentFile()->GetName() << endl; + if (fFileName.CompareTo(fTree->GetCurrentFile()->GetName())) { + fFileName = fTree->GetCurrentFile()->GetName(); + tGAliceFilename = fFileName; + tGAliceFilename.ReplaceAll("AliESDs","galice"); + cout << "Reading RunLoader from " << tGAliceFilename.Data() << endl; + if (fRunLoader) delete fRunLoader; + fRunLoader = AliRunLoader::Open(tGAliceFilename.Data()); + if (fRunLoader==0x0) + { + cout << "No Kine tree in file " << tGAliceFilename.Data() << endl; + exit(0); + } + if(fRunLoader->LoadHeader()) + { + cout << "Could not read RunLoader header in file " << tGAliceFilename.Data() << endl; + exit(0); + } + fRunLoader->LoadKinematics(); + fCurRLEvent = 0; + } + + fRunLoader->GetEvent(fCurRLEvent); AliStack* tStack = 0x0; tStack = fRunLoader->Stack(); @@ -387,36 +296,6 @@ AliFemtoEvent* AliFemtoEventReaderESDKine::ReturnHbtEvent() nofTracks=fEvent->GetNumberOfTracks(); int realnofTracks=0;//number of track which we use ina analysis - // Clear the shared cluster list - for (int tPad=0; tPadclear(); - } - for (int tPad=0; tPadclear(); - } - - - for (int i=0;iGetTrack(i);//getting next track - - list::iterator tClustIter; - - Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster]; - Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices); - for (int tNcl=0; tNcl= 0) { - tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]); - if (tClustIter == fClusterPerPadrow[tNcl]->end()) { - fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]); - } - else { - fSharedList[tNcl]->push_back(tTrackIndices[tNcl]); - } - } - } - - } - for (int i=0;iGetMagneticField())*kilogauss,(double)(trackCopy->Charge())); + AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); trackCopy->SetHelix(helix); trackCopy->SetTrackId(esdtrack->GetID()); trackCopy->SetFlags(esdtrack->GetStatus()); - //trackCopy->SetLabel(esdtrack->GetLabel()); + trackCopy->SetLabel(esdtrack->GetLabel()); //some stuff which could be useful float impact[2]; @@ -477,33 +356,9 @@ AliFemtoEvent* AliFemtoEventReaderESDKine::ReturnHbtEvent() trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma()); - // Fill cluster per padrow information - Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster]; - Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices); - for (int tNcl=0; tNcl 0) - trackCopy->SetTPCcluster(tNcl, 1); - else - trackCopy->SetTPCcluster(tNcl, 0); - } - - // Fill shared cluster information - list::iterator tClustIter; - - for (int tNcl=0; tNcl 0) { - tClustIter = find(fSharedList[tNcl]->begin(), fSharedList[tNcl]->end(), tTrackIndices[tNcl]); - if (tClustIter != fSharedList[tNcl]->end()) { - trackCopy->SetTPCshared(tNcl, 1); - cout << "Event next" << endl; - cout << "Track: " << i << endl; - cout << "Shared cluster: " << tNcl << " " << tTrackIndices[tNcl] << endl; - } - else { - trackCopy->SetTPCshared(tNcl, 0); - } - } - } + trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap()); + trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap()); + // Fill the hidden information with the simulated data TParticle *tPart = tStack->Particle(TMath::Abs(esdtrack->GetLabel())); AliAODParticle* tParticle= new AliAODParticle(*tPart,i); @@ -537,12 +392,12 @@ AliFemtoEvent* AliFemtoEventReaderESDKine::ReturnHbtEvent() hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event fCurEvent++; + fCurRLEvent++; cout<<"end of reading nt "<Reset(); delete fTree; - fEsdFile->Close(); } return hbtEvent; } diff --git a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoEventReaderESDKine.h b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoEventReaderESDKine.h index 013470e25b4..5a8acf70985 100644 --- a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoEventReaderESDKine.h +++ b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoEventReaderESDKine.h @@ -12,6 +12,9 @@ /* *$Id$ *$Log$ + *Revision 1.1 2007/05/25 12:42:54 akisiel + *Adding a reader for the Kine information + * *Revision 1.1 2007/05/16 10:22:11 akisiel *Making the directory structure of AliFemto flat. All files go into one common directory * @@ -27,16 +30,15 @@ */ -#ifndef ALIFEMTOEVENTREADERESD_H -#define ALIFEMTOEVENTREADERESD_H +#ifndef ALIFEMTOEVENTREADERESDKINE_H +#define ALIFEMTOEVENTREADERESDKINE_H #include "AliFemtoEventReader.h" #include "AliFemtoEnumeration.h" #include #include -#include "TTree.h" -#include "AliESD.h" -#include "AliESDfriend.h" +#include "TChain.h" +#include "AliESDEvent.h" #include #include "AliRunLoader.h" #include "AliFemtoModelHiddenInfo.h" @@ -53,7 +55,7 @@ class AliFemtoEventReaderESDKine : public AliFemtoEventReader AliFemtoEventReaderESDKine& operator=(const AliFemtoEventReaderESDKine& aReader); AliFemtoEvent* ReturnHbtEvent(); - AliFemtoString Report(); + AliFemtoString Report() const; //void SetFileName(const char* fileName); void SetInputFile(const char* inputFile); void SetConstrained(const bool constrained); @@ -62,23 +64,15 @@ class AliFemtoEventReaderESDKine : public AliFemtoEventReader protected: private: - bool GetNextFile(); // setting next file to read - - string fInputFile; // name of input file with ESD filenames - string fFileName; // name of current ESD file + TString fInputFile; // name of input file with ESD filenames + TString fFileName; // name of current ESD file bool fConstrained; // flag to set which momentum from ESD file will be use int fNumberofEvent; // number of Events in ESD file int fCurEvent; // number of current event - unsigned int fCurFile; // number of current file - vector fListOfFiles; // list of ESD files - TTree* fTree; // ESD tree - AliESD* fEvent; // ESD event - TFile* fEsdFile; // ESD file - AliESDfriend* fEventFriend; // ESD friend informaion + int fCurRLEvent; // Current simulated event + TChain* fTree; // ESD tree + AliESDEvent* fEvent; // ESD event AliRunLoader* fRunLoader; // Run loader for kine reading - - list **fSharedList; //! Table (one list per padrow) of clusters which are shared - list **fClusterPerPadrow; //! Table (one list per padrow) of clusters in each padrow #ifdef __ROOT__ ClassDef(AliFemtoEventReaderESDKine, 1) diff --git a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelBPLCMSCorrFctn.cxx b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelBPLCMSCorrFctn.cxx index 0c68d6d3077..628f9bfe5a5 100644 --- a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelBPLCMSCorrFctn.cxx +++ b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelBPLCMSCorrFctn.cxx @@ -7,6 +7,8 @@ /// /// //////////////////////////////////////////////////////////////////////////////// #include "AliFemtoModelBPLCMSCorrFctn.h" +#include "AliFemtoPair.h" +#include "AliFemtoModelManager.h" #include #ifdef __ROOT__ @@ -16,42 +18,57 @@ ClassImp(AliFemtoModelBPLCMSCorrFctn) //____________________________ AliFemtoModelBPLCMSCorrFctn::AliFemtoModelBPLCMSCorrFctn(char* title, const int& nbins, const float& QLo, const float& QHi) : - fNumerator(0), - fDenominator(0), + AliFemtoModelCorrFctn(title, nbins, QLo, QHi), + fNumerator3DTrue(0), + fNumerator3DFake(0), + fDenominator3D(0), fQinvHisto(0) { - // set up numerator - char TitNum[100] = "Num"; - strcat(TitNum,title); - fNumerator = new TH3D(TitNum,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi); + // set up true numerator + char TitNumT[100] = "Num3DTrue"; + strcat(TitNumT,title); + fNumerator3DTrue = new TH3D(TitNumT,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi); + // set up fake numerator + char TitNumF[100] = "Num3DFake"; + strcat(TitNumF,title); + fNumerator3DFake = new TH3D(TitNumF,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi); // set up denominator - char TitDen[100] = "Den"; + char TitDen[100] = "Den3D"; strcat(TitDen,title); - fDenominator = new TH3D(TitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi); + fDenominator3D = new TH3D(TitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi); // set up ave qInv char TitQinv[100] = "Qinv"; strcat(TitQinv,title); fQinvHisto = new TH3D(TitQinv,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi); // to enable error bar calculation... - fNumerator->Sumw2(); - fDenominator->Sumw2(); + fNumerator3DTrue->Sumw2(); + fNumerator3DFake->Sumw2(); + fDenominator3D->Sumw2(); } AliFemtoModelBPLCMSCorrFctn::AliFemtoModelBPLCMSCorrFctn(const AliFemtoModelBPLCMSCorrFctn& aCorrFctn) : - fNumerator(0), - fDenominator(0), + AliFemtoModelCorrFctn(aCorrFctn), + fNumerator3DTrue(0), + fNumerator3DFake(0), + fDenominator3D(0), fQinvHisto(0) { - fNumerator = new TH3D(*aCorrFctn.fNumerator); - fDenominator = new TH3D (*aCorrFctn.fDenominator); - fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto); + fNumerator3DTrue = new TH3D(*aCorrFctn.fNumerator3DTrue); + fNumerator3DFake = new TH3D(*aCorrFctn.fNumerator3DFake); + fDenominator3D = new TH3D(*aCorrFctn.fDenominator3D); + fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto); } //____________________________ -AliFemtoModelBPLCMSCorrFctn::~AliFemtoModelBPLCMSCorrFctn(){ - delete fNumerator; - delete fDenominator; +AliFemtoModelBPLCMSCorrFctn::~AliFemtoModelBPLCMSCorrFctn() +{ + if (fNumeratorTrue) delete fNumeratorTrue; + if (fNumeratorFake) delete fNumeratorFake; + if (fDenominator) delete fDenominator; + delete fNumerator3DTrue; + delete fNumerator3DFake; + delete fDenominator3D; delete fQinvHisto; } //_________________________ @@ -59,10 +76,12 @@ AliFemtoModelBPLCMSCorrFctn& AliFemtoModelBPLCMSCorrFctn::operator=(const AliFem { if (this == &aCorrFctn) return *this; - if (fNumerator) delete fNumerator; - fNumerator = new TH3D(*aCorrFctn.fNumerator); - if (fDenominator) delete fDenominator; - fDenominator = new TH3D(*aCorrFctn.fDenominator); + if (fNumerator3DTrue) delete fNumerator3DTrue; + fNumerator3DTrue = new TH3D(*aCorrFctn.fNumerator3DTrue); + if (fNumerator3DFake) delete fNumerator3DFake; + fNumerator3DFake = new TH3D(*aCorrFctn.fNumerator3DFake); + if (fDenominator3D) delete fDenominator3D; + fDenominator3D = new TH3D(*aCorrFctn.fDenominator3D); if (fQinvHisto) delete fQinvHisto; fQinvHisto = new TH3D(*aCorrFctn.fQinvHisto); @@ -70,10 +89,11 @@ AliFemtoModelBPLCMSCorrFctn& AliFemtoModelBPLCMSCorrFctn::operator=(const AliFem } //_________________________ -void AliFemtoModelBPLCMSCorrFctn::WriteOutHistos(){ - - fNumerator->Write(); - fDenominator->Write(); +void AliFemtoModelBPLCMSCorrFctn::Write(){ + AliFemtoModelCorrFctn::Write(); + fNumerator3DTrue->Write(); + fNumerator3DFake->Write(); + fDenominator3D->Write(); fQinvHisto->Write(); } @@ -84,20 +104,12 @@ void AliFemtoModelBPLCMSCorrFctn::Finish(){ //____________________________ AliFemtoString AliFemtoModelBPLCMSCorrFctn::Report(){ - string stemp = "LCMS Frame Bertsch-Pratt 3D Correlation Function Report:\n"; + string stemp = "LCMS Frame Bertsch-Pratt 3D Model Correlation Function Report:\n"; char ctemp[100]; - sprintf(ctemp,"Number of entries in numerator:\t%E\n",fNumerator->GetEntries()); + sprintf(ctemp,"Number of entries in numerator:\t%E\n",fNumeratorTrue->GetEntries()); stemp += ctemp; sprintf(ctemp,"Number of entries in denominator:\t%E\n",fDenominator->GetEntries()); stemp += ctemp; - sprintf(ctemp,"Number of entries in ratio:\t%E\n",fRatio->GetEntries()); - stemp += ctemp; - sprintf(ctemp,"Normalization region in Qinv was:\t%E\t%E\n",fQinvNormLo,fQinvNormHi); - stemp += ctemp; - sprintf(ctemp,"Number of pairs in Normalization region was:\n"); - stemp += ctemp; - sprintf(ctemp,"In numerator:\t%lu\t In denominator:\t%lu\n",fNumRealsNorm,fNumMixedNorm); - stemp += ctemp; /* if (fCorrection) { float radius = fCorrection->GetRadius(); @@ -110,79 +122,40 @@ AliFemtoString AliFemtoModelBPLCMSCorrFctn::Report(){ stemp += ctemp; */ - if (fPairCut){ - sprintf(ctemp,"Here is the PairCut specific to this CorrFctn\n"); - stemp += ctemp; - stemp += fPairCut->Report(); - } - else{ - sprintf(ctemp,"No PairCut specific to this CorrFctn\n"); - stemp += ctemp; - } - // AliFemtoString returnThis = stemp; return returnThis; } //____________________________ void AliFemtoModelBPLCMSCorrFctn::AddRealPair( AliFemtoPair* pair){ + Double_t weight = fManager->GetWeight(pair); - if (fPairCut){ - if (!(fPairCut->Pass(pair))) return; - } - - double Qinv = fabs(pair->qInv()); // note - qInv() will be negative for identical pairs... - if ((Qinv < fQinvNormHi) && (Qinv > fQinvNormLo)) fNumRealsNorm++; - double qOut = fabs(pair->qOutCMS()); - double qSide = fabs(pair->qSideCMS()); - double qLong = fabs(pair->qLongCMS()); + double qOut = fabs(pair->QOutCMS()); + double qSide = fabs(pair->QSideCMS()); + double qLong = fabs(pair->QLongCMS()); - fNumerator->Fill(qOut,qSide,qLong); + fNumerator3DTrue->Fill(qOut, qSide, qLong, weight); + fNumeratorTrue->Fill(pair->QInv(), weight); } //____________________________ void AliFemtoModelBPLCMSCorrFctn::AddMixedPair( AliFemtoPair* pair){ + Double_t weight = fManager->GetWeight(pair); - if (fPairCut){ - if (!(fPairCut->Pass(pair))) return; - } - - // double CoulombWeight = (fCorrection ? fCorrection->CoulombCorrect(pair) : 1.0); - double CoulombWeight = 1.0; + double qOut = fabs(pair->QOutCMS()); + double qSide = fabs(pair->QSideCMS()); + double qLong = fabs(pair->QLongCMS()); - double Qinv = fabs(pair->qInv()); // note - qInv() will be negative for identical pairs... - if ((Qinv < fQinvNormHi) && (Qinv > fQinvNormLo)) fNumMixedNorm++; - double qOut = fabs(pair->qOutCMS()); - double qSide = fabs(pair->qSideCMS()); - double qLong = fabs(pair->qLongCMS()); + fNumerator3DFake->Fill(qOut, qSide, qLong, weight); + fDenominator3D->Fill(qOut, qSide, qLong, 1.0); + fNumeratorFake->Fill(pair->QInv(), weight); + fDenominator->Fill(pair->QInv(), 1.0); - fDenominator->Fill(qOut,qSide,qLong,CoulombWeight); - // fUncorrectedDenominator->Fill(qOut,qSide,qLong,1.0); - fQinvHisto->Fill(qOut,qSide,qLong,Qinv); - - /* - // now for the momentum resolution stuff... - if (fSmearPair){ - double CorrWeight = 1.0 + - fLambda*exp((-qOut*qOut*fRout2 -qSide*qSide*fRside2 -qLong*qLong*fRlong2)/0.038936366329); - CorrWeight *= CoulombWeight; // impt. - - fIDNumHisto->Fill(qOut,qSide,qLong,CorrWeight); - fIDDenHisto->Fill(qOut,qSide,qLong,CoulombWeight); - - fSmearPair->SetUnsmearedPair(pair); - double qOut_prime = fabs(fSmearPair->SmearedPair().qOutCMS()); - double qSide_prime = fabs(fSmearPair->SmearedPair().qSideCMS()); - double qLong_prime = fabs(fSmearPair->SmearedPair().qLongCMS()); - - fSMNumHisto->Fill(qOut_prime,qSide_prime,qLong_prime,CorrWeight); - - double SmearedCoulombWeight = ( fCorrection ? - fCorrection->CoulombCorrect(&(fSmearPair->SmearedPair())) : - 1.0); - - fSMDenHisto->Fill(qOut_prime,qSide_prime,qLong_prime,SmearedCoulombWeight); - } - */ } - - +//_______________________ +AliFemtoModelCorrFctn* AliFemtoModelBPLCMSCorrFctn::Clone() +{ + // Clone the correlation function + AliFemtoModelBPLCMSCorrFctn *tCopy = new AliFemtoModelBPLCMSCorrFctn(*this); + + return tCopy; +} diff --git a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelBPLCMSCorrFctn.h b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelBPLCMSCorrFctn.h index 56fdb4d13d2..47a5bd6ec14 100644 --- a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelBPLCMSCorrFctn.h +++ b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelBPLCMSCorrFctn.h @@ -6,15 +6,15 @@ /// Authors: Adam Kisiel, kisiel@mps.ohio-state.edu /// /// /// //////////////////////////////////////////////////////////////////////////////// -#ifndef AliFemtoModelBPLCMS3DCorrFctn_hh -#define AliFemtoModelBPLCMS3DCorrFctn_hh +#ifndef AliFemtoModelBPLCMSCorrFctn_hh +#define AliFemtoModelBPLCMSCorrFctn_hh #include "AliFemtoCorrFctn.h" #include "AliFemtoModelCorrFctn.h" #include "AliFemtoPairCut.h" #include "TH3D.h" -class AliFemtoModelBPLCMS3DCorrFctn : public AliFemtoModelCorrFctn { +class AliFemtoModelBPLCMSCorrFctn : public AliFemtoModelCorrFctn { public: AliFemtoModelBPLCMSCorrFctn(); AliFemtoModelBPLCMSCorrFctn(char* title, const int& nbins, const float& QLo, const float& QHi); @@ -31,26 +31,19 @@ class AliFemtoModelBPLCMS3DCorrFctn : public AliFemtoModelCorrFctn { virtual void Write(); - virtual AliFemtoModelCorrFctnSource* Clone(); + virtual AliFemtoModelCorrFctn* Clone(); - TH3D* Numerator(); - TH3D* Denominator(); - TH3D* QinvHisto(); +protected: + TH3D* fNumerator3DTrue; // 3D Numerator with pairs from same event only + TH3D* fNumerator3DFake; // 3D Numerator with pairs from mixed events + TH3D* fDenominator3D; // 3D Denominator with the weight of 1.0 -private: - TH3D* fNumerator; - TH3D* fDenominator; - - TH3D* fQinvHisto; + TH3D* fQinvHisto; // Averag qinv histogram #ifdef __ROOT__ ClassDef(AliFemtoModelBPLCMSCorrFctn, 1) #endif }; -inline TH3D* AliFemtoModelBPLCMSCorrFctn::Numerator(){return fNumerator;} -inline TH3D* AliFemtoModelBPLCMSCorrFctn::Denominator(){return fDenominator;} -inline TH3D* AliFemtoModelBPLCMSCorrFctn::QinvHisto(){return fQinvHisto;} - #endif diff --git a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelGausRinvFreezeOutGenerator.h b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelGausRinvFreezeOutGenerator.h index 224f68d6109..281a380ff92 100644 --- a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelGausRinvFreezeOutGenerator.h +++ b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelGausRinvFreezeOutGenerator.h @@ -18,7 +18,7 @@ class AliFemtoModelGausRinvFreezeOutGenerator : public AliFemtoModelFreezeOutGen AliFemtoModelGausRinvFreezeOutGenerator(); AliFemtoModelGausRinvFreezeOutGenerator(const AliFemtoModelGausRinvFreezeOutGenerator &aModel); virtual ~AliFemtoModelGausRinvFreezeOutGenerator(); - virtual void GenerateFreezeOut(AliFemtoPair *aPair);; + virtual void GenerateFreezeOut(AliFemtoPair *aPair); void SetSizeInv(Double_t aSizeInv); diff --git a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoQPairCut.cxx b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoQPairCut.cxx index 22919bb1934..8854c7fc83a 100644 --- a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoQPairCut.cxx +++ b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoQPairCut.cxx @@ -1,3 +1,9 @@ +///////////////////////////////////////////////////////////////////////////// +// // +// AliFemtoQPairCut - a simple cut which selects pairs based on the values // +// of their respective q components / +// // +///////////////////////////////////////////////////////////////////////////// /*************************************************************************** * * $Id$ @@ -7,6 +13,12 @@ *************************************************************************** * * $Log$ + * Revision 1.2.6.1 2007/11/01 17:10:38 akisiel + * Fix code rule conformace + * + * Revision 1.2 2007/05/22 09:01:42 akisiel + * Add the possibiloity to save cut settings in the ROOT file + * * Revision 1.1 2007/05/16 10:25:06 akisiel * Making the directory structure of AliFemtoUser flat. All files go into one common directory * @@ -38,6 +50,7 @@ AliFemtoQPairCut::AliFemtoQPairCut(): fNPairsPassed(0), fNPairsFailed(0) { + // Default constructor fNPairsPassed = fNPairsFailed = 0; fQlong[0]=-1.0; fQlong[1]=100.0; fQout[0]=-1.0; fQout[1]=100.0; @@ -52,6 +65,7 @@ AliFemtoQPairCut::~AliFemtoQPairCut() //__________________ bool AliFemtoQPairCut::Pass(const AliFemtoPair* pair) { + // Select pairs based on their q values //bool temp = true; //temp ? fNPairsPassed++ : fNPairsFailed++; if ((fabs(pair->QLongCMS())QLongCMS())>fQlong[1])) @@ -80,11 +94,12 @@ bool AliFemtoQPairCut::Pass(const AliFemtoPair* pair) //__________________ AliFemtoString AliFemtoQPairCut::Report() { - string Stemp = "AliFemtoQ Pair Cut \n"; - char Ctemp[100]; - sprintf(Ctemp,"Number of pairs which passed:\t%ld Number which failed:\t%ld\n",fNPairsPassed,fNPairsFailed); - Stemp += Ctemp; - AliFemtoString returnThis = Stemp; + // Prepare a report + string stemp = "AliFemtoQ Pair Cut \n"; + char ctemp[100]; + sprintf(ctemp,"Number of pairs which passed:\t%ld Number which failed:\t%ld\n",fNPairsPassed,fNPairsFailed); + stemp += ctemp; + AliFemtoString returnThis = stemp; return returnThis; } //__________________ diff --git a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoQPairCut.h b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoQPairCut.h index 90fe35a46aa..5400876c966 100644 --- a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoQPairCut.h +++ b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoQPairCut.h @@ -1,3 +1,10 @@ +///////////////////////////////////////////////////////////////////////////// +// // +// AliFemtoQPairCut - a simple cut which selects pairs based on the values // +// of their respective q components / +// // +///////////////////////////////////////////////////////////////////////////// + /*************************************************************************** * * $Id $ @@ -8,6 +15,12 @@ *************************************************************************** * * $Log$ + * Revision 1.2.6.1 2007/11/01 17:10:38 akisiel + * Fix code rule conformace + * + * Revision 1.2 2007/05/22 09:01:42 akisiel + * Add the possibiloity to save cut settings in the ROOT file + * * Revision 1.1 2007/05/16 10:25:06 akisiel * Making the directory structure of AliFemtoUser flat. All files go into one common directory * @@ -21,8 +34,8 @@ **************************************************************************/ -#ifndef AliFemtoQPairCut_hh -#define AliFemtoQPairCut_hh +#ifndef ALIFEMTOQPAIRCUT_H +#define ALIFEMTOQPAIRCUT_H // do I need these lines ? //#ifndef StMaker_H @@ -36,7 +49,7 @@ public: AliFemtoQPairCut(); ~AliFemtoQPairCut(); - virtual bool Pass(const AliFemtoPair*); + virtual bool Pass(const AliFemtoPair* pair); virtual AliFemtoString Report(); virtual TList *ListSettings(); @@ -48,12 +61,12 @@ public: private: - long fNPairsPassed; - long fNPairsFailed; - float fQlong[2]; - float fQout[2]; - float fQside[2]; - float fQinv[2]; + long fNPairsPassed; // Number of pairs that passed the cut + long fNPairsFailed; // Number of pairs that failed the cut + float fQlong[2]; // Qlong range + float fQout[2]; // Qout range + float fQside[2]; // Qside range + float fQinv[2]; // Qinv range #ifdef __ROOT__ diff --git a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityCorrFctn.cxx b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityCorrFctn.cxx index c9dc4f73077..0effcc86d4e 100644 --- a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityCorrFctn.cxx +++ b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityCorrFctn.cxx @@ -25,21 +25,21 @@ AliFemtoShareQualityCorrFctn::AliFemtoShareQualityCorrFctn(char* title, const in // title = "Num Qinv (MeV/c)"; char tTitNum[100] = "NumShare"; strcat(tTitNum,title); - fShareNumerator = new TH2D(tTitNum,title,nbins,QinvLo,QinvHi,50,0.0,1.00001); + fShareNumerator = new TH2D(tTitNum,title,nbins,QinvLo,QinvHi,100,0.0,1.00001); // set up denominator //title = "Den Qinv (MeV/c)"; char tTitDen[100] = "DenShare"; strcat(tTitDen,title); - fShareDenominator = new TH2D(tTitDen,title,nbins,QinvLo,QinvHi,50,0.0,1.00001); + fShareDenominator = new TH2D(tTitDen,title,nbins,QinvLo,QinvHi,100,0.0,1.00001); char tTit2Num[100] = "NumQuality"; strcat(tTit2Num,title); - fQualityNumerator = new TH2D(tTit2Num,title,nbins,QinvLo,QinvHi,75,-0.500001,1.000001); + fQualityNumerator = new TH2D(tTit2Num,title,nbins,QinvLo,QinvHi,150,-0.500001,1.000001); // set up denominator //title = "Den Qinv (MeV/c)"; char tTit2Den[100] = "DenQuality"; strcat(tTit2Den,title); - fQualityDenominator = new TH2D(tTit2Den,title,nbins,QinvLo,QinvHi,75,-0.500001,1.000001); + fQualityDenominator = new TH2D(tTit2Den,title,nbins,QinvLo,QinvHi,150,-0.500001,1.000001); // set up ratio //title = "Ratio Qinv (MeV/c)"; // this next bit is unfortunately needed so that we can have many histos of same "title" @@ -135,7 +135,7 @@ AliFemtoString AliFemtoShareQualityCorrFctn::Report(){ //____________________________ void AliFemtoShareQualityCorrFctn::AddRealPair( AliFemtoPair* pair){ // add real (effect) pair - double tQinv = fabs(pair->QInv()); // note - qInv() will be negative for identical pairs... + // double tQinv = fabs(pair->QInv()); // note - qInv() will be negative for identical pairs... Int_t nh = 0; Int_t an = 0; Int_t ns = 0; @@ -148,9 +148,9 @@ void AliFemtoShareQualityCorrFctn::AddRealPair( AliFemtoPair* pair){ if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap) && pair->Track2()->Track()->TPCsharing().TestBitNumber(imap)) { - if (tQinv < 0.01) { - cout << "Shared cluster in row " << imap << endl; - } +// if (tQinv < 0.01) { +// cout << "Shared cluster in row " << imap << endl; +// } an++; nh+=2; ns+=2; @@ -169,23 +169,23 @@ void AliFemtoShareQualityCorrFctn::AddRealPair( AliFemtoPair* pair){ nh++; } } - if (tQinv < 0.01) { - cout << "Qinv of the pair is " << tQinv << endl; - cout << "Clusters: " << endl; - for (unsigned int imap=0; imapTrack1()->Track()->TPCclusters().GetNbits(); imap++) { - cout << imap ; - if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 "; - else cout << " 0 " ; - if (pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 "; - else cout << " 0 " ; - cout << " "; - if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S "; - else cout << " X "; - if (pair->Track2()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S "; - else cout << " X "; - cout << endl; - } - } +// if (tQinv < 0.01) { +// cout << "Qinv of the pair is " << tQinv << endl; +// cout << "Clusters: " << endl; +// for (unsigned int imap=0; imapTrack1()->Track()->TPCclusters().GetNbits(); imap++) { +// cout << imap ; +// if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 "; +// else cout << " 0 " ; +// if (pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 "; +// else cout << " 0 " ; +// cout << " "; +// if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S "; +// else cout << " X "; +// if (pair->Track2()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S "; +// else cout << " X "; +// cout << endl; +// } +// } Float_t hsmval = 0.0; Float_t hsfval = 0.0; @@ -195,12 +195,49 @@ void AliFemtoShareQualityCorrFctn::AddRealPair( AliFemtoPair* pair){ hsfval = ns*1.0/nh; } - if (tQinv < 0.01) { - cout << "Quality Sharity " << hsmval << " " << hsfval << " " << pair->Track1()->Track() << " " << pair->Track2()->Track() << endl; - } +// if ((tQinv < 0.01) && (hsmval<-0.46)) { +// cout << "Quality Sharity " << hsmval << " " << hsfval << " " << pair->Track1()->Track() << " " << pair->Track2()->Track() << endl; +// cout << "Qinv of the pair is " << tQinv << endl; +// cout << "Clusters: " << endl; +// for (unsigned int imap=0; imapTrack1()->Track()->TPCclusters().GetNbits(); imap++) { +// cout << imap ; +// if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 "; +// else cout << " 0 " ; +// if (pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) cout << " 1 "; +// else cout << " 0 " ; +// cout << " "; +// if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S "; +// else cout << " X "; +// if (pair->Track2()->Track()->TPCsharing().TestBitNumber(imap)) cout << " S "; +// else cout << " X "; +// cout << endl; +// } +// cout << "Momentum1 " +// << pair->Track1()->Track()->P().x() << " " +// << pair->Track1()->Track()->P().y() << " " +// << pair->Track1()->Track()->P().z() << " " +// << pair->Track1()->Track()->Label() << " " +// << pair->Track1()->Track()->TrackId() << " " +// << pair->Track1()->Track()->Flags() << " " +// << pair->Track1()->Track()->KinkIndex(0) << " " +// << pair->Track1()->Track()->KinkIndex(1) << " " +// << pair->Track1()->Track()->KinkIndex(2) << " " +// << endl; +// cout << "Momentum2 " +// << pair->Track2()->Track()->P().x() << " " +// << pair->Track2()->Track()->P().y() << " " +// << pair->Track2()->Track()->P().z() << " " +// << pair->Track2()->Track()->Label() << " " +// << pair->Track2()->Track()->TrackId() << " " +// << pair->Track2()->Track()->Flags() << " " +// << pair->Track2()->Track()->KinkIndex(0) << " " +// << pair->Track2()->Track()->KinkIndex(1) << " " +// << pair->Track2()->Track()->KinkIndex(2) << " " +// << endl; +// } - fShareNumerator->Fill(tQinv, hsfval); - fQualityNumerator->Fill(tQinv, hsmval); +// fShareNumerator->Fill(tQinv, hsfval); +// fQualityNumerator->Fill(tQinv, hsmval); // cout << "AliFemtoShareQualityCorrFctn::AddRealPair : " << pair->qInv() << " " << tQinv << //" " << pair->Track1().FourMomentum() << " " << pair->Track2().FourMomentum() << endl; } @@ -218,7 +255,7 @@ void AliFemtoShareQualityCorrFctn::AddMixedPair( AliFemtoPair* pair){ if (pair->Track1()->Track()->TPCclusters().TestBitNumber(imap) && pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) { // Do they share it ? - if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap) || + if (pair->Track1()->Track()->TPCsharing().TestBitNumber(imap) && pair->Track2()->Track()->TPCsharing().TestBitNumber(imap)) { // cout << "A shared cluster !!!" << endl; diff --git a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityPairCut.cxx b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityPairCut.cxx index 22d51007f2d..a25d59fd96b 100644 --- a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityPairCut.cxx +++ b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityPairCut.cxx @@ -1,3 +1,9 @@ +///////////////////////////////////////////////////////////////////////////// +// // +// AliFemtoShareQualityPairCut - a pair cut which checks for some pair // +// qualities that attempt to identify slit/doubly reconstructed tracks // +// // +///////////////////////////////////////////////////////////////////////////// /*************************************************************************** * * $Id$ @@ -25,8 +31,12 @@ ClassImp(AliFemtoShareQualityPairCut) AliFemtoShareQualityPairCut::AliFemtoShareQualityPairCut(): fNPairsPassed(0), fNPairsFailed(0), - fShareQualityMax(1.0) + fShareQualityMax(1.0), + fShareFractionMax(1.0), + fRemoveSameLabel(0) { + // Default constructor + // Nothing to do } //__________________ AliFemtoShareQualityPairCut::~AliFemtoShareQualityPairCut(){ @@ -34,6 +44,7 @@ AliFemtoShareQualityPairCut::~AliFemtoShareQualityPairCut(){ } //__________________ bool AliFemtoShareQualityPairCut::Pass(const AliFemtoPair* pair){ + // Check for pairs that are possibly shared/double reconstruction bool temp; Int_t nh = 0; @@ -79,26 +90,35 @@ bool AliFemtoShareQualityPairCut::Pass(const AliFemtoPair* pair){ hsfval = ns*1.0/nh; } // if (hsmval > -0.4) { - cout << "Pair quality: " << hsmval << " " << an << " " << nh << " " - << (pair->Track1()->Track()) << " " - << (pair->Track2()->Track()) << endl; - cout << "Bits: " << pair->Track1()->Track()->TPCclusters().GetNbits() << endl; +// cout << "Pair quality: " << hsmval << " " << an << " " << nh << " " +// << (pair->Track1()->Track()) << " " +// << (pair->Track2()->Track()) << endl; +// cout << "Bits: " << pair->Track1()->Track()->TPCclusters().GetNbits() << endl; // } - if (hsfval > 0.0) { - cout << "Pair sharity: " << hsfval << " " << ns << " " << nh << " " << hsmval << " " << an << " " << nh << endl; - } +// if (hsfval > 0.0) { +// cout << "Pair sharity: " << hsfval << " " << ns << " " << nh << " " << hsmval << " " << an << " " << nh << endl; +// } - temp = hsmval < fShareQualityMax; + temp = (hsmval < fShareQualityMax) && (hsfval < fShareFractionMax); + + if (fRemoveSameLabel) { + if (abs(pair->Track1()->Track()->Label()) == abs(pair->Track2()->Track()->Label())) { + cout << "Found a pair with same label " << pair->Track1()->Track()->Label() << endl; + cout << "Quality Sharity Passed " << hsmval << " " << hsfval << " " << pair->QInv() << " " << temp << endl; + temp = kFALSE; + } + } temp ? fNPairsPassed++ : fNPairsFailed++; return temp; } //__________________ AliFemtoString AliFemtoShareQualityPairCut::Report(){ - string Stemp = "AliFemtoShareQuality Pair Cut - remove shared and split pairs\n"; char Ctemp[100]; - sprintf(Ctemp,"Number of pairs which passed:\t%ld Number which failed:\t%ld\n",fNPairsPassed,fNPairsFailed); - Stemp += Ctemp; - AliFemtoString returnThis = Stemp; + // Prepare the report from the execution + string stemp = "AliFemtoShareQuality Pair Cut - remove shared and split pairs\n"; char ctemp[100]; + sprintf(ctemp,"Number of pairs which passed:\t%ld Number which failed:\t%ld\n",fNPairsPassed,fNPairsFailed); + stemp += ctemp; + AliFemtoString returnThis = stemp; return returnThis;} //__________________ @@ -106,17 +126,30 @@ void AliFemtoShareQualityPairCut::SetShareQualityMax(Double_t aShareQualityMax) fShareQualityMax = aShareQualityMax; } -Double_t AliFemtoShareQualityPairCut::GetAliFemtoShareQualityMax() { +Double_t AliFemtoShareQualityPairCut::GetAliFemtoShareQualityMax() const { return fShareQualityMax; } +void AliFemtoShareQualityPairCut::SetShareFractionMax(Double_t aShareFractionMax) { + fShareFractionMax = aShareFractionMax; +} +Double_t AliFemtoShareQualityPairCut::GetAliFemtoShareFractionMax() const { + return fShareFractionMax; +} + TList *AliFemtoShareQualityPairCut::ListSettings() { // return a list of settings in a writable form TList *tListSetttings = new TList(); char buf[200]; snprintf(buf, 200, "AliFemtoShareQualityPairCut.sharequalitymax=%lf", fShareQualityMax); + snprintf(buf, 200, "AliFemtoShareQualityPairCut.sharefractionmax=%lf", fShareFractionMax); tListSetttings->AddLast(new TObjString(buf)); return tListSetttings; } + +void AliFemtoShareQualityPairCut::SetRemoveSameLabel(Bool_t aRemove) +{ + fRemoveSameLabel = aRemove; +} diff --git a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityPairCut.h b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityPairCut.h index 33aa755626a..e355fe1e729 100644 --- a/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityPairCut.h +++ b/PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityPairCut.h @@ -1,3 +1,9 @@ +///////////////////////////////////////////////////////////////////////////// +// // +// AliFemtoShareQualityPairCut - a pair cut which checks for some pair // +// qualities that attempt to identify slit/doubly reconstructed tracks // +// // +///////////////////////////////////////////////////////////////////////////// /*************************************************************************** * * $Id$ @@ -14,8 +20,8 @@ **************************************************************************/ -#ifndef AliFemtoShareQualityPairCut_hh -#define AliFemtoShareQualityPairCut_hh +#ifndef ALIFEMTOSHAREQUALITYPAIRCUT_H +#define ALIFEMTOSHAREQUALITYPAIRCUT_H // do I need these lines ? //#ifndef StMaker_H @@ -27,20 +33,28 @@ class AliFemtoShareQualityPairCut : public AliFemtoPairCut{ public: AliFemtoShareQualityPairCut(); - AliFemtoShareQualityPairCut(const AliFemtoShareQualityPairCut&); - ~AliFemtoShareQualityPairCut(); - - virtual bool Pass(const AliFemtoPair*); + AliFemtoShareQualityPairCut(const AliFemtoShareQualityPairCut& cut); + virtual ~AliFemtoShareQualityPairCut(); + + virtual bool Pass(const AliFemtoPair* pair); virtual AliFemtoString Report(); virtual TList *ListSettings(); AliFemtoShareQualityPairCut* Clone(); void SetShareQualityMax(Double_t aAliFemtoShareQualityMax); - Double_t GetAliFemtoShareQualityMax(); + Double_t GetAliFemtoShareQualityMax() const; + void SetShareFractionMax(Double_t aAliFemtoShareFractionMax); + Double_t GetAliFemtoShareFractionMax() const; + void SetRemoveSameLabel(Bool_t aRemove); + + protected: + long fNPairsPassed; // Number of pairs consideered that passed the cut + long fNPairsFailed; // Number of pairs consideered that failed the cut + + private: + Double_t fShareQualityMax; // Maximum allowed pair quality + Double_t fShareFractionMax; // Maximum allowed share fraction + Bool_t fRemoveSameLabel; // If 1 pairs with two tracks with the same label will be removed -private: - long fNPairsPassed; - long fNPairsFailed; - Double_t fShareQualityMax; #ifdef __ROOT__ ClassDef(AliFemtoShareQualityPairCut, 0) @@ -51,7 +65,8 @@ inline AliFemtoShareQualityPairCut::AliFemtoShareQualityPairCut(const AliFemtoSh AliFemtoPairCut(c), fNPairsPassed(0), fNPairsFailed(0), - fShareQualityMax(1.0) // no cut + fShareQualityMax(1.0), + fShareFractionMax(1.0)// no cut { /* no-op */ } inline AliFemtoShareQualityPairCut* AliFemtoShareQualityPairCut::Clone() { AliFemtoShareQualityPairCut* c = new AliFemtoShareQualityPairCut(*this); return c;} -- 2.43.5