[u/mrichter/AliRoot.git] / ANALYSIS / TenderSupplies / AliTPCTenderSupply.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


///////////////////////////////////////////////////////////////////////////////
//                                                                           //
// TPC tender: reapply pid on the fly                                        //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

#include <TList.h>
#include <TObjArray.h>
#include <TObjString.h>
#include <TChain.h>
#include <TFile.h>
#include <TString.h>
#include <TPRegexp.h>
#include <TGraphErrors.h>

#include <AliDCSSensor.h>
#include <AliGRPObject.h>
#include <AliESDpid.h>
#include <AliLog.h>
#include <AliESDEvent.h>
#include <AliExternalTrackParam.h>
#include <AliESDtrack.h>
#include <AliESDInputHandler.h>
#include <AliAnalysisManager.h>
#include <AliSplineFit.h>
#include <AliCDBId.h>
#include <AliCDBManager.h>
#include <AliCDBEntry.h>
#include <AliCDBRunRange.h>
#include <AliTender.h>
#include <AliTPCcalibDButil.h>
#include <AliPID.h>

#include "AliTPCTenderSupply.h"

ClassImp(AliTPCTenderSupply)

AliTPCTenderSupply::AliTPCTenderSupply() :
AliTenderSupply(),
fESDpid(0x0),
fGainNew(0x0),
fGainOld(0x0),
fGainAttachment(0x0),
fIsMC(kFALSE),
fGainCorrection(kTRUE),
fAttachmentCorrection(kFALSE),
fPcorrection(kFALSE),
fMultiCorrection(kFALSE),
fArrPidResponseMaster(0x0),
fMultiCorrMean(0x0),
fMultiCorrSigma(0x0),
fSpecificStorages(0x0),
fDebugLevel(0),
fMip(50),
fGRP(0x0),
fBeamType("PP"),
fLHCperiod(),
fMCperiod(),
fRecoPass(0)
{
  //
  // default ctor
  //
}

//_____________________________________________________
AliTPCTenderSupply::AliTPCTenderSupply(const char *name, const AliTender *tender) :
AliTenderSupply(name,tender),
fESDpid(0x0),
fGainNew(0x0),
fGainOld(0x0),
fGainAttachment(0x0),
fIsMC(kFALSE),
fGainCorrection(kTRUE),
fAttachmentCorrection(kFALSE),
fPcorrection(kFALSE),
fMultiCorrection(kFALSE),
fArrPidResponseMaster(0x0),
fMultiCorrMean(0x0),
fMultiCorrSigma(0x0),
fSpecificStorages(0x0),
fDebugLevel(0),
fMip(50),
fGRP(0x0),
fBeamType("PP"),
fLHCperiod(),
fMCperiod(),
fRecoPass(0)
{
  //
  // named ctor
  //
}

//_____________________________________________________
void AliTPCTenderSupply::Init()
{
  //
  // Initialise TPC tender
  //
  
  AliLog::SetClassDebugLevel("AliTPCTenderSupply",10);
  AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
  //
  // Setup PID object
  //
  fIsMC=mgr->GetMCtruthEventHandler();
  
  // Check if another detector already created the esd pid object
  // if not we create it and set it to the ESD input handler
  fESDpid=fTender->GetESDhandler()->GetESDpid();
  if (!fESDpid) {
    fESDpid=new AliESDpid(fIsMC);
    fTender->GetESDhandler()->SetESDpid(fESDpid);
  } 
  
  //setup specific storages
  if (fSpecificStorages){
    TNamed *storage;
    TIter nextStorage(fSpecificStorages);;
    while ( (storage=(TNamed*)nextStorage()) ){
      fTender->GetCDBManager()->SetSpecificStorage(storage->GetName(),storage->GetTitle());
      AliInfo(Form("Setting specific storage: %s (%s)",storage->GetName(), storage->GetTitle()));
    }
  }

  if (fIsMC){
    //force no gain and P correction in MC
    fGainCorrection=kFALSE;
    fPcorrection=kFALSE;
    fAttachmentCorrection=kFALSE;
  }
}

//_____________________________________________________
void AliTPCTenderSupply::ProcessEvent()
{
  //
  // Reapply pid information
  //
  
  AliESDEvent *event=fTender->GetEvent();
  if (!event) return;
  
  //load gain correction if run has changed
  if (fTender->RunChanged()){
    SetBeamType();
    SetRecoInfo();
    if ( fBeamType == "PBPB" ) fMultiCorrection=kTRUE;

    if (fDebugLevel>0) AliInfo(Form("Run Changed (%d)",fTender->GetRun()));
    SetParametrisation();
    if (fGainCorrection) SetSplines();
  }
  
  //
  // get gain correction factor
  //
  Double_t corrFactor = GetGainCorrection();
  Double_t corrAttachSlope = 0;
  Double_t corrGainMultiplicityPbPb=1;
  if (fAttachmentCorrection && fGainAttachment) corrAttachSlope = fGainAttachment->Eval(event->GetTimeStamp());
  if (fMultiCorrection&&fMultiCorrMean) corrGainMultiplicityPbPb = fMultiCorrMean->Eval(GetTPCMultiplicityBin());
  
  //
  // - correct TPC signals
  // - recalculate PID probabilities for TPC
  // - correct TPC signal multiplicity dependence

  Int_t ntracks=event->GetNumberOfTracks();
  for(Int_t itrack = 0; itrack < ntracks; itrack++){
    AliESDtrack *track=event->GetTrack(itrack);
    const AliExternalTrackParam *inner=track->GetInnerParam();
    
    // skip tracks without TPC information
    if (!inner) continue;

    //calculate total gain correction factor given by
    // o gain calibration factor
    // o attachment correction
    // o multiplicity correction in PbPb
    Float_t meanDrift= 250. - 0.5*TMath::Abs(2*inner->GetZ() + (247-83)*inner->GetTgl());
    Double_t corrGainTotal=corrFactor*(1 + corrAttachSlope*180.)/(1 + corrAttachSlope*meanDrift)/corrGainMultiplicityPbPb;

    // apply gain correction
    track->SetTPCsignal(track->GetTPCsignal()*corrGainTotal ,track->GetTPCsignalSigma(), track->GetTPCsignalN());

    // recalculate pid probabilities
    fESDpid->MakeTPCPID(track);
  }
}

//_____________________________________________________
Double_t AliTPCTenderSupply::GetTPCMultiplicityBin()
{
  //
  // Get TPC multiplicity in bins of 150
  //

  AliESDEvent *event=fTender->GetEvent();
  const AliESDVertex* vertexTPC = event->GetPrimaryVertexTPC();
  Double_t tpcMulti=0.;
  if(vertexTPC){
    Double_t vertexContribTPC=vertexTPC->GetNContributors();
    tpcMulti=vertexContribTPC/150.;
    if (tpcMulti>20.) tpcMulti=20.;
  }
  return tpcMulti;
}

//_____________________________________________________
Double_t AliTPCTenderSupply::GetMultiplicityCorrectionMean(Double_t tpcMulti)
{
  //
  // calculate correction factor for dEdx mean
  //
  Double_t meancorrection;
  if (fIsMC)
  {
    // MC data
    meancorrection=1.00054 + (0.00189566)*tpcMulti + (2.07777e-05)*tpcMulti*tpcMulti;
  }
  else
  {
    // real data
    meancorrection=0.999509 + (-0.00271488)*tpcMulti + (-2.98873e-06)*tpcMulti*tpcMulti;
  }
  
  return meancorrection;
  
}

//_____________________________________________________
Double_t AliTPCTenderSupply::GetMultiplicityCorrectionSigma(Double_t tpcMulti)
{
  //
  // calculate correction factor for dEdx sigma
  //
  Double_t sigmacorrection;
  if (fIsMC)
  {
    // MC data
    sigmacorrection=0.95972 + 0.0103721*tpcMulti;
  }
  else
  {
  // real data
    sigmacorrection=1.01817 + 0.0143673*tpcMulti;
  }
  return sigmacorrection;
  
}

//_____________________________________________________
void AliTPCTenderSupply::SetSplines()
{
  //
  // Get Gain splines from OCDB
  //
  
  AliInfo("Update Gain splines");
  
  //
  // Get GPR info for pressure correction
  //
  fPcorrection=kFALSE;
  
  AliCDBEntry *entryGRP=fTender->GetCDBManager()->Get("GRP/GRP/Data",fTender->GetRun());
  if (!entryGRP) {
    AliError("No new GRP entry found");
  } else {
    fGRP = (AliGRPObject*)entryGRP->GetObject();
  }
  if (fDebugLevel>1) AliInfo(Form("GRP entry used: %s",entryGRP->GetId().ToString().Data()));
  
  fGainNew=0x0;
  fGainOld=0x0;
  //
  //find previous entry from the UserInfo
  //
//   TTree *tree=((TChain*)fTender->GetInputData(0))->GetTree();
  AliAnalysisManager*mgr = AliAnalysisManager::GetAnalysisManager();
  AliAnalysisTaskSE *task = (AliAnalysisTaskSE*)mgr->GetTasks()->First();
  if (!task) return;
  TTree *tree=((TChain*)task->GetInputData(0))->GetTree();
  if (!tree) {
    AliError("Tree not found in ESDhandler");
    return;
  }
  
  TList *userInfo=(TList*)tree->GetUserInfo();
  if (!userInfo) {
    AliError("No UserInfo found in tree");
    return;
  }
  
  TList *cdbList=(TList*)userInfo->FindObject("cdbList");
  if (!cdbList) {
    AliError("No cdbList found in UserInfo");
    if (AliLog::GetGlobalLogLevel()>=AliLog::kError) userInfo->Print();
    return;
  }
  
  TIter nextCDB(cdbList);
  TObjString *os=0x0;
  while ( (os=(TObjString*)nextCDB()) ){
    if (!(os->GetString().Contains("TPC/Calib/TimeGain"))) continue;
    AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
    
    AliCDBEntry *entry=fTender->GetCDBManager()->Get(*id);
    if (!entry) {
      AliError("No previous gain calibration entry found");
      return;
    }
    
    if (fDebugLevel>1) AliInfo(Form("Used old Gain entry: %s",entry->GetId().ToString().Data()));
    
    TObjArray *arr=(TObjArray *)entry->GetObject();
    if (!arr) {
      AliError("Gain Splines array not found in calibration entry");
      return;
    }
    
    AliSplineFit *fit=(AliSplineFit*)arr->At(0);
    if (!fit) {
      AliError("Spline fit not found in array");
      return;
    }
    
    fGainOld = fit;
    delete id;
    break;
  }
  
  //
  //new gain correction
  //

  // This is in principle sill needed for the 2009 data and the 2010 b+c pass1 data
  //   however not supported any longer.
  // For the LHC10c pass2 there is an exception and we still load the splines, but a defined version
  // In case the attachment correction should be check again, this part of the code needs to be changed
  //   in order to load the gain entry again
  Bool_t special10cPass2=fLHCperiod=="LHC10C" && fRecoPass==2;
              
  AliCDBEntry *entryNew=0x0;
  if (special10cPass2) {
    entryNew=fTender->GetCDBManager()->Get("TPC/Calib/TimeGain",fTender->GetRun(),8);
  }
  if (!entryNew) {
    AliError("No new gain calibration entry found");
    return;
  }
  if (fDebugLevel>1) AliInfo(Form("Used new Gain entry: %s",entryNew->GetId().ToString().Data()));
  
  if (entryNew->GetId().GetLastRun()==AliCDBRunRange::Infinity()){
    if (fDebugLevel>0) AliInfo("Use P correction");
    fPcorrection=kTRUE;
  }
  
  TObjArray *arrSplines=(TObjArray *)entryNew->GetObject();
  if (!arrSplines) {
    AliError("Gain Splines array not found in new calibration entry");
    return;
  }
  
  fGainNew        = (AliSplineFit*)arrSplines->At(0);
  fGainAttachment = (TGraphErrors*)arrSplines->FindObject("TGRAPHERRORS_MEAN_ATTACHMENT_BEAM_ALL");
  
  if (!fGainNew) AliError("No recent spline fit object found");
}

//_____________________________________________________
Double_t AliTPCTenderSupply::GetGainCorrection()
{
  //
  // Calculate gain correction factor
  //
  AliESDEvent *event=fTender->GetEvent();
  UInt_t time=event->GetTimeStamp();
  
  Double_t gain=1;
  
  
  if (fGainOld){
    //TODO TODO TODO
    //first correction for the eval const problem
    // needs to be removed when the fix is in AliROOT and
    // the production was done with EvalGraphConst
    // This should be the case from V4-20-Rev11 on
    //TODO TODO TODO
    if ( fLHCperiod.Contains("LHC09") ||
         ((fLHCperiod=="LHC10B" || fLHCperiod=="LHC10C" || fLHCperiod=="LHC10D") && fRecoPass==2) ||
         (fLHCperiod=="LHC10E" && fRecoPass==1) ||
         (fLHCperiod=="LHC10H" && fRecoPass==1)
         ) {
      Double_t valDefault = fGainOld->Eval(time);
      Double_t valConst   = AliTPCcalibDButil::EvalGraphConst(fGainOld, time);
      gain = valDefault/valConst;
    }
    
    if (fGainNew){
      gain *= AliTPCcalibDButil::EvalGraphConst(fGainOld,time)/AliTPCcalibDButil::EvalGraphConst(fGainNew,time);
    }
  }
  
  //If there is only the default calibration, at least apply correction for pressure
  if (fPcorrection){
    if (fGRP) {
      Double_t pressure=fGRP->GetCavernAtmosPressure()->GetValue(time);
//       gain=fGainOld->Eval(time)/(7.03814-0.00459798*pressure)/49.53*fMip;
//       gain=fGainOld->Eval(time)/(7.03814-0.00459798*pressure)/51.51*fMip;
      gain=AliTPCcalibDButil::EvalGraphConst(fGainOld,time)/(7.03814-0.00459798*pressure)/51.51*fMip;
    }
  }
  
  return gain;
}

//_____________________________________________________
void AliTPCTenderSupply::SetBeamType()
{
  //
  // Set the beam type
  //

  fBeamType=fTender->GetEvent()->GetBeamType();
  if (fBeamType.IsNull()||fBeamType.Contains("No Beam")) fBeamType="p-p";
  fBeamType.ToUpper();
  fBeamType.ReplaceAll("-","");
}

//_____________________________________________________
void AliTPCTenderSupply::SetRecoInfo()
{
  //
  // Set reconstruction information
  //

  //reset information
  fRecoPass=0;
  fLHCperiod="";
  fMCperiod="";
  
  //Get the current file to check the reconstruction pass (UGLY, but not stored in ESD... )
  AliESDInputHandler *esdIH = dynamic_cast<AliESDInputHandler*> (fTender->GetESDhandler());
  if (!esdIH) return;
  TTree *tree= (TTree*)esdIH->GetTree();
  TFile *file= (TFile*)tree->GetCurrentFile();
  if (!file) {
    AliError("Current file not found, cannot set reconstruction information");
    return;
  }
  
  TString fileName(file->GetName());
  
  Int_t run=fTender->GetRun();
  
  
  TPRegexp reg(".*(LHC11[a-z]+[0-9]+[a-z]*)/.*");
  //find the period by run number (UGLY, but not stored in ESD... )
  if (run>=114737&&run<=117223)      { fLHCperiod="LHC10B"; fMCperiod="LHC10D1";  }
  else if (run>=118503&&run<=121040) { fLHCperiod="LHC10C"; fMCperiod="LHC10D1";  }
  else if (run>=122195&&run<=126437) { fLHCperiod="LHC10D"; fMCperiod="LHC10F6A"; }
  else if (run>=127719&&run<=130850) { fLHCperiod="LHC10E"; fMCperiod="LHC10F6A"; }
  else if (run>=133004&&run<=135029) { fLHCperiod="LHC10F"; fMCperiod="LHC10F6A"; }
  else if (run>=135654&&run<=136377) { fLHCperiod="LHC10G"; fMCperiod="LHC10F6A"; }
  else if (run>=136851&&run<=139517) { 
    fLHCperiod="LHC10H"; 
    fMCperiod="LHC10H8";  
    if (reg.MatchB(fileName)) fMCperiod="LHC11A10";
  }
  else if ( (run>=144871 && run <=146459 ) || ( run >=146686 && run<= 146860) ) {
   // low energy: 146686 - 146860
    fLHCperiod="LHC11A"; fMCperiod="LHC10F6A";
  }
  else if ( run>=148541  ){
    fLHCperiod="LHC11B"; fMCperiod="LHC10F6A";
  }
  
  //exception new pp MC productions from 2011
  if (fBeamType=="PP" && reg.MatchB(fileName)) fMCperiod="LHC11B2";
  
  //find pass from file name (UGLY, but not stored in ESD... )
  if (fileName.Contains("/pass1")) {
    fRecoPass=1;
  } else if (fileName.Contains("/pass2")) {
    fRecoPass=2;
  } else if (fileName.Contains("/pass3")) {
    fRecoPass=3;
  }
  
}

//_____________________________________________________
void AliTPCTenderSupply::SetParametrisation()
{
  //
  // Change PbPb multiplicity gain correction factor
  //
  
  if (fLHCperiod.IsNull()) {
    AliError("No period set, not changing parametrisation");
    return;
  }
  
  //Get CDB Entry with pid response parametrisations
  AliCDBEntry *pidCDB=fTender->GetCDBManager()->Get("TPC/Calib/PidResponse",fTender->GetRun());
  if (!fArrPidResponseMaster && pidCDB){
    fArrPidResponseMaster=dynamic_cast<TObjArray*>(pidCDB->GetObject());
    AliInfo(Form("Using pid response objects: %s",pidCDB->GetId().ToString().Data()));
  }

  if (!fArrPidResponseMaster){
    AliError("No valid PidResponse master found in OCDB");
    return;
  }
  //data type
  TString datatype="DATA";
  TString period=fLHCperiod;
  //in case of mc fRecoPass is per default 1
  if (fIsMC) {
    datatype="MC";
    fRecoPass=1;
    period=fMCperiod;
  }

  // Set PbPb correction
  fMultiCorrMean=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));

  if (fMultiCorrMean) AliInfo(Form("Setting multiplicity correction function: %s",fMultiCorrMean->GetName()));
  
}

//____________________________________________________________
void AliTPCTenderSupply::AddSpecificStorage(const char* cdbPath, const char* storage)
{
  //
  // Add a specific storage to be set up in Init
  //
  if (!fSpecificStorages) fSpecificStorages=new TObjArray;
  fSpecificStorages->Add(new TNamed(cdbPath, storage));
}
Commit	Line	Data
e75408ba	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16
	17	///////////////////////////////////////////////////////////////////////////////
	18	// //
	19	// TPC tender: reapply pid on the fly //
	20	// //
	21	///////////////////////////////////////////////////////////////////////////////
	22
	23	#include <TList.h>
ee981ab3	24	#include <TObjArray.h>
e75408ba	25	#include <TObjString.h>
e75408ba	26	#include <TChain.h>
ee981ab3	27	#include <TFile.h>
	28	#include <TString.h>
	29	#include <TPRegexp.h>
5dea385b	30	#include <TGraphErrors.h>
e75408ba	31
	32	#include <AliDCSSensor.h>
	33	#include <AliGRPObject.h>
	34	#include <AliESDpid.h>
	35	#include <AliLog.h>
	36	#include <AliESDEvent.h>
5dea385b	37	#include <AliExternalTrackParam.h>
e75408ba	38	#include <AliESDtrack.h>
	39	#include <AliESDInputHandler.h>
	40	#include <AliAnalysisManager.h>
	41	#include <AliSplineFit.h>
	42	#include <AliCDBId.h>
	43	#include <AliCDBManager.h>
	44	#include <AliCDBEntry.h>
f19a0c5e	45	#include <AliCDBRunRange.h>
e75408ba	46	#include <AliTender.h>
ee981ab3	47	#include <AliTPCcalibDButil.h>
ee981ab3	48	#include <AliPID.h>
e75408ba	49
	50	#include "AliTPCTenderSupply.h"
	51
e26aa0bb	52	ClassImp(AliTPCTenderSupply)
e75408ba	53
e75408ba	54	AliTPCTenderSupply::AliTPCTenderSupply() :
ee981ab3	55	AliTenderSupply(),
	56	fESDpid(0x0),
	57	fGainNew(0x0),
	58	fGainOld(0x0),
5dea385b	59	fGainAttachment(0x0),
5dea385b	60	fIsMC(kFALSE),
ee981ab3	61	fGainCorrection(kTRUE),
5dea385b	62	fAttachmentCorrection(kFALSE),
ee981ab3	63	fPcorrection(kFALSE),
5dea385b	64	fMultiCorrection(kFALSE),
ee981ab3	65	fArrPidResponseMaster(0x0),
5dea385b	66	fMultiCorrMean(0x0),
	67	fMultiCorrSigma(0x0),
	68	fSpecificStorages(0x0),
ee981ab3	69	fDebugLevel(0),
ee981ab3	70	fMip(50),
5dea385b	71	fGRP(0x0),
	72	fBeamType("PP"),
	73	fLHCperiod(),
	74	fMCperiod(),
	75	fRecoPass(0)
e75408ba	76	{
	77	//
	78	// default ctor
	79	//
	80	}
	81
	82	//_____________________________________________________
	83	AliTPCTenderSupply::AliTPCTenderSupply(const char name, const AliTender tender) :
ee981ab3	84	AliTenderSupply(name,tender),
	85	fESDpid(0x0),
	86	fGainNew(0x0),
	87	fGainOld(0x0),
5dea385b	88	fGainAttachment(0x0),
5dea385b	89	fIsMC(kFALSE),
ee981ab3	90	fGainCorrection(kTRUE),
5dea385b	91	fAttachmentCorrection(kFALSE),
ee981ab3	92	fPcorrection(kFALSE),
5dea385b	93	fMultiCorrection(kFALSE),
ee981ab3	94	fArrPidResponseMaster(0x0),
5dea385b	95	fMultiCorrMean(0x0),
	96	fMultiCorrSigma(0x0),
	97	fSpecificStorages(0x0),
ee981ab3	98	fDebugLevel(0),
ee981ab3	99	fMip(50),
5dea385b	100	fGRP(0x0),
	101	fBeamType("PP"),
	102	fLHCperiod(),
	103	fMCperiod(),
	104	fRecoPass(0)
e75408ba	105	{
	106	//
	107	// named ctor
	108	//
	109	}
	110
	111	//_____________________________________________________
	112	void AliTPCTenderSupply::Init()
	113	{
	114	//
	115	// Initialise TPC tender
	116	//
5dea385b	117
ee981ab3	118	AliLog::SetClassDebugLevel("AliTPCTenderSupply",10);
e75408ba	119	AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
	120	//
	121	// Setup PID object
	122	//
5dea385b	123	fIsMC=mgr->GetMCtruthEventHandler();
e75408ba	124
	125	// Check if another detector already created the esd pid object
	126	// if not we create it and set it to the ESD input handler
	127	fESDpid=fTender->GetESDhandler()->GetESDpid();
	128	if (!fESDpid) {
5dea385b	129	fESDpid=new AliESDpid(fIsMC);
e75408ba	130	fTender->GetESDhandler()->SetESDpid(fESDpid);
5dea385b	131	}
e75408ba	132
5dea385b	133	//setup specific storages
	134	if (fSpecificStorages){
	135	TNamed *storage;
	136	TIter nextStorage(fSpecificStorages);;
	137	while ( (storage=(TNamed*)nextStorage()) ){
	138	fTender->GetCDBManager()->SetSpecificStorage(storage->GetName(),storage->GetTitle());
	139	AliInfo(Form("Setting specific storage: %s (%s)",storage->GetName(), storage->GetTitle()));
	140	}
	141	}
	142
	143	if (fIsMC){
f19a0c5e	144	//force no gain and P correction in MC
e75408ba	145	fGainCorrection=kFALSE;
f19a0c5e	146	fPcorrection=kFALSE;
5dea385b	147	fAttachmentCorrection=kFALSE;
e75408ba	148	}
e75408ba	149	}
	150
	151	//_____________________________________________________
	152	void AliTPCTenderSupply::ProcessEvent()
	153	{
	154	//
	155	// Reapply pid information
	156	//
ee981ab3	157
e75408ba	158	AliESDEvent *event=fTender->GetEvent();
	159	if (!event) return;
	160
	161	//load gain correction if run has changed
	162	if (fTender->RunChanged()){
5dea385b	163	SetBeamType();
	164	SetRecoInfo();
	165	if ( fBeamType == "PBPB" ) fMultiCorrection=kTRUE;
	166
	167	if (fDebugLevel>0) AliInfo(Form("Run Changed (%d)",fTender->GetRun()));
ee981ab3	168	SetParametrisation();
e75408ba	169	if (fGainCorrection) SetSplines();
	170	}
	171
	172	//
	173	// get gain correction factor
	174	//
5dea385b	175	Double_t corrFactor = GetGainCorrection();
	176	Double_t corrAttachSlope = 0;
	177	Double_t corrGainMultiplicityPbPb=1;
	178	if (fAttachmentCorrection && fGainAttachment) corrAttachSlope = fGainAttachment->Eval(event->GetTimeStamp());
	179	if (fMultiCorrection&&fMultiCorrMean) corrGainMultiplicityPbPb = fMultiCorrMean->Eval(GetTPCMultiplicityBin());
e75408ba	180
	181	//
	182	// - correct TPC signals
	183	// - recalculate PID probabilities for TPC
5dea385b	184	// - correct TPC signal multiplicity dependence
5dea385b	185
e75408ba	186	Int_t ntracks=event->GetNumberOfTracks();
	187	for(Int_t itrack = 0; itrack < ntracks; itrack++){
	188	AliESDtrack *track=event->GetTrack(itrack);
5dea385b	189	const AliExternalTrackParam *inner=track->GetInnerParam();
	190
	191	// skip tracks without TPC information
	192	if (!inner) continue;
	193
	194	//calculate total gain correction factor given by
	195	// o gain calibration factor
	196	// o attachment correction
	197	// o multiplicity correction in PbPb
	198	Float_t meanDrift= 250. - 0.5TMath::Abs(2inner->GetZ() + (247-83)*inner->GetTgl());
	199	Double_t corrGainTotal=corrFactor(1 + corrAttachSlope180.)/(1 + corrAttachSlope*meanDrift)/corrGainMultiplicityPbPb;
	200
	201	// apply gain correction
	202	track->SetTPCsignal(track->GetTPCsignal()*corrGainTotal ,track->GetTPCsignalSigma(), track->GetTPCsignalN());
	203
	204	// recalculate pid probabilities
e75408ba	205	fESDpid->MakeTPCPID(track);
e75408ba	206	}
5dea385b	207	}
	208
	209	//_____________________________________________________
	210	Double_t AliTPCTenderSupply::GetTPCMultiplicityBin()
	211	{
	212	//
	213	// Get TPC multiplicity in bins of 150
	214	//
	215
	216	AliESDEvent *event=fTender->GetEvent();
	217	const AliESDVertex* vertexTPC = event->GetPrimaryVertexTPC();
	218	Double_t tpcMulti=0.;
	219	if(vertexTPC){
	220	Double_t vertexContribTPC=vertexTPC->GetNContributors();
	221	tpcMulti=vertexContribTPC/150.;
	222	if (tpcMulti>20.) tpcMulti=20.;
	223	}
	224	return tpcMulti;
	225	}
	226
	227	//_____________________________________________________
	228	Double_t AliTPCTenderSupply::GetMultiplicityCorrectionMean(Double_t tpcMulti)
	229	{
	230	//
	231	// calculate correction factor for dEdx mean
	232	//
	233	Double_t meancorrection;
	234	if (fIsMC)
	235	{
	236	// MC data
	237	meancorrection=1.00054 + (0.00189566)tpcMulti + (2.07777e-05)tpcMulti*tpcMulti;
	238	}
	239	else
	240	{
	241	// real data
	242	meancorrection=0.999509 + (-0.00271488)tpcMulti + (-2.98873e-06)tpcMulti*tpcMulti;
	243	}
	244
	245	return meancorrection;
	246
	247	}
	248
	249	//_____________________________________________________
	250	Double_t AliTPCTenderSupply::GetMultiplicityCorrectionSigma(Double_t tpcMulti)
	251	{
	252	//
	253	// calculate correction factor for dEdx sigma
	254	//
	255	Double_t sigmacorrection;
	256	if (fIsMC)
	257	{
	258	// MC data
	259	sigmacorrection=0.95972 + 0.0103721*tpcMulti;
	260	}
	261	else
	262	{
	263	// real data
	264	sigmacorrection=1.01817 + 0.0143673*tpcMulti;
	265	}
	266	return sigmacorrection;
e75408ba	267
	268	}
	269
	270	//_____________________________________________________
	271	void AliTPCTenderSupply::SetSplines()
	272	{
	273	//
	274	// Get Gain splines from OCDB
	275	//
ee981ab3	276
e75408ba	277	AliInfo("Update Gain splines");
ee981ab3	278
e75408ba	279	//
	280	// Get GPR info for pressure correction
	281	//
f19a0c5e	282	fPcorrection=kFALSE;
f19a0c5e	283
e75408ba	284	AliCDBEntry *entryGRP=fTender->GetCDBManager()->Get("GRP/GRP/Data",fTender->GetRun());
	285	if (!entryGRP) {
	286	AliError("No new GRP entry found");
	287	} else {
	288	fGRP = (AliGRPObject*)entryGRP->GetObject();
	289	}
5dea385b	290	if (fDebugLevel>1) AliInfo(Form("GRP entry used: %s",entryGRP->GetId().ToString().Data()));
e75408ba	291
	292	fGainNew=0x0;
	293	fGainOld=0x0;
	294	//
	295	//find previous entry from the UserInfo
	296	//
0e373bac	297	// TTree tree=((TChain)fTender->GetInputData(0))->GetTree();
	298	AliAnalysisManager*mgr = AliAnalysisManager::GetAnalysisManager();
	299	AliAnalysisTaskSE task = (AliAnalysisTaskSE)mgr->GetTasks()->First();
5dea385b	300	if (!task) return;
0e373bac	301	TTree tree=((TChain)task->GetInputData(0))->GetTree();
e75408ba	302	if (!tree) {
	303	AliError("Tree not found in ESDhandler");
	304	return;
	305	}
	306
	307	TList userInfo=(TList)tree->GetUserInfo();
	308	if (!userInfo) {
	309	AliError("No UserInfo found in tree");
	310	return;
	311	}
ee981ab3	312
e75408ba	313	TList cdbList=(TList)userInfo->FindObject("cdbList");
	314	if (!cdbList) {
	315	AliError("No cdbList found in UserInfo");
	316	if (AliLog::GetGlobalLogLevel()>=AliLog::kError) userInfo->Print();
	317	return;
	318	}
ee981ab3	319
e75408ba	320	TIter nextCDB(cdbList);
	321	TObjString *os=0x0;
	322	while ( (os=(TObjString*)nextCDB()) ){
	323	if (!(os->GetString().Contains("TPC/Calib/TimeGain"))) continue;
	324	AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
	325
	326	AliCDBEntry entry=fTender->GetCDBManager()->Get(id);
	327	if (!entry) {
	328	AliError("No previous gain calibration entry found");
	329	return;
	330	}
	331
5dea385b	332	if (fDebugLevel>1) AliInfo(Form("Used old Gain entry: %s",entry->GetId().ToString().Data()));
f19a0c5e	333
e75408ba	334	TObjArray arr=(TObjArray )entry->GetObject();
	335	if (!arr) {
	336	AliError("Gain Splines array not found in calibration entry");
	337	return;
	338	}
	339
	340	AliSplineFit fit=(AliSplineFit)arr->At(0);
	341	if (!fit) {
	342	AliError("Spline fit not found in array");
	343	return;
	344	}
	345
	346	fGainOld = fit;
	347	delete id;
	348	break;
	349	}
ee981ab3	350
e75408ba	351	//
	352	//new gain correction
	353	//
5dea385b	354
	355	// This is in principle sill needed for the 2009 data and the 2010 b+c pass1 data
	356	// however not supported any longer.
	357	// For the LHC10c pass2 there is an exception and we still load the splines, but a defined version
	358	// In case the attachment correction should be check again, this part of the code needs to be changed
	359	// in order to load the gain entry again
	360	Bool_t special10cPass2=fLHCperiod=="LHC10C" && fRecoPass==2;
	361
	362	AliCDBEntry *entryNew=0x0;
	363	if (special10cPass2) {
	364	entryNew=fTender->GetCDBManager()->Get("TPC/Calib/TimeGain",fTender->GetRun(),8);
	365	}
e75408ba	366	if (!entryNew) {
	367	AliError("No new gain calibration entry found");
	368	return;
	369	}
5dea385b	370	if (fDebugLevel>1) AliInfo(Form("Used new Gain entry: %s",entryNew->GetId().ToString().Data()));
f19a0c5e	371
f19a0c5e	372	if (entryNew->GetId().GetLastRun()==AliCDBRunRange::Infinity()){
5dea385b	373	if (fDebugLevel>0) AliInfo("Use P correction");
f19a0c5e	374	fPcorrection=kTRUE;
f19a0c5e	375	}
e75408ba	376
	377	TObjArray arrSplines=(TObjArray )entryNew->GetObject();
	378	if (!arrSplines) {
	379	AliError("Gain Splines array not found in new calibration entry");
	380	return;
	381	}
	382
5dea385b	383	fGainNew = (AliSplineFit*)arrSplines->At(0);
5dea385b	384	fGainAttachment = (TGraphErrors*)arrSplines->FindObject("TGRAPHERRORS_MEAN_ATTACHMENT_BEAM_ALL");
e75408ba	385
	386	if (!fGainNew) AliError("No recent spline fit object found");
	387	}
	388
	389	//_____________________________________________________
	390	Double_t AliTPCTenderSupply::GetGainCorrection()
	391	{
	392	//
	393	// Calculate gain correction factor
	394	//
	395	AliESDEvent *event=fTender->GetEvent();
	396	UInt_t time=event->GetTimeStamp();
5dea385b	397
ee981ab3	398	Double_t gain=1;
e75408ba	399
ee981ab3	400
	401	if (fGainOld){
	402	//TODO TODO TODO
5dea385b	403	//first correction for the eval const problem
ee981ab3	404	// needs to be removed when the fix is in AliROOT and
ee981ab3	405	// the production was done with EvalGraphConst
5dea385b	406	// This should be the case from V4-20-Rev11 on
ee981ab3	407	//TODO TODO TODO
5dea385b	408	if ( fLHCperiod.Contains("LHC09") \|\|
	409	((fLHCperiod=="LHC10B" \|\| fLHCperiod=="LHC10C" \|\| fLHCperiod=="LHC10D") && fRecoPass==2) \|\|
	410	(fLHCperiod=="LHC10E" && fRecoPass==1) \|\|
	411	(fLHCperiod=="LHC10H" && fRecoPass==1)
	412	) {
ee981ab3	413	Double_t valDefault = fGainOld->Eval(time);
	414	Double_t valConst = AliTPCcalibDButil::EvalGraphConst(fGainOld, time);
	415	gain = valDefault/valConst;
	416	}
5dea385b	417
ee981ab3	418	if (fGainNew){
	419	gain *= AliTPCcalibDButil::EvalGraphConst(fGainOld,time)/AliTPCcalibDButil::EvalGraphConst(fGainNew,time);
	420	}
	421	}
e75408ba	422
f19a0c5e	423	//If there is only the default calibration, at least apply correction for pressure
f19a0c5e	424	if (fPcorrection){
e75408ba	425	if (fGRP) {
e75408ba	426	Double_t pressure=fGRP->GetCavernAtmosPressure()->GetValue(time);
ee981ab3	427	// gain=fGainOld->Eval(time)/(7.03814-0.00459798pressure)/49.53fMip;
	428	// gain=fGainOld->Eval(time)/(7.03814-0.00459798pressure)/51.51fMip;
	429	gain=AliTPCcalibDButil::EvalGraphConst(fGainOld,time)/(7.03814-0.00459798pressure)/51.51fMip;
e75408ba	430	}
	431	}
	432
	433	return gain;
	434	}
	435
ee981ab3	436	//_____________________________________________________
5dea385b	437	void AliTPCTenderSupply::SetBeamType()
ee981ab3	438	{
ee981ab3	439	//
5dea385b	440	// Set the beam type
ee981ab3	441	//
ee981ab3	442
5dea385b	443	fBeamType=fTender->GetEvent()->GetBeamType();
	444	if (fBeamType.IsNull()\|\|fBeamType.Contains("No Beam")) fBeamType="p-p";
	445	fBeamType.ToUpper();
	446	fBeamType.ReplaceAll("-","");
	447	}
	448
	449	//_____________________________________________________
	450	void AliTPCTenderSupply::SetRecoInfo()
	451	{
	452	//
	453	// Set reconstruction information
	454	//
ee981ab3	455
5dea385b	456	//reset information
	457	fRecoPass=0;
	458	fLHCperiod="";
	459	fMCperiod="";
	460
ee981ab3	461	//Get the current file to check the reconstruction pass (UGLY, but not stored in ESD... )
0e373bac	462	AliESDInputHandler esdIH = dynamic_cast<AliESDInputHandler> (fTender->GetESDhandler());
	463	if (!esdIH) return;
	464	TTree tree= (TTree)esdIH->GetTree();
5dea385b	465	TFile file= (TFile)tree->GetCurrentFile();
ee981ab3	466	if (!file) {
5dea385b	467	AliError("Current file not found, cannot set reconstruction information");
ee981ab3	468	return;
ee981ab3	469	}
5dea385b	470
f25dc3ef	471	TString fileName(file->GetName());
f25dc3ef	472
ee981ab3	473	Int_t run=fTender->GetRun();
5dea385b	474
	475
	476	TPRegexp reg(".(LHC11[a-z]+[0-9]+[a-z])/.*");
ee981ab3	477	//find the period by run number (UGLY, but not stored in ESD... )
5dea385b	478	if (run>=114737&&run<=117223) { fLHCperiod="LHC10B"; fMCperiod="LHC10D1"; }
	479	else if (run>=118503&&run<=121040) { fLHCperiod="LHC10C"; fMCperiod="LHC10D1"; }
	480	else if (run>=122195&&run<=126437) { fLHCperiod="LHC10D"; fMCperiod="LHC10F6A"; }
	481	else if (run>=127719&&run<=130850) { fLHCperiod="LHC10E"; fMCperiod="LHC10F6A"; }
	482	else if (run>=133004&&run<=135029) { fLHCperiod="LHC10F"; fMCperiod="LHC10F6A"; }
	483	else if (run>=135654&&run<=136377) { fLHCperiod="LHC10G"; fMCperiod="LHC10F6A"; }
	484	else if (run>=136851&&run<=139517) {
	485	fLHCperiod="LHC10H";
	486	fMCperiod="LHC10H8";
	487	if (reg.MatchB(fileName)) fMCperiod="LHC11A10";
	488	}
	489	else if ( (run>=144871 && run <=146459 ) \|\| ( run >=146686 && run<= 146860) ) {
	490	// low energy: 146686 - 146860
	491	fLHCperiod="LHC11A"; fMCperiod="LHC10F6A";
	492	}
	493	else if ( run>=148541 ){
	494	fLHCperiod="LHC11B"; fMCperiod="LHC10F6A";
	495	}
	496
	497	//exception new pp MC productions from 2011
	498	if (fBeamType=="PP" && reg.MatchB(fileName)) fMCperiod="LHC11B2";
ee981ab3	499
ee981ab3	500	//find pass from file name (UGLY, but not stored in ESD... )
ee981ab3	501	if (fileName.Contains("/pass1")) {
5dea385b	502	fRecoPass=1;
ee981ab3	503	} else if (fileName.Contains("/pass2")) {
5dea385b	504	fRecoPass=2;
	505	} else if (fileName.Contains("/pass3")) {
	506	fRecoPass=3;
ee981ab3	507	}
ee981ab3	508
5dea385b	509	}
	510
	511	//_____________________________________________________
	512	void AliTPCTenderSupply::SetParametrisation()
	513	{
ee981ab3	514	//
5dea385b	515	// Change PbPb multiplicity gain correction factor
ee981ab3	516	//
5dea385b	517
	518	if (fLHCperiod.IsNull()) {
	519	AliError("No period set, not changing parametrisation");
	520	return;
ee981ab3	521	}
ee981ab3	522
5dea385b	523	//Get CDB Entry with pid response parametrisations
	524	AliCDBEntry *pidCDB=fTender->GetCDBManager()->Get("TPC/Calib/PidResponse",fTender->GetRun());
	525	if (!fArrPidResponseMaster && pidCDB){
	526	fArrPidResponseMaster=dynamic_cast<TObjArray*>(pidCDB->GetObject());
	527	AliInfo(Form("Using pid response objects: %s",pidCDB->GetId().ToString().Data()));
	528	}
ee981ab3	529
f25dc3ef	530	if (!fArrPidResponseMaster){
	531	AliError("No valid PidResponse master found in OCDB");
	532	return;
	533	}
ee981ab3	534	//data type
ee981ab3	535	TString datatype="DATA";
5dea385b	536	TString period=fLHCperiod;
	537	//in case of mc fRecoPass is per default 1
	538	if (fIsMC) {
ee981ab3	539	datatype="MC";
5dea385b	540	fRecoPass=1;
5dea385b	541	period=fMCperiod;
ee981ab3	542	}
ee981ab3	543
5dea385b	544	// Set PbPb correction
5dea385b	545	fMultiCorrMean=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
ee981ab3	546
5dea385b	547	if (fMultiCorrMean) AliInfo(Form("Setting multiplicity correction function: %s",fMultiCorrMean->GetName()));
	548
	549	}
ee981ab3	550
5dea385b	551	//____________________________________________________________
	552	void AliTPCTenderSupply::AddSpecificStorage(const char* cdbPath, const char* storage)
	553	{
	554	//
	555	// Add a specific storage to be set up in Init
	556	//
	557	if (!fSpecificStorages) fSpecificStorages=new TObjArray;
	558	fSpecificStorages->Add(new TNamed(cdbPath, storage));
ee981ab3	559	}
ee981ab3	560