***************************************************************************
*
* $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
*
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),
fminTPCclsF(0),
fminITScls(0),
fNTracksPassed(0),
- fNTracksFailed(0)
+ fNTracksFailed(0),
+ fRemoveKinks(kFALSE),
+ fMostProbable(0)
{
// Default constructor
fNTracksPassed = fNTracksFailed = 0;
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"<<endl;
//cout<<fPidProbPion[0]<<" < pi ="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
//cout<<" status "<<track->Label()<<" "<<track->Flags()<<" "<<track->TPCnclsF()<<" "<<track->ITSncls()<<endl;
if ((track->Flags()&fStatus)!=fStatus)
{
- // cout<<track->Flags()<<" "<<fStatus<<" no go through status"<<endl;
+ // cout<<track->Flags()<<" "<<fStatus<<" no go through status"<<endl;
return false;
}
//cout<<fPt[0]<<" < Pt ="<<Pt<<" <"<<fPt[1]<<endl;
return false;
}
+// cout << "Track has pids: "
+// << track->PidProbElectron() << " "
+// << track->PidProbMuon() << " "
+// << track->PidProbPion() << " "
+// << track->PidProbKaon() << " "
+// << track->PidProbProton() << " "
+// << track->PidProbElectron()+track->PidProbMuon()+track->PidProbPion()+track->PidProbKaon()+track->PidProbProton() << endl;
+
+
if ((track->PidProbElectron()<fPidProbElectron[0])||(track->PidProbElectron()>fPidProbElectron[1]))
{
fNTracksFailed++;
//cout<<fPidProbMuon[0]<<" < mi="<<track->PidProbMuon()<<" <"<<fPidProbMuon[1]<<endl;
return false;
}
+ if (fRemoveKinks) {
+ if ((track->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"<<endl;
// cout<<fLabel<<" Label="<<track->Label()<<endl;
//------------------------------
AliFemtoString AliFemtoESDTrackCut::Report()
{
+ // Prepare report from the execution
string tStemp;
char tCtemp[100];
sprintf(tCtemp,"Particle mass:\t%E\n",this->Mass());
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);
+}
public:
AliFemtoESDTrackCut();
- //~AliFemtoESDTrackCut();
+ virtual ~AliFemtoESDTrackCut();
virtual bool Pass(const AliFemtoTrack* aTrack);
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);
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
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)
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
/*
*$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
*
#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"
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<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
- for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
- fClusterPerPadrow[tPad] = new list<Int_t>();
- }
- fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
- for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
- fSharedList[tPad] = new list<Int_t>();
- }
}
AliFemtoEventReaderESDKine::AliFemtoEventReaderESDKine(const AliFemtoEventReaderESDKine &aReader) :
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;
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<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
- for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
- fClusterPerPadrow[tPad] = new list<Int_t>();
- list<Int_t>::iterator iter;
- for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
- fClusterPerPadrow[tPad]->push_back(*iter);
- }
- }
- fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
- for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
- fSharedList[tPad] = new list<Int_t>();
- list<Int_t>::iterator iter;
- for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
- fSharedList[tPad]->push_back(*iter);
- }
- }
- for (unsigned int veciter = 0; veciter<aReader.fListOfFiles.size(); veciter++)
- fListOfFiles.push_back(aReader.fListOfFiles[veciter]);
+ fEvent = new AliESDEvent();
}
//__________________
//Destructor
//delete fListOfFiles;
delete fTree;
delete fEvent;
- delete fEsdFile;
if (fRunLoader) delete fRunLoader;
-
- for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
- fClusterPerPadrow[tPad]->clear();
- delete fClusterPerPadrow[tPad];
- }
- delete [] fClusterPerPadrow;
- for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
- fSharedList[tPad]->clear();
- delete fSharedList[tPad];
- }
- delete [] fSharedList;
}
//__________________
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; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
- fClusterPerPadrow[tPad]->clear();
- delete fClusterPerPadrow[tPad];
- }
- delete [] fClusterPerPadrow;
- }
-
- if (fSharedList) {
- for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
- fSharedList[tPad]->clear();
- delete fSharedList[tPad];
- }
- delete [] fSharedList;
- }
-
- fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
- for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
- fClusterPerPadrow[tPad] = new list<Int_t>();
- list<Int_t>::iterator iter;
- for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
- fClusterPerPadrow[tPad]->push_back(*iter);
- }
- }
- fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
- for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
- fSharedList[tPad] = new list<Int_t>();
- list<Int_t>::iterator iter;
- for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
- fSharedList[tPad]->push_back(*iter);
- }
- }
- for (unsigned int veciter = 0; veciter<aReader.fListOfFiles.size(); veciter++)
- fListOfFiles.push_back(aReader.fListOfFiles[veciter]);
-
return *this;
}
//__________________
-AliFemtoString AliFemtoEventReaderESDKine::Report()
+AliFemtoString AliFemtoEventReaderESDKine::Report() const
{
// create reader report
AliFemtoString temp = "\n This is the AliFemtoEventReaderESDKine\n";
fInputFile=string(inputFile);
cout<<"Input File set on "<<fInputFile<<endl;
ifstream infile(inputFile);
+
+ fTree = new TChain("esdTree");
+
if(infile.good()==true)
{
//checking if all give files have good tree inside
if (tree!=0x0)
{
cout<<"putting file "<<string(buffer)<<" into analysis"<<endl;
- fListOfFiles.push_back(string(buffer));
+ fTree->AddFile(buffer);
delete tree;
}
esdFile->Close();
}
}
-//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 "<<fFileName<<" "<<fCurFile<<endl;
-
- fCurFile++;
- return true;
-}
void AliFemtoEventReaderESDKine::SetConstrained(const bool constrained)
{
fConstrained=constrained;
// convert it to AliFemtoEvent and return
// for further analysis
AliFemtoEvent *hbtEvent = 0;
- TString tFriendFileName;
TString tGAliceFilename;
if (fCurEvent==fNumberofEvent)//open next file
{
- cout<<"next file"<<endl;
- if(GetNextFile())
- {
- delete fEventFriend;
- fEventFriend = 0;
- delete fEvent;//added 1.04.2007
- fEvent=new AliESD();
- // delete fTree;
- //fTree=0;
- delete fEsdFile;
+ if (fNumberofEvent == 0) {
+ fEvent=new AliESDEvent();
//ESD data
- fEsdFile=TFile::Open(fFileName.c_str(),"READ");
- fTree = (TTree*) fEsdFile->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 "<<fNumberofEvent<<endl;
- fCurEvent=0;
- // simulation data reading setup
- tGAliceFilename = fFileName.c_str();
- tGAliceFilename.ReplaceAll("AliESDs","galice");
- 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();
-
+// fTree->SetBranchStatus("fTracks.fCalibContainer",0);
+
+
+ fNumberofEvent=fTree->GetEntries();
+
+ if (fNumberofEvent == 0) {
+ cout<<"no event in input "<<endl;
+ fReaderStatus=1;
+ return hbtEvent;
}
+
+ cout<<"Number of Entries in the input "<<fNumberofEvent<<endl;
+ fCurEvent=0;
+ // simulation data reading setup
+
+ }
else //no more data to read
{
cout<<"no more files "<<hbtEvent<<endl;
}
cout<<"starting to read event "<<fCurEvent<<endl;
fTree->GetEvent(fCurEvent);//getting next event
- fEvent->SetESDfriend(fEventFriend);
// vector<int> 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();
nofTracks=fEvent->GetNumberOfTracks();
int realnofTracks=0;//number of track which we use ina analysis
- // Clear the shared cluster list
- for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
- fClusterPerPadrow[tPad]->clear();
- }
- for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
- fSharedList[tPad]->clear();
- }
-
-
- for (int i=0;i<nofTracks;i++) {
- const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
-
- list<Int_t>::iterator tClustIter;
-
- Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
- Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
- for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
- if (tTrackIndices[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;i<nofTracks;i++)
{
bool tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
delete trackCopy;
continue;
}
- const AliFmThreeVectorD origin(fV1[0],fV1[1],fV1[2]);
+ const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
//setting helix I do not if it is ok
- AliFmPhysicalHelixD helix(ktP,origin,(double)(fEvent->GetMagneticField())*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];
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<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
- if (tTrackIndices[tNcl] > 0)
- trackCopy->SetTPCcluster(tNcl, 1);
- else
- trackCopy->SetTPCcluster(tNcl, 0);
- }
-
- // Fill shared cluster information
- list<Int_t>::iterator tClustIter;
-
- for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
- if (tTrackIndices[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);
hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
fCurEvent++;
+ fCurRLEvent++;
cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
if (fCurEvent== fNumberofEvent)//if end of current file close all
{
fTree->Reset();
delete fTree;
- fEsdFile->Close();
}
return hbtEvent;
}
/*
*$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
*
*/
-#ifndef ALIFEMTOEVENTREADERESD_H
-#define ALIFEMTOEVENTREADERESD_H
+#ifndef ALIFEMTOEVENTREADERESDKINE_H
+#define ALIFEMTOEVENTREADERESDKINE_H
#include "AliFemtoEventReader.h"
#include "AliFemtoEnumeration.h"
#include <string>
#include <vector>
-#include "TTree.h"
-#include "AliESD.h"
-#include "AliESDfriend.h"
+#include "TChain.h"
+#include "AliESDEvent.h"
#include <list>
#include "AliRunLoader.h"
#include "AliFemtoModelHiddenInfo.h"
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);
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<string> 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<Int_t> **fSharedList; //! Table (one list per padrow) of clusters which are shared
- list<Int_t> **fClusterPerPadrow; //! Table (one list per padrow) of clusters in each padrow
#ifdef __ROOT__
ClassDef(AliFemtoEventReaderESDKine, 1)
/// ///
////////////////////////////////////////////////////////////////////////////////
#include "AliFemtoModelBPLCMSCorrFctn.h"
+#include "AliFemtoPair.h"
+#include "AliFemtoModelManager.h"
#include <cstdio>
#ifdef __ROOT__
//____________________________
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;
}
//_________________________
{
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);
}
//_________________________
-void AliFemtoModelBPLCMSCorrFctn::WriteOutHistos(){
-
- fNumerator->Write();
- fDenominator->Write();
+void AliFemtoModelBPLCMSCorrFctn::Write(){
+ AliFemtoModelCorrFctn::Write();
+ fNumerator3DTrue->Write();
+ fNumerator3DFake->Write();
+ fDenominator3D->Write();
fQinvHisto->Write();
}
//____________________________
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();
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;
+}
/// 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);
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
AliFemtoModelGausRinvFreezeOutGenerator();
AliFemtoModelGausRinvFreezeOutGenerator(const AliFemtoModelGausRinvFreezeOutGenerator &aModel);
virtual ~AliFemtoModelGausRinvFreezeOutGenerator();
- virtual void GenerateFreezeOut(AliFemtoPair *aPair);;
+ virtual void GenerateFreezeOut(AliFemtoPair *aPair);
void SetSizeInv(Double_t aSizeInv);
+/////////////////////////////////////////////////////////////////////////////
+// //
+// AliFemtoQPairCut - a simple cut which selects pairs based on the values //
+// of their respective q components /
+// //
+/////////////////////////////////////////////////////////////////////////////
/***************************************************************************
*
* $Id$
***************************************************************************
*
* $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
*
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;
//__________________
bool AliFemtoQPairCut::Pass(const AliFemtoPair* pair)
{
+ // Select pairs based on their q values
//bool temp = true;
//temp ? fNPairsPassed++ : fNPairsFailed++;
if ((fabs(pair->QLongCMS())<fQlong[0])||(fabs(pair->QLongCMS())>fQlong[1]))
//__________________
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;
}
//__________________
+/////////////////////////////////////////////////////////////////////////////
+// //
+// AliFemtoQPairCut - a simple cut which selects pairs based on the values //
+// of their respective q components /
+// //
+/////////////////////////////////////////////////////////////////////////////
+
/***************************************************************************
*
* $Id $
***************************************************************************
*
* $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
*
**************************************************************************/
-#ifndef AliFemtoQPairCut_hh
-#define AliFemtoQPairCut_hh
+#ifndef ALIFEMTOQPAIRCUT_H
+#define ALIFEMTOQPAIRCUT_H
// do I need these lines ?
//#ifndef StMaker_H
AliFemtoQPairCut();
~AliFemtoQPairCut();
- virtual bool Pass(const AliFemtoPair*);
+ virtual bool Pass(const AliFemtoPair* pair);
virtual AliFemtoString Report();
virtual TList *ListSettings();
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__
// 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"
//____________________________
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;
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;
nh++;
}
}
- if (tQinv < 0.01) {
- cout << "Qinv of the pair is " << tQinv << endl;
- cout << "Clusters: " << endl;
- for (unsigned int imap=0; imap<pair->Track1()->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; imap<pair->Track1()->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;
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; imap<pair->Track1()->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;
}
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;
+/////////////////////////////////////////////////////////////////////////////
+// //
+// AliFemtoShareQualityPairCut - a pair cut which checks for some pair //
+// qualities that attempt to identify slit/doubly reconstructed tracks //
+// //
+/////////////////////////////////////////////////////////////////////////////
/***************************************************************************
*
* $Id$
AliFemtoShareQualityPairCut::AliFemtoShareQualityPairCut():
fNPairsPassed(0),
fNPairsFailed(0),
- fShareQualityMax(1.0)
+ fShareQualityMax(1.0),
+ fShareFractionMax(1.0),
+ fRemoveSameLabel(0)
{
+ // Default constructor
+ // Nothing to do
}
//__________________
AliFemtoShareQualityPairCut::~AliFemtoShareQualityPairCut(){
}
//__________________
bool AliFemtoShareQualityPairCut::Pass(const AliFemtoPair* pair){
+ // Check for pairs that are possibly shared/double reconstruction
bool temp;
Int_t nh = 0;
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;}
//__________________
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;
+}
+/////////////////////////////////////////////////////////////////////////////
+// //
+// AliFemtoShareQualityPairCut - a pair cut which checks for some pair //
+// qualities that attempt to identify slit/doubly reconstructed tracks //
+// //
+/////////////////////////////////////////////////////////////////////////////
/***************************************************************************
*
* $Id$
**************************************************************************/
-#ifndef AliFemtoShareQualityPairCut_hh
-#define AliFemtoShareQualityPairCut_hh
+#ifndef ALIFEMTOSHAREQUALITYPAIRCUT_H
+#define ALIFEMTOSHAREQUALITYPAIRCUT_H
// do I need these lines ?
//#ifndef StMaker_H
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)
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;}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff