Separation Cut in Pixels implemented
authorskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 13 Jun 2004 16:41:28 +0000 (16:41 +0000)
committerskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 13 Jun 2004 16:41:28 +0000 (16:41 +0000)
HBTAN/AliHBTPairCut.cxx
HBTAN/AliHBTPairCut.h
HBTAN/AliHBTParticle.cxx
HBTAN/AliHBTParticle.h
HBTAN/AliHBTReaderESD.cxx
HBTAN/AliHBTReaderESD.h
HBTAN/AliHBTTrackPoints.cxx
HBTAN/AliHBTTrackPoints.h
HBTAN/HBTAnalysisLinkDef.h

index c873783b31085f65c59476d6605d7de3aa31cc16..ebd1925d4b469c0c44249d3d81375907fc51454b 100644 (file)
@@ -131,6 +131,16 @@ Bool_t AliHBTPairCut::PassPairProp(AliHBTPair* pair) const
 }
 /**********************************************************/
 
+void AliHBTPairCut::Print()
+{
+ //Prints the cut
+  for (Int_t i = 0;i<fNCuts;i++)
+    {
+      fCuts[i]->Dump();
+    }
+}
+/**********************************************************/
+
 void AliHBTPairCut::SetFirstPartCut(AliHBTParticleCut* cut)
 {
   // set cut for the first particle
@@ -229,6 +239,7 @@ void AliHBTPairCut::SetKStarRange(Double_t min, Double_t max)
   else fCuts[fNCuts++] = new AliHBTKStarCut(min,max);
 }
 /**********************************************************/
+
 void AliHBTPairCut::SetAvSeparationRange(Double_t min, Double_t max)
 {
   //sets avarage separation cut ->Anti-Merging cut
@@ -238,6 +249,15 @@ void AliHBTPairCut::SetAvSeparationRange(Double_t min, Double_t max)
 }
 /**********************************************************/
 
+void AliHBTPairCut::SetPixelSeparation(Double_t drphi, Double_t dz)
+{
+  //Anti-Merging Cut for first pixel layer
+  AliHbtBasePairCut* cut= FindCut(kHbtPairCutPropPixelSepar);
+  if(cut) cut->SetRange(drphi,dz);//In this cut fMin is drphi, and fMax dz
+  else fCuts[fNCuts++] = new AliHBTPixelSeparationCut(drphi,dz);
+}
+/**********************************************************/
+
 void AliHBTPairCut::SetClusterOverlapRange(Double_t min,Double_t max)
 {
   //sets cluster overlap factor cut ->Anti-Splitting cut
@@ -354,6 +374,54 @@ Double_t AliHBTAvSeparationCut::GetValue(AliHBTPair* pair) const
 }
 /******************************************************************/
 
+ClassImp(AliHBTPixelSeparationCut)
+
+Bool_t AliHBTPixelSeparationCut::Pass(AliHBTPair* pair) const
+{
+ //Checks if two tracks do not cross first pixels too close to each other
+ //If two tracks use the same cluster in pixels they are given
+ //the same position what skews theta angles (both are the same)
+ //These guys create artificial correlation in non-id analyses
+ //which is positive for identical polar angles (Qlong=0) 
+ //and negative for a little bit different theta angle (Qlong=epsilon)
+ //Such tracks "attracks" each other.
+  AliHBTTrackPoints* tpts1 = pair->Particle1()->GetITSTrackPoints();
+  if ( tpts1 == 0x0)
+   {//it could be simulated pair
+     Warning("Pass","Track 1 does not have ITS Track Points. Pair NOT Passed.");
+     return kTRUE;//reject 
+   }
+
+  AliHBTTrackPoints* tpts2 = pair->Particle2()->GetITSTrackPoints();
+  if ( tpts2 == 0x0)
+   {
+     Warning("Pass","Track 2 does not have ITS Track Points. Pair NOT Passed.");
+     return kTRUE;//reject 
+   }
+  Float_t  x1=0.0,y1=0.0,z1=0.0,x2=0.0,y2=0.0,z2=0.0;
+  tpts1->PositionAt(0,x1,y1,z1);
+  tpts2->PositionAt(0,x2,y2,z2);
+  
+//  Info("Pass","rphi %f z %f",fMin,fMax);
+//  Info("Pass","P1: %f %f %f", x1,y1,z1);
+//  Info("Pass","P2: %f %f %f", x2,y2,z2);
+  
+  Double_t dz = TMath::Abs(z1-z2);
+  
+  //fMax encodes treshold valaue of distance in Z
+  if (dz > fMax) return kFALSE;//pair accepted
+  
+  Double_t drphi = TMath::Hypot(x1-x2,y1-y2);
+  
+  //fMin encodes treshold valaue of distance in r-phi
+  if (drphi > fMin) return kFALSE;
+  
+  Info("Pass","Rejected !!!!!");
+  return kTRUE;//they are too close, rejected
+}
+/******************************************************************/
+
 ClassImp(AliHBTCluterOverlapCut)
 
 Double_t  AliHBTCluterOverlapCut::GetValue(AliHBTPair* pair) const
index 565f35f7690ca3a0b52c71aed0a61b8888821180..907115585950dab49094b7998c1c251515ee6b1e 100644 (file)
@@ -27,6 +27,7 @@ enum AliHBTPairCutProperty
   kHbtPairCutPropDeltaPt,
   kHbtPairCutPropAvSepar,
   kHbtPairCutPropClOverlap,
+  kHbtPairCutPropPixelSepar,
   kHbtPairCutPropNone
 };
 /******************************************************************/
@@ -50,6 +51,8 @@ class AliHBTPairCut: public TNamed
   
   virtual void AddBasePairCut(AliHbtBasePairCut* cut);
   
+  virtual void Print();
+  
   void SetQInvRange(Double_t min, Double_t max);
   void SetKtRange(Double_t min, Double_t max);
   void SetKStarRange(Double_t min, Double_t max);
@@ -57,6 +60,7 @@ class AliHBTPairCut: public TNamed
   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);//Anti-Merging Cut
+  void SetPixelSeparation(Double_t drphi=0.01,Double_t dz = 0.08);//Anti-Merging Cut for first pixel layer
   void SetClusterOverlapRange(Double_t min,Double_t max);//Anti-Splitting Max range -0.5 1.0
       
   AliHBTParticleCut* GetFirstPartCut() const {return fFirstPartCut;}
@@ -281,6 +285,21 @@ class AliHBTAvSeparationCut: public AliHbtBasePairCut
   ClassDef(AliHBTAvSeparationCut,1)
 };
 /******************************************************************/
+  
+class AliHBTPixelSeparationCut: public AliHbtBasePairCut
+{
+//Anti merging cut for the first layer of pixels
+ public:
+  AliHBTPixelSeparationCut(Double_t deltarphi = 0.01, Double_t deltaz = 0.08):
+    AliHbtBasePairCut(deltarphi,deltaz,kHbtPairCutPropPixelSepar){}
+  virtual ~AliHBTPixelSeparationCut(){}
+  Bool_t   Pass(AliHBTPair* pair) const;
+
+ protected:
+  virtual Double_t  GetValue(AliHBTPair* /*pair*/) const {return 0.0;}//not used
+  ClassDef(AliHBTPixelSeparationCut,1)
+};
+/******************************************************************/
 
 class AliHBTOutSideSameSignCut: public AliHbtBasePairCut
 {
index eeced89f53dd225c010f3aa1baea97a7dbd254b7..2760b5cc0406b7b4f0cc6a9780c0a47c97a0aa24 100644 (file)
@@ -25,7 +25,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),fClusterMap(0x0)
+ fTrackPoints(0x0),fITSTrackPoints(0x0),fClusterMap(0x0)
 {//empty particle
 }
 //______________________________________________________________________________
@@ -37,7 +37,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),fClusterMap(0x0)
+  fTrackPoints(0x0),fITSTrackPoints(0x0),fClusterMap(0x0)
 {
 //mormal constructor
   SetPdgCode(pdg);
@@ -58,7 +58,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),fClusterMap(0x0)
+  fTrackPoints(0x0),fITSTrackPoints(0x0),fClusterMap(0x0)
 {
 //mormal constructor
   SetPdgCode(pdg,prob);
@@ -78,7 +78,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), fClusterMap(0x0)
+   fTrackPoints(0x0), fITSTrackPoints(0x0), fClusterMap(0x0)
 {
  //Copy constructor
  for(Int_t i = 0; i<fNPids; i++)
@@ -89,6 +89,8 @@ AliHBTParticle::AliHBTParticle(const AliHBTParticle& in):
  
  if (in.fTrackPoints)
    fTrackPoints = (AliHBTTrackPoints*)in.fTrackPoints->Clone();
+ if (in.fITSTrackPoints)
+   fITSTrackPoints = (AliHBTTrackPoints*)in.fITSTrackPoints->Clone();
  if (in.fClusterMap)
    fClusterMap = (AliHBTClusterMap*)in.fClusterMap->Clone();
 }
@@ -100,7 +102,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),fClusterMap(0x0)
+   fTrackPoints(0x0), fITSTrackPoints(0x0), fClusterMap(0x0)
 {
  //all copied in the initialization
  SetPdgCode(p.GetPdgCode());
@@ -113,6 +115,7 @@ AliHBTParticle::~AliHBTParticle()
   delete [] fPids;
   delete [] fPidProb;
   delete fTrackPoints;
+  delete fITSTrackPoints;
   delete fClusterMap;
 }
 //______________________________________________________________________________
@@ -145,7 +148,10 @@ AliHBTParticle& AliHBTParticle::operator=(const AliHBTParticle& in)
   fVt = in.T();
   
   delete fTrackPoints;
-  fTrackPoints = (in.fTrackPoints)?(AliHBTTrackPoints*)fTrackPoints->Clone():0x0;
+  fTrackPoints = (in.fTrackPoints)?(AliHBTTrackPoints*)in.fTrackPoints->Clone():0x0;
+
+  delete fITSTrackPoints;
+  fITSTrackPoints = (in.fTrackPoints)?(AliHBTTrackPoints*)in.fITSTrackPoints->Clone():0x0;
   
   delete fClusterMap;
   fClusterMap =  (in.fClusterMap)?(AliHBTClusterMap*)in.fClusterMap->Clone():0x0;
index cd126ea3bd78b05d6dd24624a769b42c66fa7d79..4c6e3d2639a2e24e3da7b82d825194dbde202e89 100644 (file)
@@ -118,6 +118,10 @@ public:
   
   void           SetTrackPoints(AliHBTTrackPoints* tpts){fTrackPoints = tpts;}
   AliHBTTrackPoints* GetTrackPoints() const {return fTrackPoints;}
+
+  void           SetITSTrackPoints(AliHBTTrackPoints* tpts){fITSTrackPoints = tpts;}
+  AliHBTTrackPoints* GetITSTrackPoints() const {return fITSTrackPoints;}
+  
   void           SetClusterMap(AliHBTClusterMap* cm){fClusterMap = cm;}
   AliHBTClusterMap* GetClusterMap() const {return fClusterMap;}
   
@@ -146,7 +150,9 @@ private:
   Double_t       fVz;                   // z of production vertex
   Double_t       fVt;                   // t of production vertex
 
-  AliHBTTrackPoints* fTrackPoints;      // track positions along trajectory - used by anti-merging cut
+  AliHBTTrackPoints* fTrackPoints;      // track positions along trajectory - used by anti-merging cut 
+  AliHBTTrackPoints* fITSTrackPoints;   // track position at first pixels
+  
   AliHBTClusterMap*  fClusterMap;       // bit map of cluters occupation; 1 if has cluter on given layer/padrow/...
     
   static Int_t   fgDebug; //debug printout level
index 717a5426093d257bf4c416d6e0cc1a10b6292231..53988462feb01f2bc5a07fc6b153ca20fea13a5c 100644 (file)
@@ -10,6 +10,8 @@
 //                                                                  //
 //////////////////////////////////////////////////////////////////////
 
+#include <Riostream.h>
+
 #include <TPDGCode.h>
 #include <TString.h>
 #include <TObjString.h>
@@ -23,6 +25,7 @@
 #include <AliRunLoader.h>
 #include <AliStack.h>
 #include <AliESDtrack.h>
+#include <AliKalmanTrack.h>
 #include <AliESD.h>
 
 #include "AliHBTRun.h"
@@ -45,6 +48,7 @@ AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename, const Char_t* galfil
  fNTrackPoints(0),
  fdR(0.0),
  fClusterMap(kFALSE),
+ fITSTrackPoints(kFALSE),
  fMustTPC(kFALSE),
  fNTPCClustMin(0),
  fNTPCClustMax(150),
@@ -92,6 +96,7 @@ AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename, cons
  fNTrackPoints(0),
  fdR(0.0),
  fClusterMap(kFALSE),
+ fITSTrackPoints(kFALSE),
  fMustTPC(kFALSE),
  fNTPCClustMin(0),
  fNTPCClustMax(150),
@@ -257,12 +262,17 @@ Int_t AliHBTReaderESD::ReadESD(AliESD* esd)
    
   Float_t mf = esd->GetMagneticField(); 
 
-  if ( (mf == 0.0) && (fNTrackPoints > 0) )
+  if ( (mf == 0.0) && ((fNTrackPoints > 0) || fITSTrackPoints) )
    {
       Error("ReadESD","Magnetic Field is 0 and Track Points Demended. Skipping to next event.");
       return 1;
    }
-
+  if (fITSTrackPoints)
+   {
+     Info("ReadESD","Magnetic Field is %f",mf/10.);
+     AliKalmanTrack::SetMagneticField(mf/10.);
+   }
+   
   AliStack* stack = 0x0;
   if (fReadParticles && fRunLoader)
    {
@@ -410,6 +420,13 @@ Int_t AliHBTReaderESD::ReadESD(AliESD* esd)
        {
          tpts = new AliHBTTrackPoints(fNTrackPoints,esdtrack,mf,fdR);
        }
+      
+      AliHBTTrackPoints* itstpts = 0x0;
+      if (fITSTrackPoints) 
+       {
+         itstpts = new AliHBTTrackPoints(AliHBTTrackPoints::kITS,esdtrack);
+       }
+      
 
       AliHBTClusterMap* cmap = 0x0; 
       if (  fClusterMap ) 
@@ -464,6 +481,11 @@ Int_t AliHBTReaderESD::ReadESD(AliESD* esd)
            track->SetTrackPoints(tpts); 
          }
 
+        if (itstpts)
+         {
+           track->SetITSTrackPoints(itstpts); 
+         }
+
         if (cmap) 
          { 
            track->SetClusterMap(cmap);
@@ -486,8 +508,26 @@ Int_t AliHBTReaderESD::ReadESD(AliESD* esd)
       {
         delete particle;//particle was not stored in event
         delete tpts;
+        delete itstpts;
         delete cmap;
       }
+     else
+      {
+        if (particle->P() < 0.00001)
+         {
+           Info("ReadNext","###################################");
+           Info("ReadNext","###################################");
+           Info("ReadNext","Track Label %d",esdtrack->GetLabel());
+           TParticle *p = stack->Particle(esdtrack->GetLabel());
+           Info("ReadNext","");
+           p->Print();
+           Info("ReadNext","");
+           particle->Print();
+           TString command("touch BadInput");
+           ofstream sfile("BadEvent",ios::out);
+           sfile<<fCurrentDir<<endl;
+         }
+      } 
 
    }//for (Int_t i = 0;i<ntr; i++)  -- loop over tracks
 
index 5d5c0333c2d638e52f89e7129394fe87c884999d..30731c3364e349063aad884d81f2db3dc6b77af7 100644 (file)
@@ -50,6 +50,7 @@ class AliHBTReaderESD: public AliHBTReader
     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;}
+    void          SetITSTrackPoints(Bool_t flag = kTRUE){fITSTrackPoints = flag;}
     void          MustTPC(Bool_t flag){fMustTPC = flag;}
     
     enum ESpecies {kESDElectron = 0, kESDMuon, kESDPion, kESDKaon, kESDProton, kNSpecies};
@@ -76,6 +77,10 @@ class AliHBTReaderESD: public AliHBTReader
     
     Bool_t        fClusterMap;//Flag indicating if Claster Map should be created for each track
                               //Claster map is needed for Anti-Splitting Cut
+    
+    Bool_t        fITSTrackPoints;//Flag indicalting if track positions in ITS are to be read
+                                  //currently we use only position at first pixels wich are
+                     //used by anti-merging cut in non-id analysis
 
     Bool_t        fMustTPC;// must be reconstructed in TPC -> reject tracks reconstructed ITS stand alone
     
index 673aa8bb156812cbeb1a51333b9b6ba5daa5164d..7fe5222146d499f4b7085efeacd62c3de4243032 100644 (file)
@@ -1,3 +1,4 @@
+#include "AliHBTTrackPoints.h"
 //_________________________________
 ////////////////////////////////////////////////////////////
 //                                                        //
@@ -17,9 +18,9 @@
 #include <TMath.h>
 
 #include "AliESDtrack.h"
-#include "AliHBTTrackPoints.h"
 #include "AliTPCtrack.h"
 #include "AliTrackReference.h"
+#include "AliITStrackV2.h"
 
 ClassImp(AliHBTTrackPoints)
 
@@ -34,6 +35,32 @@ AliHBTTrackPoints::AliHBTTrackPoints():
 }
 /***************************************************************/
 
+AliHBTTrackPoints::AliHBTTrackPoints(AliHBTTrackPoints::ETypes type, AliESDtrack* track):
+ fN(0),
+ fX(0x0),
+ fY(0x0),
+ fZ(0x0)
+{
+  //constructor
+  switch (type)
+   {
+     case kITS:
+       //Up to now we keep only point in pixels
+       //Used only in non-id analysis
+       fN = 1;
+       fX = new Float_t[fN];
+       fY = new Float_t[fN];
+       fZ = new Float_t[fN];
+       MakeITSPoints(track);
+       break;
+
+     default:
+       Info("AliHBTTrackPoints","Not recognized type");
+   }
+   
+}
+/***************************************************************/
+
 AliHBTTrackPoints::AliHBTTrackPoints(Int_t n, AliESDtrack* track, Float_t mf, Float_t dr, Float_t r0):
  fN(n),
  fX(new Float_t[fN]),
@@ -131,6 +158,16 @@ AliHBTTrackPoints::AliHBTTrackPoints(Int_t n, AliTPCtrack* track, Float_t dr, Fl
   Double_t c=track->GetC();
   MakePoints(dr,r0,x,par,c,alpha);
 }  
+/***************************************************************/
+
+AliHBTTrackPoints::~AliHBTTrackPoints()
+{
+  //destructor
+  delete [] fX;
+  delete [] fY;
+  delete [] fZ;
+}
+/***************************************************************/
 
 void AliHBTTrackPoints::MakePoints( Float_t dr, Float_t r0, Double_t x, Double_t* par, Double_t c, Double_t alpha)
 {
@@ -235,15 +272,23 @@ void AliHBTTrackPoints::MakePoints( Float_t dr, Float_t r0, Double_t x, Double_t
 }
 /***************************************************************/
 
-AliHBTTrackPoints::~AliHBTTrackPoints()
+void AliHBTTrackPoints::MakeITSPoints(AliESDtrack* track)
 {
-  //destructor
-  delete [] fX;
-  delete [] fY;
-  delete [] fZ;
+//Calculates points in ITS
+// z=R*Pz/Pt
+ AliITStrackV2 itstrack(*track,kTRUE);
+ Double_t x,y,z;
+ itstrack.GetGlobalXYZat(4.0,x,y,z);
+ fX[0] = x;
+ fY[0] = y;
+ fZ[0] = z;
+//  Info("MakeITSPoints","X %f Y %f Z %f R asked %f R obtained %f",
+//             fX[0],fY[0],fZ[0],4.0,TMath::Hypot(fX[0],fY[0]));
 }
-/***************************************************************/
 
+/***************************************************************/
 void AliHBTTrackPoints::PositionAt(Int_t n, Float_t &x,Float_t &y,Float_t &z)
 {
   //returns position at point n
index d0751a998a15e816ba89f2b8d284e9a4a9d54b2d..cd7450ba45f876a5293f2f23c17cfd5474c3d9e0 100644 (file)
@@ -21,9 +21,12 @@ class AliESDtrack;
 class AliHBTTrackPoints: public TObject
 {
   public:
+    enum ETypes{kITS = 1};
+
     AliHBTTrackPoints();
     AliHBTTrackPoints(Int_t n, AliTPCtrack* track, Float_t dr=30, Float_t r0 = 84.1); //min TPC R  = 84.1; max TPC R =  246.6cm, 
     AliHBTTrackPoints(Int_t n, AliESDtrack* track, Float_t mf, Float_t dr=30,Float_t r0 = 84.1); //min TPC R  = 84.1; max TPC R =  246.6cm, 
+    AliHBTTrackPoints(AliHBTTrackPoints::ETypes type, AliESDtrack* track);
     AliHBTTrackPoints(const AliHBTTrackPoints& in);
     
     virtual ~AliHBTTrackPoints();
@@ -35,8 +38,11 @@ class AliHBTTrackPoints: public TObject
     void  SetDebug(Int_t deblevel){fgDebug = deblevel;} 
     static void testtpc(Int_t entr);
     static void testesd(Int_t entr,const char* fname = "AliESDs.root");
+
   protected:
     void MakePoints( Float_t dr, Float_t r0, Double_t x, Double_t* par, Double_t c, Double_t alpha);
+    void MakeITSPoints(AliESDtrack* track);
+    
   private:
     Int_t    fN;//number of points
     Float_t* fX;//[fN]positions at x
index cb86324964bebcd0ca7266095c075b96ab69d6b8..7e27427f5835e3992b96d049aca3639d744caf72 100644 (file)
@@ -54,6 +54,7 @@
 #pragma link C++ class AliHBTDeltaPhiCut+;
 #pragma link C++ class AliHBTDeltaThetaCut+;
 #pragma link C++ class AliHBTAvSeparationCut+;
+#pragma link C++ class AliHBTPixelSeparationCut+;
 #pragma link C++ class AliHBTCluterOverlapCut+;
 #pragma link C++ class AliHBTOutSideSameSignCut+;
 #pragma link C++ class AliHBTOutSideDiffSignCut+;