1. Adding function to find cosmic track pairs
authormarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 21 Jun 2008 16:14:03 +0000 (16:14 +0000)
committermarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 21 Jun 2008 16:14:03 +0000 (16:14 +0000)
2. Prpagating the pair of tracks to the DCA 0,0
3. Selection function isPair

(Marian)

TPC/AliTPCcalibCosmic.cxx
TPC/AliTPCcalibCosmic.h

index 772f2a1..8b8d7d8 100644 (file)
@@ -44,12 +44,15 @@ ClassImp(AliTPCcalibCosmic)
 
 AliTPCcalibCosmic::AliTPCcalibCosmic() 
   :AliTPCcalibBase(),
-  fHistNTracks(0),
-  fClusters(0),
-  fModules(0),
-  fHistPt(0),
-  fPtResolution(0),
-  fDeDx(0)
+   fHistNTracks(0),
+   fClusters(0),
+   fModules(0),
+   fHistPt(0),
+   fPtResolution(0),
+   fDeDx(0),
+   fCutMaxD(5),        // maximal distance in rfi ditection
+   fCutTheta(0.03),    // maximal distan theta
+   fCutMinDir(-0.99)   // direction vector products
 {  
   AliInfo("Defualt Constructor");  
 }
@@ -57,12 +60,15 @@ AliTPCcalibCosmic::AliTPCcalibCosmic()
 
 AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title) 
   :AliTPCcalibBase(),
-  fHistNTracks(0),
-  fClusters(0),
-  fModules(0),
-  fHistPt(0),
-  fPtResolution(0),
-  fDeDx(0)
+   fHistNTracks(0),
+   fClusters(0),
+   fModules(0),
+   fHistPt(0),
+   fPtResolution(0),
+   fDeDx(0),
+   fCutMaxD(5),        // maximal distance in rfi ditection
+   fCutTheta(0.03),    // maximal distan theta
+   fCutMinDir(-0.99)   // direction vector products
 {  
   SetName(name);
   SetTitle(title);
@@ -75,7 +81,7 @@ AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
   fPtResolution = new TH1F("PtResolution","Pt resolution",100,-50,50);
   fDeDx = new TH2F("DeDx","dEdx",500,0.01,20.,500,0.,500);
   BinLogX(fDeDx);
-  AliInfo("Non Defualt Constructor");  
+  AliInfo("Non Default Constructor");  
 }
 
 AliTPCcalibCosmic::~AliTPCcalibCosmic(){
@@ -85,53 +91,57 @@ AliTPCcalibCosmic::~AliTPCcalibCosmic(){
 }
 
 
+
+
+
 void AliTPCcalibCosmic::Process(AliESDEvent *event) {
-    
+  //
+  //
+  //
   if (!event) {
     Printf("ERROR: ESD not available");
     return;
-  }
-  
+  }  
   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
   if (!ESDfriend) {
    Printf("ERROR: ESDfriend not available");
    return;
   }
+  FindPairs(event);
 
-  printf("Hallo world: Im here\n");
-  
-  Int_t n=event->GetNumberOfTracks(); 
-  fHistNTracks->Fill(n);
-  
+  if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
+  Int_t ntracks=event->GetNumberOfTracks(); 
+  fHistNTracks->Fill(ntracks);
+  TObjArray  tpcSeeds(ntracks);
+  if (ntracks==0) return;
+  //
   //track loop
-  for (Int_t i=0;i<n;++i) { 
+  //
+  for (Int_t i=0;i<ntracks;++i) { 
    AliESDtrack *track = event->GetTrack(i); 
-   fClusters->Fill(track->GetTPCNcls());
-   
+   fClusters->Fill(track->GetTPCNcls());   
    AliExternalTrackParam * trackIn = new AliExternalTrackParam(*track->GetInnerParam());
    
    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
    TObject *calibObject;
    AliTPCseed *seed = 0;
-   for (Int_t l=0;calibObject=friendTrack->GetCalibObject(l);++l) {
-    if (seed=dynamic_cast<AliTPCseed*>(calibObject)) break;
+   for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
+     if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
    }
-   if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0));
-
-    
+   if (seed) tpcSeeds.AddAt(seed,i);
+   if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0)); 
   }
-  
-  // dE/dx,pt and ACORDE study --> studies which need the pair selection  
-  if (n > 2) return;
-  
-  for (Int_t i=0;i<n;++i) {
+  if (ntracks<2) return;
+
+
+  // dE/dx,pt and ACORDE study --> studies which need the pair selection    
+  for (Int_t i=0;i<ntracks;++i) {
     AliESDtrack *track1 = event->GetTrack(i);
      
     Double_t d1[3];
     track1->GetDirection(d1);
     
-    for (Int_t j=i+1;j<n;++j) {
+    for (Int_t j=i+1;j<ntracks;++j) {
      AliESDtrack *track2 = event->GetTrack(j);   
      Double_t d2[3];
      track2->GetDirection(d2);
@@ -187,18 +197,6 @@ void AliTPCcalibCosmic::Process(AliESDEvent *event) {
        fModules->Fill(z, x);
       }
       
-      printf("My stream level=%d\n",fStreamLevel);
-
-      if (fStreamLevel>0){
-       TTreeSRedirector * cstream =  GetDebugStreamer();
-       printf("My stream=%p\n",cstream);
-       if (cstream) {
-         (*cstream) << "Track" <<
-           "Track1.=" << track1 <<      // original track 1
-           "Track2.=" << track2 <<      // original track2
-           "\n";
-       }
-      }
   //     AliExternalTrackParam * trackOut = new AliExternalTrackParam(*track2->GetOuterParam());
 //       AliTracker::PropagateTrackTo(trackOut,850.,105.658,30);
 //       delete trackOut;
@@ -217,11 +215,174 @@ void AliTPCcalibCosmic::Process(AliESDEvent *event) {
   
 }    
 
+void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
+  //
+  // Find cosmic pairs
+  //
+  //
+  if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
+  AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
+  Int_t ntracks=event->GetNumberOfTracks(); 
+  fHistNTracks->Fill(ntracks);
+  TObjArray  tpcSeeds(ntracks);
+  if (ntracks==0) return;
+  Double_t vtxx[3]={0,0,0};
+  Double_t svtxx[3]={0.000001,0.000001,100.};
+  AliESDVertex vtx(vtxx,svtxx);
+  //
+  //track loop
+  //
+  for (Int_t i=0;i<ntracks;++i) { 
+   AliESDtrack *track = event->GetTrack(i); 
+   fClusters->Fill(track->GetTPCNcls());   
+   AliExternalTrackParam * trackIn = new AliExternalTrackParam(*track->GetInnerParam());
+   
+   AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
+   TObject *calibObject;
+   AliTPCseed *seed = 0;
+   for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
+     if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
+   }
+   if (seed) tpcSeeds.AddAt(seed,i);
+   if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0)); 
+  }
+  if (ntracks<2) return;
+  //
+  // Find pairs
+  //
+  for (Int_t i=0;i<ntracks;++i) {
+    AliESDtrack *track0 = event->GetTrack(i);     
+    Double_t d1[3];
+    track0->GetDirection(d1);    
+    for (Int_t j=i+1;j<ntracks;++j) {
+     AliESDtrack *track1 = event->GetTrack(j);   
+     Double_t d2[3];
+     track1->GetDirection(d2);
+      printf("My stream level=%d\n",fStreamLevel);
+      AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
+      AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
+      if (! seed0) continue;
+      if (! seed1) continue;
+      Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0);
+      Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0);
+      Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
+      Float_t d0  = track0->GetLinearD(0,0);
+      Float_t d1  = track1->GetLinearD(0,0);
+      //
+      // conservative cuts - convergence to be guarantied
+      // applying before track propagation
+      if (TMath::Abs(d0+d1)>fCutMaxD) continue;   // distance to the 0,0
+      if (dir>fCutMinDir) continue;               // direction vector product
+      Float_t bz = AliTracker::GetBz();
+      Float_t dvertex0[2];   //distance to 0,0
+      Float_t dvertex1[2];   //distance to 0,0 
+      track0->GetDZ(0,0,0,bz,dvertex0);
+      track1->GetDZ(0,0,0,bz,dvertex1);
+      if (TMath::Abs(dvertex0[1])>250) continue;
+      if (TMath::Abs(dvertex1[1])>250) continue;
+      //
+      //
+      //
+      Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
+      AliExternalTrackParam param0(*track0);
+      AliExternalTrackParam param1(*track1);
+      //
+      // Propagate using Magnetic field and correct fo material budget
+      //
+      AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
+      AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
+      //
+      // Propagate rest to the 0,0 DCA - z should be ignored
+      //
+      Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
+      Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
+      //      
+      param0.GetDZ(0,0,0,bz,dvertex0);
+      param1.GetDZ(0,0,0,bz,dvertex1);
+      //
+      Double_t xyz0[3];//,pxyz0[3];
+      Double_t xyz1[3];//,pxyz1[3];
+      param0.GetXYZ(xyz0);
+      param1.GetXYZ(xyz1);
+      Bool_t isPair = IsPair(&param0,&param1);
+      //
+      if (fStreamLevel>0){
+       TTreeSRedirector * cstream =  GetDebugStreamer();
+       printf("My stream=%p\n",(void*)cstream);
+       if (cstream) {
+         (*cstream) << "Track0" <<
+           "dir="<<dir<<               //  direction
+           "OK="<<isPair<<             // will be accepted
+           "b0="<<b0<<                 //  propagate status
+           "b1="<<b1<<                 //  propagate status
+           "Orig0.=" << track0 <<      //  original track  0
+           "Orig1.=" << track1 <<      //  original track  1
+           "Tr0.="<<&param0<<          //  track propagated to the DCA 0,0
+           "Tr1.="<<&param1<<          //  track propagated to the DCA 0,0        
+           "v00="<<dvertex0[0]<<       //  distance using kalman
+           "v01="<<dvertex0[1]<<       // 
+           "v10="<<dvertex1[0]<<       //
+           "v11="<<dvertex1[1]<<       // 
+           "d0="<<d0<<                 //  linear distance to 0,0
+           "d1="<<d1<<                 //  linear distance to 0,0
+           //
+           "x00="<<xyz0[0]<<
+           "x01="<<xyz0[1]<<
+           "x02="<<xyz0[2]<<
+           //
+           "x10="<<xyz1[0]<<
+           "x11="<<xyz1[1]<<
+           "x12="<<xyz1[2]<<
+           //
+           "Seed0.=" << track0 <<      //  original seed 0
+           "Seed1.=" << track1 <<      //  original seed 1
+           "dedx0="<<dedx0<<           //  dedx0
+           "dedx1="<<dedx1<<           //  dedx1
+           "\n";
+       }
+      }      
+    }
+  }  
+}    
+
+
 
-Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {
 
+
+Long64_t AliTPCcalibCosmic::Merge(TCollection */*li*/) {
+  
 }
 
+Bool_t  AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
+  //
+  //
+  /*
+  // 0. Same direction - OPOSITE  - cutDir +cutT    
+  TCut cutDir("cutDir","dir<-0.99")
+  // 1. 
+  TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
+  //
+  // 2. The same rphi 
+  TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
+  //
+  //
+  //
+  TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");  
+  // 1/Pt diff cut
+  */
+  const Double_t *p0 = tr0->GetParameter();
+  const Double_t *p1 = tr1->GetParameter();
+  if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
+  if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
+  Double_t d0[3], d1[3];
+  tr0->GetDirection(d0);    
+  tr1->GetDirection(d1);       
+  if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
+  //
+  return kTRUE;  
+}
+
+
 
 void AliTPCcalibCosmic::BinLogX(TH1 *h) {
 
index 5b7729f..2a2de44 100644 (file)
@@ -21,8 +21,8 @@ public:
   
   virtual void     Process(AliESDEvent *event);
   virtual Long64_t Merge(TCollection *li);
-  
-  
+  void             FindPairs(AliESDEvent *event);
+  Bool_t           IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1);
 public:
   static void BinLogX(TH1 * h);   // method for correct histogram binning
 
@@ -32,7 +32,11 @@ public:
   TH1F  *fHistPt;
   TH1F  *fPtResolution;
   TH2F  *fDeDx;
-
+  // cuts
+  //
+  Float_t fCutMaxD;     // maximal distance in rfi ditection
+  Float_t fCutTheta;    // maximal distance in theta ditection
+  Float_t fCutMinDir;   // direction vector products
   AliTPCcalibCosmic(const AliTPCcalibCosmic&); 
   AliTPCcalibCosmic& operator=(const AliTPCcalibCosmic&);