ATO-525,ATO-496 - addig qptTgl correction to slove dEdx splitting at pt<0.5 GeV

STAT/AliTreePlayer.cxx

@@ -1172,7 +1172,7 @@ TObjArray  * AliTreePlayer::MakeHistograms(TTree * tree, TString hisString, TStr
   Int_t nExpressions=hisString.CountChar(':')+hisString.CountChar(';')+1;
   TObjArray * formulaArray   = new TObjArray(nExpressions);    // array of all expressions  - OWNER
   TString queryString = "";
-  Int_t hisSizeFull=0;
+  Long_t hisSizeFull=0;
   //
   //  1.) Analyze formula, book list of TObjString
   //
@@ -1307,7 +1307,7 @@ TObjArray  * AliTreePlayer::MakeHistograms(TTree * tree, TString hisString, TStr
         }
       }
       if (verbose&0x1) {
-        ::Info("AliTreePlayer::MakeHistograms","Total size=%d",hisSizeFull);
+        ::Info("AliTreePlayer::MakeHistograms","Total size=%lld",hisSizeFull);
       }
     }
     //    2.3 fill histograms
@@ -1647,6 +1647,7 @@ void AliTreePlayer::MakeCacheTreeChunk(TTree * tree, TString varList, TString ou
   for (Int_t iEntry=firstEntry; iEntry<nEntries; iEntry+=chunkSize) {
     ::Info("Processing chunk","%d",iEntry);
     if (estimate < chunkSize) tree->SetEstimate(chunkSize);
+    if (chunkSize+iEntry>=entriesAll) chunkSize=entriesAll-iEntry;   // ROOT6 does not handle properly query above limit
     Int_t entries = tree->Draw(varList.Data(), selection, "goffpara", chunkSize, iEntry);
     if (entries<=0) break;
     if (entries > estimate) {

STAT/TStatToolkit.cxx

@@ -1221,17 +1221,31 @@ TGraph * TStatToolkit::MakeGraphSparse(TTree * tree, const char * expr, const ch
     graphNew->GetXaxis()->SetBinLabel(i+1,xName);
     graphNew->GetX()[i]+=offset;
   }
-  if (tree->GetVar(1)->IsInteger() && strlen(tree->GetHistogram()->GetXaxis()->GetBinLabel(1))>0){    
-    for(Int_t i=0;i<count;i++){
-      graphNew->GetXaxis()->SetBinLabel(i+1,tree->GetHistogram()->GetXaxis()->GetBinLabel(i+1));
+  if (tree->GetV3()==NULL) {
+    if (tree->GetVar(1)->IsInteger() && strlen(tree->GetHistogram()->GetXaxis()->GetBinLabel(1)) > 0) {
+      for (Int_t i = 0; i < count; i++) {
+        graphNew->GetXaxis()->SetBinLabel(i + 1, tree->GetHistogram()->GetXaxis()->GetBinLabel(i + 1));
+      }
     }
-  }
-  if (tree->GetVar(0)->IsInteger() &&  strlen(tree->GetHistogram()->GetXaxis()->GetBinLabel(1))>0 ){
-    for(Int_t i=0;i<count;i++){
-      graphNew->GetYaxis()->SetBinLabel(i+1,tree->GetHistogram()->GetYaxis()->GetBinLabel(i+1));
+    if (tree->GetVar(0)->IsInteger() && strlen(tree->GetHistogram()->GetYaxis()->GetBinLabel(1)) > 0) {
+      for (Int_t i = 0; i < count; i++) {
+        graphNew->GetYaxis()->SetBinLabel(i + 1, tree->GetHistogram()->GetYaxis()->GetBinLabel(i + 1));
+      }
+    }
+  }else{
+    if (tree->GetVar(1)->IsInteger() && strlen(tree->GetHistogram()->GetYaxis()->GetBinLabel(1)) > 0) {
+      for (Int_t i = 0; i < count; i++) {
+        graphNew->GetXaxis()->SetBinLabel(i + 1, tree->GetHistogram()->GetYaxis()->GetBinLabel(i + 1));
+      }
+    }
+    if (tree->GetVar(0)->IsInteger() && strlen(tree->GetHistogram()->GetZaxis()->GetBinLabel(1)) > 0) {
+      for (Int_t i = 0; i < count; i++) {
+        graphNew->GetYaxis()->SetBinLabel(i + 1, tree->GetHistogram()->GetZaxis()->GetBinLabel(i + 1));
+      }
     }
   }
 
+
   graphNew->GetHistogram()->SetTitle("");
   graphNew->SetMarkerStyle(mstyle);
   graphNew->SetMarkerColor(mcolor);  graphNew->SetLineColor(mcolor);

STEER/ESD/AliESDv0.cxx

@@ -986,15 +986,23 @@ Double_t AliESDv0::GetKFInfoScale(UInt_t p1, UInt_t p2, Int_t type, Double_t d1p
   // calculate delta of energy loss
   Double_t bg1=paramP.P()  /TDatabasePDG::Instance()->GetParticle(spdg[p1])->Mass();
   Double_t bg2=paramN.P()  /TDatabasePDG::Instance()->GetParticle(spdg[p2])->Mass();
-  Double_t dP1=AliExternalTrackParam::BetheBlochGeant(bg1)/AliExternalTrackParam::BetheBlochGeant(3.); ///relative energu loss
-  Double_t dP2=AliExternalTrackParam::BetheBlochGeant(bg2)/AliExternalTrackParam::BetheBlochGeant(3.);
-
+  Double_t dE1=AliExternalTrackParam::BetheBlochGeant(bg1)/AliExternalTrackParam::BetheBlochGeant(3.); ///relative energu loss
+  Double_t dE2=AliExternalTrackParam::BetheBlochGeant(bg2)/AliExternalTrackParam::BetheBlochGeant(3.);
+  dE1*=TMath::Sqrt(1+(paramP.GetTgl()*paramP.GetTgl()));
+  dE2*=TMath::Sqrt(1+(paramN.GetTgl()*paramN.GetTgl()));
+  Double_t mass1=AliPID::ParticleMass(p1);
+  Double_t mass2=AliPID::ParticleMass(p2);
+  Double_t E1 = TMath::Sqrt(paramP.P()*paramP.P()+mass1*mass1)+eLoss*dE1;
+  Double_t E2 = TMath::Sqrt(paramN.P()*paramN.P()+mass2*mass2)+eLoss*dE2;
+  Double_t dP1=TMath::Sqrt(TMath::Max(E1*E1-mass1*mass1,0.0))-paramP.P();
+  Double_t dP2=TMath::Sqrt(TMath::Max(E2*E2-mass2*mass2,0.0))-paramN.P();
+  //printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n",dE1,dE2, E1,E2, dP1,dP2,eLoss);
   Double_t *pparam1 = (Double_t*)paramP.GetParameter();
   Double_t *pparam2 = (Double_t*)paramN.GetParameter();
   if (flag&0x1) pparam1[4]+=d1pt;
   if (flag&0x2) pparam2[4]+=d1pt;
-  if (flag&0x1) pparam1[4]*=(1+s1pt)*(1.+eLoss*dP1/paramP.P());   /// TODO   - use relative energy loss
-  if (flag&0x2) pparam2[4]*=(1+s1pt)*(1.+eLoss*dP2/paramN.P());   ///
+  if (flag&0x1) pparam1[4]*=(1+s1pt)*(1.+dP1/paramP.P());   /// TODO   - use relative energy loss
+  if (flag&0x2) pparam2[4]*=(1+s1pt)*(1.+dP2/paramN.P());   ///
 
   //
   AliKFParticle kfp1( paramP, spdg[p1] *TMath::Sign(1,p1) );

STEER/STEERBase/AliTPCPIDResponse.cxx

@@ -78,6 +78,7 @@ AliTPCPIDResponse::AliTPCPIDResponse():
   fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios),
   fOADBContainer(0x0),
   fPileupCorrection(0x0),
+  fQPtTglCorrection(0x0),
   fVoltageMap(72),
   fLowGainIROCthreshold(-40),
   fBadIROCthreshhold(-70),
@@ -105,6 +106,7 @@ AliTPCPIDResponse::AliTPCPIDResponse():
   fOROClongWeight(1.),
   fPileupCorrectionStrategy(kPileupCorrectionInExpectedSignal),
   fPileupCorrectionRequested(kFALSE),
+  fQPtTglCorrectionRequested(kFALSE),
   fSigmaParametrization(0x0),
   fMultiplicityNormalization(1),
   fRecoPassNameUsed(),
@@ -243,6 +245,7 @@ AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& that):
   fOROClongWeight(that.fOROClongWeight),
   fPileupCorrectionStrategy(that.fPileupCorrectionStrategy),
   fPileupCorrectionRequested(that.fPileupCorrectionRequested),
+  fQPtTglCorrectionRequested(that.fQPtTglCorrectionRequested),
   fSigmaParametrization(that.fSigmaParametrization),
   fMultiplicityNormalization(that.fMultiplicityNormalization),
   fRecoPassNameUsed(that.fRecoPassNameUsed),
@@ -367,7 +370,7 @@ AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that)
   fOROClongWeight            = that.fOROClongWeight;
   fPileupCorrectionStrategy  = that.fPileupCorrectionStrategy;
   fPileupCorrectionRequested = that.fPileupCorrectionRequested;
-  fPileupCorrectionRequested = that.fPileupCorrectionRequested;
+  fQPtTglCorrectionRequested = that.fQPtTglCorrectionRequested;
   fSigmaParametrization      = that.fSigmaParametrization;
   fMultiplicityNormalization = that.fMultiplicityNormalization;
   fRecoPassNameUsed          = that.fRecoPassNameUsed;
@@ -493,7 +496,7 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
                                               const TSpline3* responseFunction,
                                               Bool_t correctEta,
                                               Bool_t correctMultiplicity,
-                                              Bool_t usePileupCorrection) const 
+                                              Bool_t usePileupCorrection, Bool_t useQPtTglCorrection) const
 {
   // Calculates the expected PID signal as the function of 
   // the information stored in the track and the given parameters,
@@ -523,6 +526,11 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
   if (usePileupCorrection && fPileupCorrectionStrategy == kPileupCorrectionInExpectedSignal) {
     corrPileup = GetPileupCorrectionValue(track);
   }
+  //qPtTglCorrection
+   Double_t corrqPtTgl = 1.;
+  if (useQPtTglCorrection && fQPtTglCorrectionRequested) {
+    corrqPtTgl = GetQPtTglCorrectionValue(track);
+  }
 
   if (!correctEta && !correctMultiplicity) {
     return dEdxSplines + corrPileup;
@@ -541,7 +549,7 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
     corrFactorMultiplicity = GetMultiplicityCorrectionFast(track, dEdxSplines * corrFactorEta, fCurrentEventMultiplicity);
   }
 
-  return dEdxSplines * corrFactorEta * corrFactorMultiplicity + corrPileup;
+  return dEdxSplines * corrFactorEta * corrFactorMultiplicity * corrqPtTgl + corrPileup;
 }
 
 
@@ -551,7 +559,7 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
                                               ETPCdEdxSource dedxSource,
                                               Bool_t correctEta,
                                               Bool_t correctMultiplicity,
-                                              Bool_t usePileupCorrection) const
+                                              Bool_t usePileupCorrection, Bool_t useQPtTglCorrection) const
 {
   // Calculates the expected PID signal as the function of 
   // the information stored in the track, for the specified particle type 
@@ -582,7 +590,7 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track,
   }
 
   // Charge factor already taken into account inside the following function call
-  return GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity, usePileupCorrection);
+  return GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity, usePileupCorrection, useQPtTglCorrection);
 }
 
 //_________________________________________________________________________
@@ -646,7 +654,7 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
                                              const TSpline3* responseFunction,
                                              Bool_t correctEta,
                                              Bool_t correctMultiplicity,
-                                             Bool_t usePileupCorrection) const 
+                                             Bool_t usePileupCorrection, Bool_t useQPtTglCorrection) const
 {
   // Calculates the expected sigma of the PID signal as the function of 
   // the information stored in the track and the given parameters,
@@ -671,10 +679,10 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
   // If no sigma map is available or if no eta correction is requested (sigma maps only for corrected eta!), use the old parametrisation
   if (!fhEtaSigmaPar1 || !correctEta) {  
     if (nPoints != 0) 
-      return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity, usePileupCorrection) *
+      return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity, usePileupCorrection, useQPtTglCorrection) *
                fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints);
     else
-      return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity, usePileupCorrection)*fRes0[gainScenario];
+      return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity, usePileupCorrection, useQPtTglCorrection)*fRes0[gainScenario];
   }
 
   if (nPoints > 0) {
@@ -711,7 +719,7 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
                                              ETPCdEdxSource dedxSource,
                                              Bool_t correctEta,
                                              Bool_t correctMultiplicity,
-                                             Bool_t usePileupCorrection) const 
+                                             Bool_t usePileupCorrection, Bool_t useQPtTglCorrection) const
 {
   // Calculates the expected sigma of the PID signal as the function of 
   // the information stored in the track, for the specified particle type 
@@ -726,7 +734,7 @@ Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track,
   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
     return 999; //TODO: better handling!
 
-  return GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity, usePileupCorrection);
+  return GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity, usePileupCorrection, useQPtTglCorrection);
 }
 
 
@@ -736,7 +744,7 @@ Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track,
                              ETPCdEdxSource dedxSource,
                              Bool_t correctEta,
                              Bool_t correctMultiplicity,
-                             Bool_t usePileupCorrection) const
+                             Bool_t usePileupCorrection, Bool_t useQPtTglCorrection) const
 {
   //Calculates the number of sigmas of the PID signal from the expected value
   //for a given particle species in the presence of multiple gain scenarios
@@ -750,8 +758,8 @@ Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track,
   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
     return -999; //TODO: Better handling!
 
-  Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity, usePileupCorrection);
-  Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity, usePileupCorrection);
+  Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity, usePileupCorrection, useQPtTglCorrection);
+  Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity, usePileupCorrection, useQPtTglCorrection);
   // 999 will be returned by GetExpectedSigma e.g. in case of 0 dEdx clusters
   if (sigma >= 998) 
     return -999;
@@ -765,7 +773,7 @@ Float_t AliTPCPIDResponse::GetSignalDelta(const AliVTrack* track,
                                           ETPCdEdxSource dedxSource,
                                           Bool_t correctEta,
                                           Bool_t correctMultiplicity,
-                                          Bool_t usePileupCorrection /*= kFALSE*/,
+                                          Bool_t usePileupCorrection /*= kFALSE*/, Bool_t useQPtTglCorrection,
                                           Bool_t ratio/*=kFALSE*/)const
 {
   //Calculates the number of sigmas of the PID signal from the expected value
@@ -780,7 +788,7 @@ Float_t AliTPCPIDResponse::GetSignalDelta(const AliVTrack* track,
   if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction))
     return -9999.; //TODO: Better handling!
 
-  const Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity, usePileupCorrection);
+  const Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity, usePileupCorrection, useQPtTglCorrection);
 
   Double_t delta=-9999.;
   if (!ratio) delta=dEdx-bethe;
@@ -1837,6 +1845,28 @@ Double_t AliTPCPIDResponse::GetPileupCorrectionValue(const AliVTrack* track) con
   return corrPileup;
 }
 
+
+//_____________________________________________________________________________
+Double_t AliTPCPIDResponse::GetQPtTglCorrectionValue(const AliVTrack* track) const
+{
+  //
+  // The  correction is an additive value. The corrected dEdx is
+  // dEdx - corrQPtTgl
+  //
+  if (!fQPtTglCorrection) {
+    return 0.;
+  }
+
+  //const Double_t trackTgl = TMath::Abs(TMath::SinH(track->GetTPCTgl()));
+  const Double_t trackTgl = track->GetTgl();
+  Double_t corrVals[2] = {track->GetInnerParam()->GetSigned1Pt(), trackTgl};
+  const Double_t corrQPtTgl = (1+fQPtTglCorrection->Eval(corrVals)) ;
+
+  return corrQPtTgl;
+}
+
+
+
 //______________________________________________________________________________
 AliNDLocalRegression* AliTPCPIDResponse::GetPileupCorrectionFromFile(const TString fileName)
 {
@@ -1874,6 +1904,45 @@ AliNDLocalRegression* AliTPCPIDResponse::GetPileupCorrectionFromFile(const TStri
   return ndLocal;
 }
 
+
+//______________________________________________________________________________
+AliNDLocalRegression* AliTPCPIDResponse::GetQPtTglCorrectionFromFile(const TString fileName)
+{
+  if (fileName.Contains("alien://") && !gGrid) {
+    TGrid::Connect("alien");
+  }
+
+  TFile* f = TFile::Open(fileName);
+  if (!f || !f->IsOpen() || f->IsZombie()) {
+    printf("AliTPCPIDResponse::GetQPtTglCorrectionFromFile: Could not open QPtTgl correction file: %s", fileName.Data());
+    delete f;
+    return 0x0;
+  }
+
+  // Assume that the only object in the file is the QPtTgl correction
+  const Int_t numberOfKeys = f->GetListOfKeys()->GetEntries();
+  if ( numberOfKeys != 1) {
+    printf("AliTPCPIDResponse::GetQPtTglCorrectionFromFile: Number of objects in the file %d is not 1", numberOfKeys);
+    delete f;
+    return 0x0;
+  }
+
+  const char* objectName = f->GetListOfKeys()->At(0)->GetName();
+  TObject* o = f->Get(objectName);
+  AliNDLocalRegression* ndLocal = dynamic_cast<AliNDLocalRegression*>(o);
+  if (!ndLocal) {
+    printf("AliTPCPIDResponse::GetQPtTglCorrectionFromFile: Object in file %s with name %s is not of type AliNDLocalRegression but %s", fileName.Data(), objectName, o->IsA()->GetName());
+    delete f;
+    return 0x0;
+  }
+
+  ndLocal->SetName("QPtTglCorrection");
+
+  delete f;
+  return ndLocal;
+}
+
+
 //_____________________________________________________________________________
 Bool_t AliTPCPIDResponse::InitFromOADB(const Int_t run, const Int_t pass, TString passName,
                                        const char* oadbFile/*="$ALICE_PHYSICS/OADB/COMMON/PID/data/TPCPIDResponseOADB.root"*/,