Small update

[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEpidTRD.cxx

diff --git a/PWGHF/hfe/AliHFEpidTRD.cxx b/PWGHF/hfe/AliHFEpidTRD.cxx
index 13ff97f9480778f5e423b6f3caa7b09803aa79bc..6278ebab78229883faee3541221a0d74c2d1355e 100644 (file)
--- a/PWGHF/hfe/AliHFEpidTRD.cxx
+++ b/PWGHF/hfe/AliHFEpidTRD.cxx
@@ -53,6 +53,7 @@ AliHFEpidTRD::AliHFEpidTRD() :
   , fRunNumber(0)
   , fElectronEfficiency(0.90)
   , fTotalChargeInSlice0(kFALSE)
+  , fTRD2DPID(kFALSE)
 {
   //
   // default  constructor
@@ -70,6 +71,7 @@ AliHFEpidTRD::AliHFEpidTRD(const char* name) :
   , fRunNumber(0)
   , fElectronEfficiency(0.91)
   , fTotalChargeInSlice0(kFALSE)
+  , fTRD2DPID(kFALSE)
 {
   //
   // default  constructor
@@ -87,6 +89,7 @@ AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
   , fRunNumber(ref.fRunNumber)
   , fElectronEfficiency(ref.fElectronEfficiency)
   , fTotalChargeInSlice0(ref.fTotalChargeInSlice0)
+  , fTRD2DPID(ref.fTRD2DPID)
 {
   //
   // Copy constructor
@@ -118,6 +121,7 @@ void AliHFEpidTRD::Copy(TObject &ref) const {
   target.fCutNTracklets = fCutNTracklets;
   target.fRunNumber = fRunNumber;
   target.fTotalChargeInSlice0 = fTotalChargeInSlice0;
+  target.fTRD2DPID = fTRD2DPID;
   target.fElectronEfficiency = fElectronEfficiency;
   memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
   AliHFEpidBase::Copy(ref);
@@ -132,6 +136,21 @@ AliHFEpidTRD::~AliHFEpidTRD(){
 
 //______________________________________________________
 Bool_t AliHFEpidTRD::InitializePID(Int_t run){
+  //
+  // InitializePID call different init function depending on TRD PID method
+  //
+  //
+
+  //if(fTRD2DPID) return Initialize2D(run);
+  if(fTRD2DPID) return kTRUE;
+  else return Initialize1D(run);
+
+
+
+}
+
+//______________________________________________________
+Bool_t AliHFEpidTRD::Initialize1D(Int_t run){
   //
   // InitializePID: Load TRD thresholds and create the electron efficiency axis
   // to navigate 
@@ -147,8 +166,35 @@ Bool_t AliHFEpidTRD::InitializePID(Int_t run){
   return kFALSE;
 }
 
+/*
+//______________________________________________________
+Bool_t AliHFEpidTRD::Initialize2D(Int_t run){
+  //
+  // Initialize2DimPID
+  //
+  //
+
+    return kTRUE;
+
+}
+*/
+
 //______________________________________________________
 Int_t AliHFEpidTRD::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const {
+  //
+  // Does PID for TRD alone:
+  // PID algorithm selected according to flag
+  //
+  //
+
+    if(fTRD2DPID) return IsSelected2D(track, pidqa);
+    else return IsSelected1D(track, pidqa);
+
+
+}
+
+//______________________________________________________
+Int_t AliHFEpidTRD::IsSelected1D(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const {
   //
   // Does PID for TRD alone:
   // PID thresholds based on 90% Electron Efficiency level approximated by a linear 
@@ -207,6 +253,56 @@ Int_t AliHFEpidTRD::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager
   }
   return 211;
 
+}
+
+//______________________________________________________
+Int_t AliHFEpidTRD::IsSelected2D(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const {
+  //
+  // 2D TRD PID
+  // 
+  // 
+  //
+  AliDebug(2, "Applying TRD PID");
+  if(!fkPIDResponse){
+    AliDebug(2, "Cannot process track");
+    return 0;
+  }
+
+  AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis: AliHFEpidObject::kAODanalysis;
+  Double_t p = GetP(track->GetRecTrack(), anatype);
+  if(p < fMinP){ 
+    AliDebug(2, Form("Track momentum below %f", fMinP));
+    return 0;
+  }
+
+  if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kBeforePID); 
+
+  if(fCutNTracklets > 0){
+    AliDebug(1, Form("Number of tracklets cut applied: %d\n", fCutNTracklets));
+    Int_t ntracklets = track->GetRecTrack() ? track->GetRecTrack()->GetTRDntrackletsPID() : 0;
+    if(TestBit(kExactTrackletCut)){
+      AliDebug(1, Form("Exact cut applied: %d tracklets found\n", ntracklets));
+      if(ntracklets != fCutNTracklets) return 0;
+    } else {
+      AliDebug(1, Form("Greater Equal cut applied: %d tracklets found\n", ntracklets));
+      if(ntracklets < fCutNTracklets) return 0;
+    }
+  }
+  AliDebug(1,"Track selected\n");
+
+  Int_t centralitybin = track->IsPbPb() ? track->GetCentrality() : -2;
+  Float_t fCentralityLimitsdefault[12]= {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
+  Float_t centrality=-1;
+  if(centralitybin>=0) centrality=fCentralityLimitsdefault[centralitybin]+1;
+
+  if(fkPIDResponse->IdentifiedAsElectronTRD(track->GetRecTrack(),fElectronEfficiency,centrality,AliTRDPIDResponse::kLQ2D)){
+      AliDebug(2, Form("Electron effi: %f %i %i %f\n", fElectronEfficiency,track->GetCentrality(),centralitybin,centrality));
+      if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kAfterPID); 
+      return 11;
+  } else return 211;
+
+
+
 }
 
 //___________________________________________________________________
@@ -277,7 +373,8 @@ Double_t AliHFEpidTRD::GetElectronLikelihood(const AliVTrack *track, AliHFEpidOb
     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
     if(esdtrack) esdtrack->GetTRDpid(pidProbs);
   } else {
-    fkPIDResponse->ComputeTRDProbability(track, AliPID::kSPECIES, pidProbs);
+      if(fTRD2DPID) fkPIDResponse->ComputeTRDProbability(track, AliPID::kSPECIES, pidProbs,AliTRDPIDResponse::kLQ2D); 
+      else fkPIDResponse->ComputeTRDProbability(track, AliPID::kSPECIES, pidProbs,AliTRDPIDResponse::kLQ1D);
   }
   if(!IsRenormalizeElPi()) return pidProbs[AliPID::kElectron];
   Double_t probsNew[AliPID::kSPECIES];
@@ -323,9 +420,9 @@ Double_t AliHFEpidTRD::GetChargeLayer(const AliVParticle *track, UInt_t layer, A
     AliAODPid *aoddetpid = aodtrack ? aodtrack->GetDetPid() : NULL;
     if(aoddetpid){
       if(fTotalChargeInSlice0)
-        charge = aoddetpid->GetTRDsignal()[layer * aoddetpid->GetTRDnSlices()];
+        charge = aoddetpid->GetTRDslices()[layer * aoddetpid->GetTRDnSlices()];
       else
-       for(Int_t islice = 0; islice < aoddetpid->GetTRDnSlices(); islice++) charge += aoddetpid->GetTRDsignal()[layer * aoddetpid->GetTRDnSlices() + islice];
+       for(Int_t islice = 0; islice < aoddetpid->GetTRDnSlices(); islice++) charge += aoddetpid->GetTRDslices()[layer * aoddetpid->GetTRDnSlices() + islice];
     }
   }
   return charge;
@@ -351,20 +448,22 @@ Double_t AliHFEpidTRD::GetTRDSignalV1(const AliESDtrack *track, Float_t truncati
   // Calculate mean over the last 2/3 slices
   //
   const Int_t kNSlices = 48;
-  const Int_t kSlicePerLayer = 7;
+  const Int_t kLastSlice = 6; // Slice 7 is taken out from the truncated mean calculation
   const Double_t kVerySmall = 1e-12;
   // Weight the slice to equalize the MPV of the dQ/dl-distribution per slice to the one in the first slice
   // Pions are used as reference for the equalization
   const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
+  const Double_t kWeightSliceNo0[8] = {1., 1., 1.271, 1.451, 1.531, 1.543, 1.553, 2.163};  // Weighting factors in case slice 0 stores the total charge
+  const Double_t *kWeightFactor = fTotalChargeInSlice0 ? kWeightSliceNo0 : kWeightSlice;
   AliDebug(3, Form("Number of Tracklets: %d\n", track->GetTRDntrackletsPID()));
   Double_t trdSlices[kNSlices], tmp[kNSlices];
   Int_t indices[48];
   Int_t icnt = 0;
   for(Int_t idet = 0; idet < 6; idet++)
-    for(Int_t islice = fTotalChargeInSlice0 ? 1 : 0 ; islice < kSlicePerLayer; islice++){
+    for(Int_t islice = fTotalChargeInSlice0 ? 1 : 0 ; islice <= kLastSlice; islice++){
       AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice)));
       if(TMath::Abs(track->GetTRDslice(idet, islice)) < kVerySmall) continue;;
-      trdSlices[icnt++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
+      trdSlices[icnt++] = track->GetTRDslice(idet, islice) * kWeightFactor[islice];
     }
   AliDebug(1, Form("Number of Slices: %d\n", icnt));
   if(icnt < 6) return 0.;   // We need at least 6 Slices for the truncated mean