Anti-Splitting Cut implemented
authorskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 19 Jan 2004 09:30:11 +0000 (09:30 +0000)
committerskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 19 Jan 2004 09:30:11 +0000 (09:30 +0000)
13 files changed:
HBTAN/AliHBTAnalysis.cxx
HBTAN/AliHBTAnalysis.h
HBTAN/AliHBTClusterMap.cxx [new file with mode: 0644]
HBTAN/AliHBTClusterMap.h [new file with mode: 0644]
HBTAN/AliHBTPairCut.cxx
HBTAN/AliHBTPairCut.h
HBTAN/AliHBTParticle.cxx
HBTAN/AliHBTParticle.h
HBTAN/AliHBTReaderTPC.cxx
HBTAN/AliHBTReaderTPC.h
HBTAN/HBTAnalysisLinkDef.h
HBTAN/hbtanalysis.C
HBTAN/libHBTAN.pkg

index d68057d52dc0bfbed42ccd940bacbbb1d7333b97..80cab788d77474c222d5a3b7b3b0401fd0bc2472 100644 (file)
 #include "AliHBTRun.h"
 #include "AliHBTEvent.h"
 #include "AliHBTReader.h"
-#include "AliHBTParticle.h"
-#include "AliHBTParticleCut.h"
 #include "AliHBTPair.h"
-#include "AliHBTPairCut.h"
 #include "AliHBTFunction.h"
 #include "AliHBTMonitorFunction.h"
 #include "AliHBTEventBuffer.h"
@@ -61,12 +58,16 @@ AliHBTAnalysis::AliHBTAnalysis():
   fTrackMonitorFunctions ( new AliHBTMonOneParticleFctn* [fgkFctnArraySize]),    
   fParticleAndTrackMonitorFunctions ( new AliHBTMonTwoParticleFctn* [fgkFctnArraySize]),    
   fPairCut(new AliHBTEmptyPairCut()),//empty cut - accepts all particles
-  fAntiMergingCut(0x0),
   fBufferSize(2),
   fDisplayMixingInfo(fgkDefaultMixingInfo),
-  fIsOwner(kFALSE)
+  fIsOwner(kFALSE),
+  fPass(&AliHBTAnalysis::PassPartAndTrack), //by default perform cut on both track and particle pair 
+  fPass1(&AliHBTAnalysis::PassPartAndTrack1), //used onluy by ProcessTracksAndParticles
+  fPass2(&AliHBTAnalysis::PassPartAndTrack2),
+  fPassPairProp(&AliHBTAnalysis::PassPairPropPartAndTrack)
  {
    //default constructor
+   
  }
 /*************************************************************************************/ 
 
@@ -86,10 +87,13 @@ AliHBTAnalysis::AliHBTAnalysis(const AliHBTAnalysis& in):
   fTrackMonitorFunctions(0x0),
   fParticleAndTrackMonitorFunctions(0x0),
   fPairCut(0x0),
-  fAntiMergingCut(0x0),
   fBufferSize(fgkDefaultBufferSize),
   fDisplayMixingInfo(fgkDefaultMixingInfo),
-  fIsOwner(kFALSE)
+  fIsOwner(kFALSE),
+  fPass(&AliHBTAnalysis::PassPartAndTrack), //by default perform cut on both track and particle pair 
+  fPass1(&AliHBTAnalysis::PassPartAndTrack1),
+  fPass2(&AliHBTAnalysis::PassPartAndTrack2),
+  fPassPairProp(&AliHBTAnalysis::PassPairPropPartAndTrack)
  {
 //copy constructor
    Fatal("AliHBTAnalysis(const AliHBTAnalysis&)","Sensless");
@@ -119,7 +123,6 @@ AliHBTAnalysis::~AliHBTAnalysis()
    delete [] fParticleAndTrackMonitorFunctions; 
 
    delete fPairCut; // always have an copy of an object - we create we dstroy
-   delete fAntiMergingCut;
  }
 
 /*************************************************************************************/ 
@@ -269,7 +272,9 @@ void AliHBTAnalysis::ProcessTracksAndParticles()
 //In order to minimize calling AliRun::GetEvent (we need at one time particles from different events),
 //the loops are splited
   
-// cuta on particles only
+// cut on particles only -- why?
+// - PID: when we make resolution analysis we want to take only tracks with correct PID
+// We need cut on tracks because there are data characteristic to 
   
   AliHBTParticle * part1, * part2;
   AliHBTParticle * track1, * track2;
@@ -352,9 +357,9 @@ void AliHBTAnalysis::ProcessTracksAndParticles()
 //            return;
 //          }
          
-         Bool_t firstcut = fPairCut->GetFirstPartCut()->Pass(part1);
+         Bool_t firstcut = (this->*fPass1)(part1,track1);
          if (fBufferSize != 0) 
-           if ( (firstcut == kFALSE) || (fPairCut->GetSecondPartCut()->Pass(part1) == kFALSE) )
+           if ( (firstcut == kFALSE) || ( (this->*fPass2)(part1,track1) == kFALSE ) )
             {
               //accepted by any cut
               // we have to copy because reader keeps only one event
@@ -383,9 +388,9 @@ void AliHBTAnalysis::ProcessTracksAndParticles()
             track2= trackEvent->GetParticle(k);
             trackpair->SetParticles(track1,track2);
 
-            if(fPairCut->Pass(partpair) ) //check pair cut 
+            if( (this->*fPass)(partpair,trackpair) ) //check pair cut 
               { //do not meets crietria of the pair cut, try with swapped pairs
-                if( fPairCut->Pass(partpair->GetSwapedPair()) )
+                if( (this->*fPass)(partpair->GetSwapedPair(),trackpair->GetSwapedPair()) )
                   continue; //swaped pairs do not meet criteria of pair cut as well, take next particle
                 else 
                  { //swaped pair meets all the criteria
@@ -435,9 +440,9 @@ void AliHBTAnalysis::ProcessTracksAndParticles()
                 track2= trackEvent2->GetParticle(l);
                 trackpair->SetParticles(track1,track2);
 
-                if( fPairCut->Pass(partpair) ) //check pair cut
+                if( (this->*fPass)(partpair,trackpair) ) //check pair cut
                   { //do not meets crietria of the 
-                    if( fPairCut->Pass(partpair->GetSwapedPair()) )
+                    if( (this->*fPass)(partpair->GetSwapedPair(),trackpair->GetSwapedPair()) )
                       continue;
                     else 
                      {
@@ -451,10 +456,6 @@ void AliHBTAnalysis::ProcessTracksAndParticles()
                   tmppartpair = partpair;
                  }
 
-                //Anti merging cut is only on tracks makes the background
-                if (fAntiMergingCut) 
-                  if (fAntiMergingCut->Pass(trackpair)) continue;
-
                 for(ii = 0;ii<fNParticleFunctions;ii++)
                   fParticleFunctions[ii]->ProcessDiffEventParticles(tmppartpair);
                  
@@ -590,10 +591,6 @@ void AliHBTAnalysis::ProcessTracks()
                   tmptrackpair = trackpair;
                  }
                  
-                //Anti merging cut is only on tracks makes the background
-                if (fAntiMergingCut) 
-                  if (fAntiMergingCut->Pass(trackpair)) continue;
-
                 for(ii = 0;ii<fNTrackFunctions;ii++)
                   fTrackFunctions[ii]->ProcessDiffEventParticles(tmptrackpair);
                  
@@ -730,18 +727,6 @@ void AliHBTAnalysis::ProcessParticles()
       partEvent1 = partbuffer.Push(partEvent1); 
     }//while (fReader->Next() == kFALSE)
 }
-/*************************************************************************************/ 
-
-void AliHBTAnalysis::SetAntiMergingCut(Float_t x)
-{
-  if (x < 0)
-   {
-     Error("SetAntiMergingCut","Value less then 0");
-     return;
-   }
-  if (fAntiMergingCut == 0x0) fAntiMergingCut = new AliHBTAvSeparationCut(x);
-  else fAntiMergingCut->SetMinimum(x);
-}
 /*************************************************************************************/
 
 void AliHBTAnalysis::WriteFunctions()
@@ -1035,7 +1020,7 @@ void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
             track2= trackEvent2->GetParticle(k);
             trackpair->SetParticles(track1,track2);
 
-            if( (fPairCut->PassPairProp(trackpair)) ) //check pair cut
+            if( (this->*fPassPairProp)(partpair,trackpair) ) //check pair cut
              { //do not meets crietria of the pair cut
               continue; 
              }
@@ -1077,7 +1062,7 @@ void AliHBTAnalysis::ProcessTracksAndParticlesNonIdentAnal()
             track2= trackEvent3->GetParticle(k);
             trackpair->SetParticles(track1,track2);
 
-            if( (fPairCut->PassPairProp(trackpair)) ) //check pair cut
+            if( (this->*fPassPairProp)(partpair,trackpair) ) //check pair cut
              { //do not meets crietria of the pair cut
               continue; 
              }
@@ -1370,9 +1355,6 @@ void AliHBTAnalysis::FilterOut(AliHBTEvent* outpart1, AliHBTEvent* outpart2, Ali
   outtrack1->Reset();
   outtrack2->Reset();
   
-  AliHBTParticleCut *cut1 = fPairCut->GetFirstPartCut();
-  AliHBTParticleCut *cut2 = fPairCut->GetSecondPartCut();
-  
   Bool_t in1, in2;
   
   for (Int_t i = 0; i < inpart->GetNumberOfParticles(); i++)
@@ -1381,8 +1363,8 @@ void AliHBTAnalysis::FilterOut(AliHBTEvent* outpart1, AliHBTEvent* outpart2, Ali
      part = inpart->GetParticle(i);
      track = intrack->GetParticle(i);
      
-     if ( (cut1->Pass(track))  ) in1 = kFALSE; //if part  is rejected by cut1, in1 is false
-     if ( (cut2->Pass(track))  ) in2 = kFALSE; //if part  is rejected by cut2, in2 is false
+     if ( ((this->*fPass1)(part,track))  ) in1 = kFALSE; //if part  is rejected by cut1, in1 is false
+     if ( ((this->*fPass2)(part,track))  ) in2 = kFALSE; //if part  is rejected by cut2, in2 is false
      
      if (gDebug)//to be removed in real analysis     
      if ( in1 && in2 ) //both cuts accepted, should never happen, just in case
@@ -1479,7 +1461,44 @@ Bool_t AliHBTAnalysis::IsNonIdentAnalysis()
  
  return kTRUE;
 }
+/*************************************************************************************/ 
 
+void AliHBTAnalysis::SetCutsOnParticles()
+{
+ // -- aplies only to Process("TracksAndParticles")
+ // (ProcessTracksAndParticles and ProcessTracksAndParticlesNonIdentAnal)
+ // Only particles properties are checkes against cuts
+  fPass = &AliHBTAnalysis::PassPart;
+  fPass1 = &AliHBTAnalysis::PassPart1;
+  fPass2 = &AliHBTAnalysis::PassPart2;
+  fPassPairProp = &AliHBTAnalysis::PassPairPropPart;
+  
+}
+/*************************************************************************************/ 
+
+void AliHBTAnalysis::SetCutsOnTracks()
+{
+ // -- aplies only to Process("TracksAndParticles")
+ // (ProcessTracksAndParticles and ProcessTracksAndParticlesNonIdentAnal)
+ // Only tracks properties are checkes against cuts
+  fPass = &AliHBTAnalysis::PassTrack;
+  fPass1 = &AliHBTAnalysis::PassTrack1;
+  fPass2 = &AliHBTAnalysis::PassTrack2;
+  fPassPairProp = &AliHBTAnalysis::PassPairPropTrack;
+}
+/*************************************************************************************/ 
+
+void AliHBTAnalysis::SetCutsOnTracksAndParticles()
+{
+ // -- aplies only to Process("TracksAndParticles")
+ // (ProcessTracksAndParticles and ProcessTracksAndParticlesNonIdentAnal)
+ // Both, tracks and particles, properties are checked against cuts
+  fPass = &AliHBTAnalysis::PassPartAndTrack;
+  fPass1 = &AliHBTAnalysis::PassPartAndTrack1;
+  fPass2 = &AliHBTAnalysis::PassPartAndTrack2;
+  fPassPairProp = &AliHBTAnalysis::PassPairPropPartAndTrack;
+}
 /*************************************************************************************/ 
 
 void AliHBTAnalysis::PressAnyKey()
@@ -1495,4 +1514,5 @@ void AliHBTAnalysis::PressAnyKey()
    }
 }
 
+/*************************************************************************************/ 
 
index aac3848ace13bcd5df1696609678c95d9d29fc9a..fd0a93d8a867c25387cdbb7900c1f232d34eec5d 100644 (file)
 //_________________________________________________________
 
 #include <TObject.h>
+#include "AliHBTParticle.h"
+#include "AliHBTPairCut.h"
+#include "AliHBTParticleCut.h"
 
-class AliHBTParticleCut;
 class AliHBTCut;
-class AliHBTPairCut;
 class AliHBTPair;
 
 class AliHBTRun;
@@ -35,7 +36,6 @@ class AliHBTTwoPairFctn;
 
 class AliHBTMonOneParticleFctn;
 class AliHBTMonTwoParticleFctn;
-class AliHBTAvSeparationCut;
 
 class TList;
 
@@ -50,8 +50,6 @@ class AliHBTAnalysis: public TObject
      virtual void Process(Option_t* option = "TracksAndParticles");
      
      void SetGlobalPairCut(AliHBTPairCut* cut);
-     void SetAntiMergingCut(AliHBTAvSeparationCut *am){fAntiMergingCut = am;}
-     void SetAntiMergingCut(Float_t x);
      
      void AddTrackFunction(AliHBTOnePairFctn* f);
      void AddParticleFunction(AliHBTOnePairFctn* f);
@@ -75,6 +73,10 @@ class AliHBTAnalysis: public TObject
      void   ResetFunctions();
      void   SetDisplayInfo(Int_t howoften){fDisplayMixingInfo = howoften;}//defines every each line info about mixing is displayed
      
+     void   SetCutsOnParticles(); // -- aplies only to Process Tracks And Particles
+     void   SetCutsOnTracks();// -- aplies only to Process Tracks And Particles
+     void   SetCutsOnTracksAndParticles();// Default // -- aplies only to Process Tracks And Particles
+     
      static void PressAnyKey();//small utility function that helps to make comfortable macros
    protected:
      
@@ -118,17 +120,52 @@ class AliHBTAnalysis: public TObject
      /**********************************************/
 
      AliHBTPairCut*          fPairCut;//! Pair cut applied for all mixed particles
-     AliHBTAvSeparationCut*  fAntiMergingCut;//Anti-Splitting cut (only denominator)
       
      Int_t  fBufferSize; //!defines the size of buffer for mixed events; -1==MIX All
      Int_t  fDisplayMixingInfo;//!defines every which particle mixing info is displayed
      Bool_t fIsOwner;//!defines of all functions are supposed to be deleted while by the way of analysis defaulr false
 
    private:
+     Bool_t (AliHBTAnalysis::*fPass)(AliHBTPair* partpair, AliHBTPair* trackpair) const;//Pointer to function that performes pair cut
+     Bool_t (AliHBTAnalysis::*fPass1)(AliHBTParticle* partpair, AliHBTParticle* trackpair) const;//Pointer to function that performes cut on first particle
+     Bool_t (AliHBTAnalysis::*fPass2)(AliHBTParticle* partpair, AliHBTParticle* trackpair) const;//Pointer to function that performes cut on second particle
+     Bool_t (AliHBTAnalysis::*fPassPairProp)(AliHBTPair* partpair, AliHBTPair* trackpair) const;//Pointer to function that performes pair cut
+     
+     Bool_t PassPartAndTrack (AliHBTPair* partpair, AliHBTPair* trackpair) const {return (fPairCut->Pass(partpair))?kTRUE:fPairCut->Pass(trackpair);}
+     Bool_t PassPartAndTrack1(AliHBTParticle* part, AliHBTParticle* track) const;
+     Bool_t PassPartAndTrack2(AliHBTParticle* part, AliHBTParticle* track) const;
+     Bool_t PassPairPropPartAndTrack (AliHBTPair* partpair, AliHBTPair* trackpair) const {return (fPairCut->PassPairProp(partpair))?kTRUE:fPairCut->PassPairProp(trackpair);}
+     
+     Bool_t PassPart (AliHBTPair* partpair, AliHBTPair* /*trackpair*/) const{return fPairCut->Pass(partpair);}
+     Bool_t PassPart1(AliHBTParticle* part, AliHBTParticle* /*track*/) const{return fPairCut->GetFirstPartCut()->Pass(part);}
+     Bool_t PassPart2(AliHBTParticle* part, AliHBTParticle* /*track*/) const{return fPairCut->GetSecondPartCut()->Pass(part);}
+     Bool_t PassPairPropPart (AliHBTPair* partpair, AliHBTPair* /*trackpair*/) const{return fPairCut->PassPairProp(partpair);}
+     
+     Bool_t PassTrack (AliHBTPair* /*partpair*/, AliHBTPair* trackpair) const{return fPairCut->Pass(trackpair);}
+     Bool_t PassTrack1(AliHBTParticle* /*part*/, AliHBTParticle* track) const{return fPairCut->GetFirstPartCut()->Pass(track);}
+     Bool_t PassTrack2(AliHBTParticle* /*part*/, AliHBTParticle* track) const{return fPairCut->GetSecondPartCut()->Pass(track);}
+     Bool_t PassPairPropTrack (AliHBTPair* /*partpair*/, AliHBTPair* trackpair) const{return fPairCut->PassPairProp(trackpair);}
+
      static const UInt_t fgkFctnArraySize;//!
      static const UInt_t fgkDefaultMixingInfo;//!
      static const Int_t  fgkDefaultBufferSize;//!
 
      ClassDef(AliHBTAnalysis,0)
  };
+inline Bool_t AliHBTAnalysis::PassPartAndTrack1(AliHBTParticle* part,AliHBTParticle* track) const
+{
+//Checks first particle from both, particle and track pairs
+  AliHBTParticleCut* pc = fPairCut->GetFirstPartCut();
+  return (pc->Pass(part))?kTRUE:pc->Pass(track);
+}
+/*************************************************************************************/ 
+inline Bool_t AliHBTAnalysis::PassPartAndTrack2(AliHBTParticle* part,AliHBTParticle* track) const
+{
+//Checks second particle from both, particle and track pairs
+  AliHBTParticleCut* pc = fPairCut->GetSecondPartCut();
+  return (pc->Pass(part))?kTRUE:pc->Pass(track);
+}
+/*************************************************************************************/ 
 #endif
diff --git a/HBTAN/AliHBTClusterMap.cxx b/HBTAN/AliHBTClusterMap.cxx
new file mode 100644 (file)
index 0000000..ee4073f
--- /dev/null
@@ -0,0 +1,186 @@
+#include "AliHBTClusterMap.h"
+
+//_________________________________________________
+///////////////////////////////////////////////////
+//
+// class AliHBTClusterMap
+//
+// class that describes cluster occupation at TPC
+// Each padraw has a corresponding bit in fPadRawMap
+// 
+//
+// more info: http://aliweb.cern.ch/people/skowron/analyzer/index.html
+// Piotr.Skowronski@cern.ch
+//
+///////////////////////////////////////////////////////////////////////////
+
+
+#include "AliESDtrack.h"
+#include "AliTPCtrack.h"
+#include "AliHBTParticle.h"
+#include <TString.h>
+const Int_t AliHBTClusterMap::fNPadRows = 159;
+
+AliHBTClusterMap::AliHBTClusterMap():
+ fPadRawMap(fNPadRows)
+{
+//ctor
+
+}
+/***********************************************************************/
+AliHBTClusterMap::AliHBTClusterMap(AliESDtrack* track):
+ fPadRawMap(fNPadRows)
+{
+ //cotor
+ track->Print();//to shut up compiler warning
+} 
+/***********************************************************************/
+
+AliHBTClusterMap::AliHBTClusterMap(AliTPCtrack* track):
+ fPadRawMap(fNPadRows)
+{
+ //cotor
+ //Does not work since indeces in the claster index array 
+ //in the TPC track does not correspond to the padraw segmatation
+ if (AliHBTParticle::GetDebug() > 9) 
+   Info("AliHBTClusterMap",
+      "#####################################################################"); 
+ if (track == 0x0)
+  {
+    Error("AliHBTClusterMap","Pointer to TPC track is NULL");
+    return;
+  }
+ Int_t prevrow = -1;
+ Int_t i = 0;
+ for ( ; i < track->GetNumberOfClusters(); i++)
+  {
+    Int_t idx = track->GetClusterIndex(i);
+    Int_t sect = (idx&0xff000000)>>24;
+    Int_t row = (idx&0x00ff0000)>>16;
+    if (sect > 18) row +=63; //if it is outer sector, add number of inner sectors
+    if (AliHBTParticle::GetDebug() > 9)  
+      Info("AliHBTClusterMap","Cl.idx is %d, sect %d, row %d",idx,sect,row);
+      
+    fPadRawMap.SetBitNumber(row,kTRUE);
+    
+    //Fill the gap between previous row and this row with 0 bits
+    if (prevrow < 0) 
+     {
+       prevrow = row;//if previous bit was not assigned yet == this is the first one
+     }
+    else
+     { //we don't know the order (inner to outer or reverse)
+       //just to be save in case it is going change
+       Int_t n = 0, m = 0;
+       if (prevrow < row)
+        {
+          n = prevrow;
+          m = row;
+        }
+       else
+        {
+          n = row;
+          m = prevrow;
+        }
+       for (Int_t j = n+1; j < m; j++)
+        {
+          fPadRawMap.SetBitNumber(j,kFALSE);
+        }
+       prevrow = row; 
+     }
+  }
+  if (AliHBTParticle::GetDebug() > 2)
+   {
+     Print();
+   } 
+}
+/***********************************************************************/
+
+void AliHBTClusterMap::Print() const
+{
+//Prints the bit map 
+  TString msg;
+  for ( Int_t i = 0; i < fNPadRows; i++)
+   {
+     if ( fPadRawMap.TestBitNumber(i) )
+      {
+        msg+="1";
+      }
+     else
+      {
+        msg+="0";
+      }
+   }
+  Info("AliHBTClusterMap","BitMap is\n  %s",msg.Data());
+  
+}
+
+/***********************************************************************/
+
+Float_t AliHBTClusterMap::GetOverlapFactor(const AliHBTClusterMap& clmap) const
+{
+  //Returns quality factor FQ = Sum(An)/Sum(clusters)
+  //      | -1; if both tracks have a cluster on padrow n
+  //An = <  0; if neither track has a cluster on padrow n
+  //     |  1; if only one trackhas a cluster on padrow n
+  // Returned value ranges between 
+  //  -0.5 (low probability that these tracks are a split track)
+  //  and
+  //   1.0 (high probability that these tracks are a split track)
+  TString msg1;
+  TString msg2;
+  
+  Int_t nh = 0;
+  Int_t an = 0;
+  for ( Int_t i = 0; i < fNPadRows; i++)
+   {
+     Bool_t x = HasClAtPadRow(i);
+     Bool_t y = clmap.HasClAtPadRow(i);
+     
+     if (x) msg1+="1";else msg1+="0";
+     if (y) msg2+="1";else msg2+="0";
+     
+     if (x && y)//both have clasters
+      {
+       an--;
+       nh+=2;
+      }
+     else 
+      {
+       
+       if (x || y)//only one have cluters
+        {
+          an++;
+          nh++;
+        }
+      }
+   }
+  
+  
+  Float_t retval = 0.0;
+  if (nh > 0) retval = ((Float_t)an)/((Float_t)nh);
+  else Warning("GetOverlapFactor","Number of counted cluters is 0.");
+  
+  if (AliHBTParticle::GetDebug() > 2)
+   {
+     Info("GetOverlapFactor","Splitting Quality Factor is %f. SumAn = %d, SumClusters %d",retval,an,nh); 
+     if (retval == 1.0) 
+      { 
+        Print();
+        Info("AliHBTClusterMap","BitMap is\n  %s\n",msg1.Data());
+        clmap.Print(); 
+        Info("AliHBTClusterMap","BitMap is\n  %s\n\n\n\n",msg2.Data());
+      }
+     if (retval == -.5) 
+      { 
+        Print();
+        Info("AliHBTClusterMap","BitMap is\n  %s\n",msg1.Data());
+        clmap.Print(); 
+        Info("AliHBTClusterMap","BitMap is\n  %s\n\n\n\n",msg2.Data());
+      }
+   } 
+  return retval;
+}
diff --git a/HBTAN/AliHBTClusterMap.h b/HBTAN/AliHBTClusterMap.h
new file mode 100644 (file)
index 0000000..046c90e
--- /dev/null
@@ -0,0 +1,42 @@
+#ifndef ALIHBTCLUSTERMAP_H
+#define ALIHBTCLUSTERMAP_H
+//_________________________________________________
+///////////////////////////////////////////////////
+//
+// class AliHBTClusterMap
+//
+// class that describes cluster occupation at TPC
+// Each padraw has a corresponding bit in fPadRawMap
+// 
+//
+// more info: http://aliweb.cern.ch/people/skowron/analyzer/index.html
+// Piotr.Skowronski@cern.ch
+//
+///////////////////////////////////////////////////////////////////////////
+
+
+#include <TObject.h>
+#include <TBits.h>
+
+class AliTPCtrack;
+class AliESDtrack;
+class TBits;
+
+class AliHBTClusterMap: public TObject
+{
+  public:
+   AliHBTClusterMap();
+   AliHBTClusterMap(AliTPCtrack* track);
+   AliHBTClusterMap(AliESDtrack* track);
+   virtual ~AliHBTClusterMap(){}
+   Float_t GetOverlapFactor(const AliHBTClusterMap& clmap) const;
+   Bool_t  HasClAtPadRow(Int_t i) const { return fPadRawMap.TestBitNumber(i);}
+   void    Print() const;
+  protected:
+  private:
+   TBits    fPadRawMap;//bit vector of length 150 correspondind to total number of padraws in TPC
+   static const Int_t fNPadRows;
+   ClassDef(AliHBTClusterMap,1)
+};
+
+#endif 
index e97513143048921101af5fa2a16083fc94522848..20afdb394ca620a602775a10eaa5b76aa16854bf 100644 (file)
@@ -13,6 +13,7 @@
 #include "AliHBTPair.h"
 #include "AliHBTParticleCut.h"
 #include "AliHBTTrackPoints.h"
+#include "AliHBTClusterMap.h"
 
 ClassImp(AliHBTPairCut)
 const Int_t AliHBTPairCut::fgkMaxCuts = 50;
@@ -237,6 +238,21 @@ void AliHBTPairCut::SetAvSeparationRange(Double_t min, Double_t max)
 }
 /**********************************************************/
 
+void AliHBTPairCut::SetClusterOverlapRange(Double_t min,Double_t max)
+{
+  //sets cluster overlap factor cut ->Anti-Splitting cut
+  //cluster overlap factor ranges between 
+  // -0.5 (in all padrows both tracks have cluters) 
+  // and 1 (in all padrows one track has cluter and second has not)
+  // When Overlap Factor is 1 this pair of tracks in highly probable to be
+  // splitted track: one particle that is recontructed twise
+  // STAR uses range from -0.5 to 0.6 
+  
+  AliHbtBasePairCut* cut= FindCut(kHbtPairCutPropClOverlap);
+  if(cut) cut->SetRange(min,max);
+  else fCuts[fNCuts++] = new AliHBTCluterOverlapCut(min,max);
+}
+
 AliHbtBasePairCut* AliHBTPairCut::FindCut(AliHBTPairCutProperty property)
 {
   // Find the cut corresponding to "property"
@@ -316,21 +332,45 @@ ClassImp(AliHBTAvSeparationCut)
 Double_t AliHBTAvSeparationCut::GetValue(AliHBTPair* pair) const 
 {
   //chacks if avarage distance of two tracks is in given range
-  Warning("Pass","Checking Av Separation."); 
   AliHBTTrackPoints* tpts1 = pair->Particle1()->GetTrackPoints();
   if ( tpts1 == 0x0)
    {
-     Warning("Pass","Track 1 does not have Track Points. Pair NOT Passed.");
+     Warning("GetValue","Track 1 does not have Track Points. Pair NOT Passed.");
      return -10e5;
    }
 
   AliHBTTrackPoints* tpts2 = pair->Particle2()->GetTrackPoints();
   if ( tpts2 == 0x0)
    {
-     Warning("Pass","Track 2 does not have Track Points. Pair NOT Passed.");
+     Warning("GetValue","Track 2 does not have Track Points. Pair NOT Passed.");
      return -10e5;
    }
    
   return tpts1->AvarageDistance(*tpts2);
 }
 
+ClassImp(AliHBTCluterOverlapCut)
+
+Double_t  AliHBTCluterOverlapCut::GetValue(AliHBTPair* pair) const
+{
+  //Returns Cluter Overlap Factor
+  //It ranges between -0.5 (in all padrows both tracks have cluters) 
+  // and 1 (in all padrows one track has cluter and second has not)
+  // When Overlap Factor is 1 this pair of tracks in highly probable to be
+  // splitted track: one particle that is recontructed twise
+
+  AliHBTClusterMap* cm1 = pair->Particle1()->GetClusterMap();
+  if ( cm1 == 0x0)
+   {
+     Warning("GetValue","Track 1 does not have Cluster Map. Returning -0.5.");
+     return -.5;
+   }
+
+  AliHBTClusterMap* cm2 = pair->Particle2()->GetClusterMap();
+  if ( cm2 == 0x0)
+   {
+     Warning("GetValue","Track 2 does not have Cluster Map. Returning -0.5.");
+     return -.5;
+   }
+  return cm1->GetOverlapFactor(*cm2);
+}
index 69e86f273afd6f5375a7b3f66cf1c3fb92b272d3..ba53bde9c67e6ed950b2ec4e5f8540e95688ac2c 100644 (file)
@@ -22,6 +22,7 @@ enum AliHBTPairCutProperty
   kHbtPairCutPropQOutCMSLC,
   kHbtPairCutPropQLongCMSLC,
   kHbtPairCutPropAvSepar,
+  kHbtPairCutPropClOverlap,
   kHbtPairCutPropNone
 };
 
@@ -50,7 +51,8 @@ class AliHBTPairCut: public TNamed
   void SetQOutCMSLRange(Double_t min, Double_t max);
   void SetQSideCMSLRange(Double_t min, Double_t max);
   void SetQLongCMSLRange(Double_t min, Double_t max);
-  void SetAvSeparationRange(Double_t min,Double_t max = 10e5);
+  void SetAvSeparationRange(Double_t min,Double_t max = 10e5);//Anti-Merging Cut
+  void SetClusterOverlapRange(Double_t min,Double_t max);//Anti-Splitting Max range -0.5 1.0
       
   AliHBTParticleCut* GetFirstPartCut() const {return fFirstPartCut;}
   AliHBTParticleCut* GetSecondPartCut() const {return fSecondPartCut;}
@@ -226,5 +228,17 @@ class AliHBTAvSeparationCut: public AliHbtBasePairCut
   virtual Double_t  GetValue(AliHBTPair* pair) const;
   ClassDef(AliHBTAvSeparationCut,1)
 };
+
+class AliHBTCluterOverlapCut: public AliHbtBasePairCut
+{
+ public:
+  AliHBTCluterOverlapCut(Double_t min = 0.0, Double_t max = 1e5):
+    AliHbtBasePairCut(min,max,kHbtPairCutPropClOverlap){}
+  virtual ~AliHBTCluterOverlapCut(){}
+  
+ protected:
+  virtual Double_t  GetValue(AliHBTPair* pair) const;
+  ClassDef(AliHBTCluterOverlapCut,1)
+};
   
 #endif
index 6524d8c091492b2e31efe201c3e42ad34578ab50..23e8a2108e8930cdfdbcb1c3012eee7173756c54 100644 (file)
@@ -13,6 +13,7 @@
 #include <TParticle.h>
 #include <TClass.h>
 #include "AliHBTTrackPoints.h"
+#include "AliHBTClusterMap.h"
 
 ClassImp(AliHBTParticle)
 
@@ -21,7 +22,7 @@ Int_t AliHBTParticle::fgDebug = 0;
 AliHBTParticle::AliHBTParticle():  
  fPdgIdx(0), fIdxInEvent(0),fNPids(0),fPids(0x0),fPidProb(0x0),
  fCalcMass(0),fPx(0), fPy(0),fPz(0),fE(0), fVx(0), fVy(0),fVz(0),fVt(0),
- fTrackPoints(0x0)
+ fTrackPoints(0x0),fClusterMap(0x0)
 {//empty particle
 }
 //______________________________________________________________________________
@@ -33,7 +34,7 @@ AliHBTParticle::AliHBTParticle(Int_t pdg, Int_t idx,
   fCalcMass(0), 
   fPx(px), fPy(py),fPz(pz),fE(etot), 
   fVx(vx), fVy(vy),fVz(vz),fVt(time),
-  fTrackPoints(0x0)
+  fTrackPoints(0x0),fClusterMap(0x0)
 {
 //mormal constructor
   SetPdgCode(pdg);
@@ -54,7 +55,7 @@ AliHBTParticle::AliHBTParticle(Int_t pdg, Float_t prob, Int_t idx,
   fCalcMass(0), 
   fPx(px), fPy(py),fPz(pz),fE(etot), 
   fVx(vx), fVy(vy),fVz(vz),fVt(time),
-  fTrackPoints(0x0)
+  fTrackPoints(0x0),fClusterMap(0x0)
 {
 //mormal constructor
   SetPdgCode(pdg,prob);
@@ -74,7 +75,7 @@ AliHBTParticle::AliHBTParticle(const AliHBTParticle& in):
    fCalcMass(in.GetCalcMass()),
    fPx(in.Px()),fPy(in.Py()),fPz(in.Pz()),fE(in.Energy()), 
    fVx(in.Vx()),fVy(in.Vy()),fVz(in.Vz()),fVt(in.T()),
-   fTrackPoints(0x0) 
+   fTrackPoints(0x0), fClusterMap(0x0)
 {
  //Copy constructor
  for(Int_t i = 0; i<fNPids; i++)
@@ -83,12 +84,10 @@ AliHBTParticle::AliHBTParticle(const AliHBTParticle& in):
     fPidProb[i] = in.fPidProb[i];
   }
  
- //fTrackPoints = (in.fTrackPoints)?(AliHBTTrackPoints*)in.fTrackPoints->Clone():0x0;
-  
  if (in.fTrackPoints)
    fTrackPoints = (AliHBTTrackPoints*)in.fTrackPoints->Clone();
-   
-  
+ if (in.fClusterMap)
+   fClusterMap = (AliHBTClusterMap*)in.fClusterMap->Clone();
 }
 
 //______________________________________________________________________________
@@ -98,7 +97,7 @@ AliHBTParticle::AliHBTParticle(const TParticle &p,Int_t idx):
    fCalcMass(p.GetCalcMass()),
    fPx(p.Px()),fPy(p.Py()),fPz(p.Pz()),fE(p.Energy()), 
    fVx(p.Vx()),fVy(p.Vy()),fVz(p.Vz()),fVt(p.T()),
-   fTrackPoints(0x0)
+   fTrackPoints(0x0),fClusterMap(0x0)
 {
  //all copied in the initialization
  SetPdgCode(p.GetPdgCode());
@@ -111,6 +110,7 @@ AliHBTParticle::~AliHBTParticle()
   delete [] fPids;
   delete [] fPidProb;
   delete fTrackPoints;
+  delete fClusterMap;
 }
 //______________________________________________________________________________
 
@@ -144,6 +144,9 @@ AliHBTParticle& AliHBTParticle::operator=(const AliHBTParticle& in)
   delete fTrackPoints;
   fTrackPoints = (in.fTrackPoints)?(AliHBTTrackPoints*)fTrackPoints->Clone():0x0;
   
+  delete fClusterMap;
+  fClusterMap =  (in.fClusterMap)?(AliHBTClusterMap*)in.fClusterMap->Clone():0x0;
+  
   return *this;
 }
 //______________________________________________________________________________
index 9904aca2d041e65960298c9f3d2363b5e01aad32..9a4d6cee57871942009e6009654159a26b663c82 100644 (file)
@@ -21,6 +21,7 @@
 
 class TParticle;
 class AliHBTTrackPoints;
+class AliHBTClusterMap;
 
 class AliHBTParticle : public TObject
 {
@@ -115,7 +116,9 @@ public:
   
   void           SetTrackPoints(AliHBTTrackPoints* tpts){fTrackPoints = tpts;}
   AliHBTTrackPoints* GetTrackPoints() const {return fTrackPoints;}
-      
+  void           SetClusterMap(AliHBTClusterMap* cm){fClusterMap = cm;}
+  AliHBTClusterMap* GetClusterMap() const {return fClusterMap;}
+  
   static void    SetDebug(Int_t dbg=1){fgDebug=dbg;}
   static Int_t   GetDebug(){return fgDebug;}
   static Int_t   fgDebug; //debug printout level
@@ -143,6 +146,8 @@ private:
   Double_t       fVt;                   // t of production vertex
 
   AliHBTTrackPoints* fTrackPoints;      // track positions along trajectory - used by anti-merging cut
+  AliHBTClusterMap*  fClusterMap;       // bit map of cluters occupation; 1 if has cluter on given layer/padrow/...
+    
   ClassDef(AliHBTParticle,3)  // TParticle vertex particle information
 };
 
index 8a643511b09b3ee6fe3a919b24956c06c817ae8b..d37498309101e0deb8e2acacd00cdaf98837492c 100644 (file)
@@ -31,7 +31,7 @@
 #include "AliHBTParticle.h"
 #include "AliHBTParticleCut.h"
 #include "AliHBTTrackPoints.h"
-
+#include "AliHBTClusterMap.h"
 
 
 ClassImp(AliHBTReaderTPC)
@@ -44,6 +44,7 @@ AliHBTReaderTPC::AliHBTReaderTPC():
  fUseMagFFromRun(kTRUE),
  fNTrackPoints(0),
  fdR(0.0),
+ fClusterMap(kFALSE),
  fNClustMin(0),
  fNClustMax(150),
  fNChi2PerClustMin(0.0),
@@ -101,6 +102,7 @@ AliHBTReaderTPC::AliHBTReaderTPC(TObjArray* dirs, const Char_t* galicefilename):
  fUseMagFFromRun(kTRUE),
  fNTrackPoints(0),
  fdR(0.0),
+ fClusterMap(kFALSE),
  fNClustMin(0),
  fNClustMax(150),
  fNChi2PerClustMin(0.0),
@@ -285,7 +287,12 @@ Int_t AliHBTReaderTPC::ReadNext()
           AliHBTTrackPoints* tpts = new AliHBTTrackPoints(fNTrackPoints,iotrack,fdR);
           track->SetTrackPoints(tpts);
         }
-        
+       if (  fClusterMap ) 
+        {
+          AliHBTClusterMap* cmap = new AliHBTClusterMap(iotrack);
+          track->SetClusterMap(cmap);
+        }
+    
        fParticlesEvent->AddParticle(part);
        fTracksEvent->AddParticle(track);
       }
index e61c3adf1c318a7fe41ff053a31481d8a0b078e3..2d2ac50ac47630efa095775fb45eba1b19850c78 100644 (file)
@@ -43,8 +43,9 @@ class AliHBTReaderTPC: public AliHBTReader
     void          SetC44Range(Float_t min, Float_t max);
     void          SetNumberOfTrackPoints(Int_t n = 5,Float_t dr = 30.0) {fNTrackPoints = n; fdR = dr;}
     Int_t         GetNumberOfTrackPoints() const {return fNTrackPoints;}
+    void          SetClusterMap(Bool_t flag = kTRUE){fClusterMap = flag;}
+
   protected:
-    //in the future this class is will read global tracking
     Int_t         ReadNext();
     Int_t         OpenNextSession();
     void          DoOpenError(const char* msgfmt, ...);
@@ -56,8 +57,14 @@ class AliHBTReaderTPC: public AliHBTReader
     Bool_t        fUseMagFFromRun;//flag indicating if using field specified in gAlice (kTRUE)
                                // or enforece other defined by fMagneticField
     
-    Int_t         fNTrackPoints;//number of track points
+    Int_t         fNTrackPoints;//number of track points; if==0 track points are not created
     Float_t       fdR;//spacing between points (along radius) in cm
+                      //Track Points are needed for Anti-Merging Cut
+    
+    Bool_t        fClusterMap;//Flag indicating if Claster Map should be created for each track
+                              //Claster map is needed for Anti-Splitting Cut
+
+    //Cut Parameters specific to TPC tracks
         
     Int_t         fNClustMin;//Number of clusters min value
     Int_t         fNClustMax;//Number of clusters max value
@@ -83,6 +90,7 @@ class AliHBTReaderTPC: public AliHBTReader
   private:
     
     Bool_t CheckTrack(AliTPCtrack* t);
+
   public:
     ClassDef(AliHBTReaderTPC,3)
 };
index 0fad6c972383feaabb241f250c851b4f5d0c4048..b2724a112602b9260e1d403c3c8e3d4d76263591 100644 (file)
@@ -6,6 +6,7 @@
  
 #pragma link C++ class AliHBTAnalysis+;
 #pragma link C++ class AliHBTParticle+;
+#pragma link C++ class AliHBTPair+;
 #pragma link C++ class AliHBTEvent+;
 #pragma link C++ class AliHBTRun+;
 #pragma link C++ class AliHBTEventBuffer+;
@@ -34,7 +35,6 @@
 #pragma link C++ class AliHBTMonTwoParticleFctn2D+;
 #pragma link C++ class AliHBTMonTwoParticleFctn3D+;
 
-#pragma link C++ class AliHBTPair+;
 #pragma link C++ class AliHBTParticleCut-;
 #pragma link C++ class AliHBTEmptyParticleCut-;
 #pragma link C++ class AliHBTPairCut-;
@@ -49,6 +49,7 @@
 #pragma link C++ class AliHBTQOutCMSLCCut+;
 #pragma link C++ class AliHBTQLongCMSLCCut+;
 #pragma link C++ class AliHBTAvSeparationCut+;
+#pragma link C++ class AliHBTCluterOverlapCut+;
     
 #pragma link C++ class AliHBTMomentumCut+;
 #pragma link C++ class AliHBTPtCut+;
@@ -77,6 +78,7 @@
 #pragma link C++ class AliHBTReaderInternal+;
 
 #pragma link C++ class AliHBTTrackPoints+;
+#pragma link C++ class AliHBTClusterMap+;
     
 #pragma link C++ class AliHBTQInvCorrelFctn+;
 #pragma link C++ class AliHBTTwoKStarCorrelFctn+;
index 14ae5ddd7cdf4b4ba22f649c85d9562abac1fb45..7fcc2e439172a5972cee6cdc4d3cab8e8b2ff955 100644 (file)
@@ -25,7 +25,9 @@ void hbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles",
                 char *outfile = "hbtanalysis.root")
  {
    
-  AliHBTTrackPoints::SetDebug(1);
+//  AliHBTTrackPoints::SetDebug(1);
+  AliHBTParticle::SetDebug(0);
+  
 //HBT Anlysis Macro
 //Anlyzes TPC recontructed tracks and simulated particles that corresponds to them
 
@@ -68,6 +70,7 @@ void hbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles",
   /***********************************************************/
   //Create analysis and reader
   AliHBTAnalysis * analysis = new AliHBTAnalysis();
+  analysis->SetCutsOnTracks();
   
   AliHBTReader* reader;
   Int_t kine = strcmp(datatype,"Kine");
@@ -91,7 +94,8 @@ void hbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles",
   else if(!TPC)
    {
     reader = new AliHBTReaderTPC();
-    ((AliHBTReaderTPC*)reader)->SetNumberOfTrackPoints(5,30.);
+    //((AliHBTReaderTPC*)reader)->SetNumberOfTrackPoints(5,30.);
+    ((AliHBTReaderTPC*)reader)->SetClusterMap();
    }
   else if(!ITSv1)
    {
@@ -146,7 +150,8 @@ void hbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles",
   Float_t qinvmax = 0.05;//50MeV
   paircut->SetQInvRange(qinvmin,qinvmax);  
   
- // paircut->SetAvSeparationRange(10.,10000.);  
+ // paircut->SetAvSeparationRange(10.); //AntiMerging
+  paircut->SetClusterOverlapRange(-.5,1.0);//AntiSplitting
   
 //  AliHBTParticleCut* partcut= new AliHBTParticleCut();
 //  partcut->SetPID(kPiPlus);
@@ -155,7 +160,6 @@ void hbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles",
 //  paircut->SetSecondPartCut(partcut);
   
   analysis->SetGlobalPairCut(paircut);
-  analysis->SetAntiMergingCut(5.);//in cm
   
   /***********************************************************/    
   /*****   W E I G H T S        ******************************/    
index b5c394014b50d5f69e6339492e5393491b4eab16..4ba5422b912866466771c43583f148da024f7d03 100644 (file)
@@ -11,10 +11,10 @@ AliHBTMonitorFunction.cxx       AliHBTTwoTrackEffFctn.cxx \
 AliHBTQResolutionFctns.cxx      AliHBTQDistributionFctns.cxx \
 AliHBTMonDistributionFctns.cxx  AliHBTMonResolutionFctns.cxx \
 AliHBTLLWeights.cxx             AliHBTWeightFctn.cxx \
-AliHBTWeightsPID.cxx          AliHBTWeightTheorFctn.cxx \
-AliHBTPositionRandomizer.cxx AliHBTEventBuffer.cxx \
+AliHBTWeightsPID.cxx            AliHBTWeightTheorFctn.cxx \
+AliHBTPositionRandomizer.cxx    AliHBTEventBuffer.cxx \
 AliHBTCorrFitFctn.cxx AliHBTCorrectCorrelFctn.cxx \
-AliHBTTrackPoints.cxx
+AliHBTTrackPoints.cxx AliHBTClusterMap.cxx
 FSRCS   = fsiini.F  fsiw.F  led_bldata.F  ltran12.F
 
 HDRS:= $(SRCS:.cxx=.h)