Go back to previous PID method

diff --git a/TRD/AliTRDtracker.cxx b/TRD/AliTRDtracker.cxx
index 9efdabc4de9ea76068c72b4a84c9c459c6f1c38b..a180ce49322cd2e24c41175b0f747519d8c76d8c 100644 (file)
--- a/TRD/AliTRDtracker.cxx
+++ b/TRD/AliTRDtracker.cxx
@@ -58,9 +58,9 @@
 
 ClassImp(AliTRDtracker)
 
-const  Float_t  AliTRDtracker::fgkMinClustersInTrack =  0.5;  
-const  Float_t  AliTRDtracker::fgkLabelFraction      =  0.8;  // ??
-const  Double_t AliTRDtracker::fgkMaxChi2            = 12.0; 
+const  Float_t  AliTRDtracker::fgkMinClustersInTrack =  0.5;  //
+const  Float_t  AliTRDtracker::fgkLabelFraction      =  0.8;  //
+const  Double_t AliTRDtracker::fgkMaxChi2            = 12.0;  //
 const  Double_t AliTRDtracker::fgkMaxSnp             =  0.95; // Corresponds to tan = 3
 const  Double_t AliTRDtracker::fgkMaxStep            =  2.0;  // Maximal step size in propagation 
 
@@ -3468,65 +3468,57 @@ Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster *c)
 void AliTRDtracker::CookdEdxTimBin(AliTRDtrack &TRDtrack)
 {
   //
-  // Build the information needed for TRD PID calculations. This are:
-  // - dE/dx - total
-  // - dE/dx - slices (0-14, 15-22, NONE). The Edep per slice is
-  //    calculated as the sum of clusters charge weighted with average
-  //    value from the <PH> spectrum
-  // - Time bin for Max. Cluster
-  //
-  // Authors:
+  // This is setting fdEdxPlane and fTimBinPlane
+  // Sums up the charge in each plane for track TRDtrack and also get the 
+  // Time bin for Max. Cluster
   // Prashant Shukla (shukla@physi.uni-heidelberg.de)
-  // Alex Bercuci (a.bercuci@gsi.de)
+  //
 
   Double_t  clscharge[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
   Double_t  maxclscharge[AliESDtrack::kNPlane];
   Int_t     nCluster[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
   Int_t     timebin[AliESDtrack::kNPlane];
 
-  // Initialization of variables  
+  // Initialization of cluster charge per plane.  
   for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
     for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
-      clscharge[iPlane][iSlice] = 0.;
+      clscharge[iPlane][iSlice] = 0.0;
       nCluster[iPlane][iSlice]  = 0;
     }
+  }
+
+  // Initialization of cluster charge per plane.  
+  for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
     timebin[iPlane]      =  -1;
     maxclscharge[iPlane] = 0.0;
   }
 
-       // Average PH spectrum used for weighting the edep. It has to
-       // come eather from caligration DB or as parameter of the class. To
-       // be discuss with software responsibles
-       const Float_t fAveragePH[] = { 5.010e-03, 2.017e-02, 6.221e-02, 6.706e-02
-                                     , 5.551e-02, 4.812e-02, 4.395e-02, 4.219e-02
-                                     , 4.172e-02, 4.156e-02, 4.162e-02, 4.204e-02
-                                     , 4.278e-02, 4.367e-02, 4.460e-02, 4.554e-02
-                                     , 4.674e-02, 4.810e-02, 5.005e-02, 5.258e-02
-                                     , 5.587e-02, 5.892e-02, 5.587e-02, 5.258e-02 };
-       
   // Loop through all clusters associated to track TRDtrack
-  AliTRDcluster *pTRDcluster = 0x0;
-  for (Int_t iClus = 0; iClus < TRDtrack.GetNumberOfClusters(); iClus++) {
+  Int_t nClus = TRDtrack.GetNumberOfClusters();  // from Kalmantrack
+  for (Int_t iClus = 0; iClus < nClus; iClus++) {
     Double_t charge = TRDtrack.GetClusterdQdl(iClus);
     Int_t    index  = TRDtrack.GetClusterIndex(iClus);
-    if (!(pTRDcluster = (AliTRDcluster*)GetCluster(index))) continue;
-    
-               // check negative time bins ?!?!?
-               Int_t    tb     = pTRDcluster->GetLocalTimeBin();
-
-               // temporary until fix in GetT0
-               if(tb<0){
-                       Warning("CookdEdxTimBin()", Form("Negative time bin value in track %d cluster %d", TRDtrack.GetLabel(), iClus));
-                       continue;
-               }
-               
+    AliTRDcluster *pTRDcluster = (AliTRDcluster *) GetCluster(index); 
+    if (!pTRDcluster) {
+      continue;
+    }
+    Int_t    tb     = pTRDcluster->GetLocalTimeBin();
+    if (!tb) {
+      continue;
+    }
     Int_t detector  = pTRDcluster->GetDetector();
     Int_t iPlane    = fGeom->GetPlane(detector);
-
-               // 2 slices method
-               Int_t iSlice = (tb < AliTRDtrack::kMidTimeBin) ? 0 : 1;
-               // weighted edep with <PH>
-               clscharge[iPlane][iSlice] += charge * fAveragePH[tb];
+    if (iPlane >= AliESDtrack::kNPlane) {
+      AliError(Form("Wrong plane %d",iPlane));
+      continue;
+    }
+    Int_t iSlice    = tb * AliESDtrack::kNSlice 
+                         / AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
+    if (iSlice >= AliESDtrack::kNSlice) {
+      AliError(Form("Wrong slice %d",iSlice));
+      continue;
+    }
+    clscharge[iPlane][iSlice] = clscharge[iPlane][iSlice] + charge;
     if (charge > maxclscharge[iPlane]) {
       maxclscharge[iPlane] = charge;
       timebin[iPlane]      = tb;
@@ -3534,18 +3526,31 @@ void AliTRDtracker::CookdEdxTimBin(AliTRDtrack &TRDtrack)
     nCluster[iPlane][iSlice]++;
   } // End of loop over cluster
 
-  // Setting the fdEdxPlane and fTimBinPlane variables
-  Double_t totalCharge = 0.;
+  // Setting the fdEdxPlane and fTimBinPlane variabales 
+  Double_t totalCharge = 0.0;
+
   for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
-               for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
-      //if (nCluster[iPlane][iSlice]) {
-      //  clscharge[iPlane][iSlice] /= nCluster[iPlane][iSlice];
-      //}
+    for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
+      if (nCluster[iPlane][iSlice]) {
+        clscharge[iPlane][iSlice] /= nCluster[iPlane][iSlice];
+      }
       TRDtrack.SetPIDsignals(clscharge[iPlane][iSlice],iPlane,iSlice);
-      totalCharge += clscharge[iPlane][iSlice];
+      totalCharge = totalCharge+clscharge[iPlane][iSlice];
     }
-    TRDtrack.SetPIDTimBin(timebin[iPlane],iPlane);
+    TRDtrack.SetPIDTimBin(timebin[iPlane],iPlane);     
   }
+
+  // Still needed ????
+  //  Int_t i;
+  //  Int_t nc=TRDtrack.GetNumberOfClusters(); 
+  //  Float_t dedx=0;
+  //  for (i=0; i<nc; i++) dedx += TRDtrack.GetClusterdQdl(i);
+  //  dedx /= nc;
+  //  for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
+  //    TRDtrack.SetPIDsignals(dedx, iPlane);
+  //    TRDtrack.SetPIDTimBin(timbin[iPlane], iPlane);
+  //  }
+
 }
 
 //_____________________________________________________________________________

diff --git a/TRD/Cal/AliTRDCalPIDLQ.cxx b/TRD/Cal/AliTRDCalPIDLQ.cxx
index c8ab5ca56dd9ab14cb34826e254cdcb47c0ed452..9d4435ee507165f2d9169e95414379c105058953 100644 (file)
--- a/TRD/Cal/AliTRDCalPIDLQ.cxx
+++ b/TRD/Cal/AliTRDCalPIDLQ.cxx
@@ -61,8 +61,8 @@ Char_t* AliTRDCalPIDLQ::fpartSymb[AliPID::kSPECIES] = {"EL", "MU", "PI", "KA", "
 AliTRDCalPIDLQ::AliTRDCalPIDLQ()
   :TNamed("pid", "PID for TRD")
   ,fNMom(0)
-  ,fTrackMomentum(0x0)
   ,fNLength(0)
+  ,fTrackMomentum(0x0)
   ,fTrackSegLength(0x0)
   ,fNTimeBins(0)
   ,fMeanChargeRatio(0)
@@ -83,8 +83,8 @@ AliTRDCalPIDLQ::AliTRDCalPIDLQ()
 AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title) 
   :TNamed(name,title)
   ,fNMom(0)
-  ,fTrackMomentum(0x0)
   ,fNLength(0)
+  ,fTrackMomentum(0x0)
   ,fTrackSegLength(0x0)
   ,fNTimeBins(0)
   ,fMeanChargeRatio(0)
@@ -106,8 +106,8 @@ AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title)
 AliTRDCalPIDLQ::AliTRDCalPIDLQ(const AliTRDCalPIDLQ &c) 
   :TNamed(c)
   ,fNMom(c.fNMom)
-  ,fTrackMomentum(0x0)
   ,fNLength(c.fNLength)
+  ,fTrackMomentum(0x0)
   ,fTrackSegLength(0x0)
   ,fNTimeBins(c.fNTimeBins)
   ,fMeanChargeRatio(c.fMeanChargeRatio)
@@ -420,13 +420,15 @@ Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Double_t mom, Double_t *dedx
        }
        ax = hist->GetXaxis(); nbinsx = ax->GetNbins();
        ay = hist->GetYaxis(); nbinsy = ay->GetNbins();
-       Bool_t kX = (dedx[0] < ax->GetBinUpEdge(nbinsx));
-       Bool_t kY = (dedx[1] < ay->GetBinUpEdge(nbinsy));
+       Float_t x = dedx[0]+dedx[1], y = dedx[2];
+  Bool_t kX = (x < ax->GetBinUpEdge(nbinsx));
+       Bool_t kY = (y < ay->GetBinUpEdge(nbinsy));
        if(kX)
-               if(kY) LQ1 = hist->GetBinContent( hist->FindBin(dedx[0], dedx[1])); //fEstimator->Estimate2D2(hist, (Float_t&)dedx[0], (Float_t&)dedx[1]);
-               else LQ1 = hist->GetBinContent(ax->FindBin(dedx[0]), nbinsy);
+               if(kY) LQ1 = hist->GetBinContent( hist->FindBin(x, y)); 
+    //fEstimator->Estimate2D2(hist, x, y);
+               else LQ1 = hist->GetBinContent(ax->FindBin(x), nbinsy);
        else
-               if(kY) LQ1 = hist->GetBinContent(nbinsx, ay->FindBin(dedx[1]));
+               if(kY) LQ1 = hist->GetBinContent(nbinsx, ay->FindBin(y));
                else LQ1 = hist->GetBinContent(nbinsx, nbinsy);
 
 
@@ -436,10 +438,11 @@ Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Double_t mom, Double_t *dedx
                return LQ1;
        }
        if(kX)
-               if(kY) LQ2 = hist->GetBinContent( hist->FindBin(dedx[0], dedx[1])); //fEstimator->Estimate2D2(hist, (Float_t&)dedx[0], (Float_t&)dedx[1]);
-               else LQ2 = hist->GetBinContent(ax->FindBin(dedx[0]), nbinsy);
+               if(kY) LQ2 = hist->GetBinContent( hist->FindBin(x, y)); 
+    //fEstimator->Estimate2D2(hist, x, y);
+               else LQ2 = hist->GetBinContent(ax->FindBin(x), nbinsy);
        else
-               if(kY) LQ2 = hist->GetBinContent(nbinsx, ay->FindBin(dedx[1]));
+               if(kY) LQ2 = hist->GetBinContent(nbinsx, ay->FindBin(y));
                else LQ2 = hist->GetBinContent(nbinsx, nbinsy);