Port of changes from v4-07-Release and additional rule conformance
authorakisiel <akisiel@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 2 Nov 2007 09:37:15 +0000 (09:37 +0000)
committerakisiel <akisiel@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 2 Nov 2007 09:37:15 +0000 (09:37 +0000)
12 files changed:
PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoESDTrackCut.cxx
PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoESDTrackCut.h
PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoEventReaderESDKine.cxx
PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoEventReaderESDKine.h
PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelBPLCMSCorrFctn.cxx
PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelBPLCMSCorrFctn.h
PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelGausRinvFreezeOutGenerator.h
PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoQPairCut.cxx
PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoQPairCut.h
PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityCorrFctn.cxx
PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityPairCut.cxx
PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityPairCut.h

index 438976d..741f008 100644 (file)
@@ -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
  *
 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"<<endl;
   //cout<<fPidProbPion[0]<<" < pi ="<<track->PidProbPion()<<" <"<<fPidProbPion[1]<<endl;
@@ -90,7 +127,7 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
       //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;
        }
        
@@ -146,6 +183,15 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
       //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++;
@@ -181,6 +227,23 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
       //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;
@@ -200,6 +263,7 @@ bool AliFemtoESDTrackCut::Pass(const AliFemtoTrack* track)
 //------------------------------
 AliFemtoString AliFemtoESDTrackCut::Report()
 {
+  // Prepare report from the execution
   string tStemp;
   char tCtemp[100];
   sprintf(tCtemp,"Particle mass:\t%E\n",this->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);  
+}
index 4322b13..a90e5cf 100644 (file)
@@ -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
 
index 1714402..761795d 100644 (file)
@@ -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
  *
 #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<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) :
@@ -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<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
@@ -141,19 +103,7 @@ AliFemtoEventReaderESDKine::~AliFemtoEventReaderESDKine()
   //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;
 }
 
 //__________________
@@ -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; 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";
@@ -234,6 +145,9 @@ void AliFemtoEventReaderESDKine::SetInputFile(const char* inputFile)
   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
@@ -248,7 +162,7 @@ void AliFemtoEventReaderESDKine::SetInputFile(const char* inputFile)
              if (tree!=0x0)
                {
                  cout<<"putting file  "<<string(buffer)<<" into analysis"<<endl;
-                 fListOfFiles.push_back(string(buffer));
+                 fTree->AddFile(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 "<<fFileName<<" "<<fCurFile<<endl;
-
-  fCurFile++;
-  return true;
-}
 void AliFemtoEventReaderESDKine::SetConstrained(const bool constrained)
 {
   fConstrained=constrained;
@@ -286,64 +188,49 @@ AliFemtoEvent* AliFemtoEventReaderESDKine::ReturnHbtEvent()
   // 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;
@@ -353,10 +240,32 @@ AliFemtoEvent* AliFemtoEventReaderESDKine::ReturnHbtEvent()
     }          
   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();
        
@@ -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; 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
@@ -451,14 +330,14 @@ AliFemtoEvent* AliFemtoEventReaderESDKine::ReturnHbtEvent()
        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];
@@ -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<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);
@@ -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 "<<nofTracks<<" real number "<<realnofTracks<<endl;
   if (fCurEvent== fNumberofEvent)//if end of current file close all
     {   
       fTree->Reset(); 
       delete fTree;
-      fEsdFile->Close();
     }
   return hbtEvent; 
 }
index 013470e..5a8acf7 100644 (file)
@@ -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
  *
  */
   
 
-#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"
@@ -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<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)
index 0c68d6d..628f9bf 100644 (file)
@@ -7,6 +7,8 @@
 ///                                                                          ///
 ////////////////////////////////////////////////////////////////////////////////
 #include "AliFemtoModelBPLCMSCorrFctn.h"
+#include "AliFemtoPair.h"
+#include "AliFemtoModelManager.h"
 #include <cstdio>
 
 #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;
+}
index 56fdb4d..47a5bd6 100644 (file)
@@ -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
 
index 224f68d..281a380 100644 (file)
@@ -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);
   
index 22919bb..8854c7f 100644 (file)
@@ -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())<fQlong[0])||(fabs(pair->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;
 }
 //__________________
index 90fe35a..5400876 100644 (file)
@@ -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__
index c9dc4f7..0effcc8 100644 (file)
@@ -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; 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;
@@ -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; 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;
 }
@@ -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;
index 22d5100..a25d59f 100644 (file)
@@ -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;
+}
index 33aa755..e355fe1 100644 (file)
@@ -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
 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;}