[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 <AliDCSSensor.h>
#include <AliGRPObject.h>
#include <AliESDpid.h>
#include <AliLog.h>
#include <AliESDEvent.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"


AliTPCTenderSupply::AliTPCTenderSupply() :
AliTenderSupply(),
fESDpid(0x0),
fGainNew(0x0),
fGainOld(0x0),
fGainCorrection(kTRUE),
fPcorrection(kFALSE),
fArrPidResponseMaster(0x0),
fDebugLevel(0),
fMip(50),
fGRP(0x0)
{
  //
  // default ctor
  //
}

//_____________________________________________________
AliTPCTenderSupply::AliTPCTenderSupply(const char *name, const AliTender *tender) :
AliTenderSupply(name,tender),
fESDpid(0x0),
fGainNew(0x0),
fGainOld(0x0),
fGainCorrection(kTRUE),
fPcorrection(kFALSE),
fArrPidResponseMaster(0x0),
fDebugLevel(0),
fMip(50),
fGRP(0x0)
{
  //
  // named ctor
  //
}

//_____________________________________________________
void AliTPCTenderSupply::Init()
{
  //
  // Initialise TPC tender
  //

  AliLog::SetClassDebugLevel("AliTPCTenderSupply",10);
  AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
  //
  // Setup PID object
  //
  
  // 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;
    fTender->GetESDhandler()->SetESDpid(fESDpid);
  }
  
  //
  //set bethe bloch parameters depending on whether we have MC or real data
  //
  // for the moment we set the values hardwired. In future they should be stored either in
  // the OCDB or an equivalent calibration data base
  //
  Double_t alephParameters[5];
  // simulation
  alephParameters[0] = 2.15898e+00/50.;
  alephParameters[1] = 1.75295e+01;
  alephParameters[2] = 3.40030e-09;
  alephParameters[3] = 1.96178e+00;
  alephParameters[4] = 3.91720e+00;
  
    // assume data if there is no mc handler
  if (!mgr->GetMCtruthEventHandler()){
    alephParameters[0] = 0.0283086/0.97;
    //alephParameters[0] = 0.0283086;
    alephParameters[1] = 2.63394e+01;
    alephParameters[2] = 5.04114e-11;
    alephParameters[3] = 2.12543e+00;
    alephParameters[4] = 4.88663e+00;
    //temporary solution
    //fESDpid->GetTPCResponse().SetMip(fMip);
    if (fDebugLevel>0) AliInfo(Form("Use Data parametrisation, Mip set to: %.3f\n",fMip));
    //fESDpid->GetTPCResponse().SetMip(49.2);
  } else {
    //force no gain and P correction in MC
    fGainCorrection=kFALSE;
    fPcorrection=kFALSE;
    if (fDebugLevel>0) AliInfo("Use MC parametrisation\n");
  }
  
  fESDpid->GetTPCResponse().SetBetheBlochParameters(
    alephParameters[0],alephParameters[1],alephParameters[2],
    alephParameters[3],alephParameters[4]);
  
  //set detector resolution parametrisation
  fESDpid->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04);
}

//_____________________________________________________
void AliTPCTenderSupply::ProcessEvent()
{
  //
  // Reapply pid information
  //
  
  AliESDEvent *event=fTender->GetEvent();
  if (!event) return;
  
  //load gain correction if run has changed
  if (fTender->RunChanged()){
    if (fDebugLevel>0) AliInfo(Form("Run Changed (%d)\n",fTender->GetRun()));
    SetParametrisation();
    if (fGainCorrection) SetSplines();
  }
  
  //
  // get gain correction factor
  //
  Double_t corrFactor = 1;
  if (fGainCorrection) corrFactor=GetGainCorrection();
  
  //
  // - correct TPC signals
  // - recalculate PID probabilities for TPC
  //
  Int_t ntracks=event->GetNumberOfTracks();
  for(Int_t itrack = 0; itrack < ntracks; itrack++){
    AliESDtrack *track=event->GetTrack(itrack);
    if (fGainCorrection)
      track->SetTPCsignal(track->GetTPCsignal()*corrFactor,track->GetTPCsignalSigma(),track->GetTPCsignalN());
    fESDpid->MakeTPCPID(track);
  }
  
}

//_____________________________________________________
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\n",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();
  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\n",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
  //
  AliCDBEntry *entryNew=fTender->GetCDBManager()->Get("TPC/Calib/TimeGain",fTender->GetRun());
  if (!entryNew) {
    AliError("No new gain calibration entry found");
    return;
  }
  if (fDebugLevel>1) AliInfo(Form("Used new Gain entry: %s\n",entryNew->GetId().ToString().Data()));
  
  if (entryNew->GetId().GetLastRun()==AliCDBRunRange::Infinity()){
    if (fDebugLevel>0) AliInfo("Use P correction\n");
    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);
  
  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
    // needs to be removed when the fix is in AliROOT and
    // the production was done with EvalGraphConst
    //TODO TODO TODO
    if (fTender->GetRun()<=138154||fTender->GetRun()==138197){
      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::SetParametrisation()
{
  //
  // Change BB parametrisation for current run
  //

  //Get CDB Entry with pid response parametrisations
  AliCDBEntry *pidCDB=fTender->GetCDBManager()->Get("TPC/Calib/PidResponse",fTender->GetRun());
  if (pidCDB){
    fArrPidResponseMaster=(TObjArray*)pidCDB->GetObject();
  }

  //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("File not found, not changing parametrisation");
    return;
  }

  Int_t run=fTender->GetRun();
  AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
  Double_t alephParameters[5]={0,0,0,0,0};

  //find the period by run number (UGLY, but not stored in ESD... )
  TString period;
  if (run>=114737&&run<=117223) period="LHC10B";
  else if (run>=118503&&run<=121040) period="LHC10C";
  else if (run>=122195&&run<=126437) period="LHC10D";
  else if (run>=127719&&run<=130850) period="LHC10E";
  else if (run>=133004&&run<=135029) period="LHC10F";
  else if (run>=135654&&run<=136377) period="LHC10G";
  else if (run>=136851&&run<=139517) period="LHC10H";
  else if (run>=139699) period="LHC11A";
  
  //find pass from file name (UGLY, but not stored in ESD... )
  TString fileName(file->GetName());
  Int_t pass=0;
  if (fileName.Contains("/pass1")) {
    pass=1;
  } else if (fileName.Contains("/pass2")) {
    pass=2;
  }

  //beam type
  TString beamtype=fTender->GetEvent()->GetBeamType();
  if (beamtype.IsNull()||beamtype.Contains("No Beam")) beamtype="p-p";
  beamtype.ToUpper();
  beamtype.ReplaceAll("-","");
  
  //
  // Set default parametrisations for data and MC
  //
  Bool_t isMC=mgr->GetMCtruthEventHandler();
  if (isMC){
    //MC data
    alephParameters[0] = 2.15898e+00/50.;
    alephParameters[1] = 1.75295e+01;
    alephParameters[2] = 3.40030e-09;
    alephParameters[3] = 1.96178e+00;
    alephParameters[4] = 3.91720e+00;
    
  } else {
    //data parametrisation

    
    //use defaut data parametrisation in case no other will be selected
    alephParameters[0] = 0.0283086/0.97;
    alephParameters[1] = 2.63394e+01;
    alephParameters[2] = 5.04114e-11;
    alephParameters[3] = 2.12543e+00;
    alephParameters[4] = 4.88663e+00;

    fESDpid->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04);
    
    if (pass==1){
      
    } else if (pass==2){
      //find period
      if (run>=114737&&run<=117223){
        //LHC10b
      } else if (run>=118503&&run<=121040) {
        //LHC10c
      } else if (run>=122195){
        //LHC10d +
        // last run in LHC10d: &&run<126437
        alephParameters[0] = 1.63246/50.;
        alephParameters[1] = 2.20028e+01;
        alephParameters[2] = TMath::Exp(-2.48879e+01);
        alephParameters[3] = 2.39804e+00;
        alephParameters[4] = 5.12090e+00;
      
        //
        fESDpid->GetTPCResponse().SetSigma(2.30176e-02, 5.60422e+02);
      }
    }
    
    if ( beamtype == "PBPB" ){
      AliInfo("BETHE-BLOCH parametrization for PbPb !!!!!!!!!!!!!!!!!!!!!!\n");

      alephParameters[0] = 1.25202/50.;   //was 1.79571/55.;
      alephParameters[1] = 2.74992e+01;   //was 22.0028;
      alephParameters[2] = TMath::Exp(-3.31517e+01);  //was1.55354e-11;
      alephParameters[3] = 2.46246;       //was 2.39804; 
      alephParameters[4] = 6.78938;       //was 5.1209;
    }
  }
  
  fESDpid->GetTPCResponse().SetBetheBlochParameters(
    alephParameters[0],alephParameters[1],alephParameters[2],
    alephParameters[3],alephParameters[4]);

  //data type
  TString datatype="DATA";
  //in case of mc pass is per default 1
  if (isMC) {
    datatype="MC";
    pass=1;
  }

  //
  //set the new parametrisation
  //

  if (fArrPidResponseMaster){
    TObject *grAll=0x0;
    //for MC don't use period information
    if (isMC) period="[A-Z0-9]*";
    //pattern for the default entry (valid for all particles)
    TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),pass,beamtype.Data()));
    
    //loop over entries and filter them
    for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
      TObject *responseFunction=fArrPidResponseMaster->At(iresp);
      TString responseName=responseFunction->GetName();
      
      if (!reg.MatchB(responseName)) continue;
      
      TObjArray *arr=reg.MatchS(responseName);
      TString particleName=arr->At(1)->GetName();
      delete arr;
      if (particleName.IsNull()) continue;
      if (particleName=="ALL") grAll=responseFunction;
      else {
        //find particle id
        for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
          TString particle=AliPID::ParticleName(ispec);
          particle.ToUpper();
          if ( particle == particleName ){
            //test if there is already a function set. If yes, cleanup
            TObject *old=const_cast<TObject*>(fESDpid->GetTPCResponse().GetResponseFunction((AliPID::EParticleType)ispec));
            if (old) delete old;
            fESDpid->GetTPCResponse().SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
            fESDpid->GetTPCResponse().SetUseDatabase(kTRUE);
            AliInfo(Form("Adding graph: %d - %s\n",ispec,responseFunction->GetName()));
            break;
          }
        }
      }
    }

    //set default response function to all particles which don't have a specific one
    if (grAll){
      for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
        if (!fESDpid->GetTPCResponse().GetResponseFunction((AliPID::EParticleType)ispec)){
          fESDpid->GetTPCResponse().SetResponseFunction((AliPID::EParticleType)ispec,grAll);
          AliInfo(Form("Adding graph: %d - %s\n",ispec,grAll->GetName()));
        }
      }
    }
  }
}
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>
e75408ba	30
	31	#include <AliDCSSensor.h>
	32	#include <AliGRPObject.h>
	33	#include <AliESDpid.h>
	34	#include <AliLog.h>
	35	#include <AliESDEvent.h>
	36	#include <AliESDtrack.h>
	37	#include <AliESDInputHandler.h>
	38	#include <AliAnalysisManager.h>
	39	#include <AliSplineFit.h>
	40	#include <AliCDBId.h>
	41	#include <AliCDBManager.h>
	42	#include <AliCDBEntry.h>
f19a0c5e	43	#include <AliCDBRunRange.h>
e75408ba	44	#include <AliTender.h>
ee981ab3	45	#include <AliTPCcalibDButil.h>
ee981ab3	46	#include <AliPID.h>
e75408ba	47
	48	#include "AliTPCTenderSupply.h"
	49
	50
	51	AliTPCTenderSupply::AliTPCTenderSupply() :
ee981ab3	52	AliTenderSupply(),
	53	fESDpid(0x0),
	54	fGainNew(0x0),
	55	fGainOld(0x0),
	56	fGainCorrection(kTRUE),
	57	fPcorrection(kFALSE),
	58	fArrPidResponseMaster(0x0),
	59	fDebugLevel(0),
	60	fMip(50),
	61	fGRP(0x0)
e75408ba	62	{
	63	//
	64	// default ctor
	65	//
	66	}
	67
	68	//_____________________________________________________
	69	AliTPCTenderSupply::AliTPCTenderSupply(const char name, const AliTender tender) :
ee981ab3	70	AliTenderSupply(name,tender),
	71	fESDpid(0x0),
	72	fGainNew(0x0),
	73	fGainOld(0x0),
	74	fGainCorrection(kTRUE),
	75	fPcorrection(kFALSE),
	76	fArrPidResponseMaster(0x0),
	77	fDebugLevel(0),
	78	fMip(50),
	79	fGRP(0x0)
e75408ba	80	{
	81	//
	82	// named ctor
	83	//
	84	}
	85
	86	//_____________________________________________________
	87	void AliTPCTenderSupply::Init()
	88	{
	89	//
	90	// Initialise TPC tender
	91	//
	92
ee981ab3	93	AliLog::SetClassDebugLevel("AliTPCTenderSupply",10);
e75408ba	94	AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
	95	//
	96	// Setup PID object
	97	//
	98
	99	// Check if another detector already created the esd pid object
	100	// if not we create it and set it to the ESD input handler
	101	fESDpid=fTender->GetESDhandler()->GetESDpid();
	102	if (!fESDpid) {
	103	fESDpid=new AliESDpid;
	104	fTender->GetESDhandler()->SetESDpid(fESDpid);
	105	}
	106
	107	//
	108	//set bethe bloch parameters depending on whether we have MC or real data
	109	//
	110	// for the moment we set the values hardwired. In future they should be stored either in
	111	// the OCDB or an equivalent calibration data base
	112	//
	113	Double_t alephParameters[5];
	114	// simulation
	115	alephParameters[0] = 2.15898e+00/50.;
	116	alephParameters[1] = 1.75295e+01;
	117	alephParameters[2] = 3.40030e-09;
	118	alephParameters[3] = 1.96178e+00;
	119	alephParameters[4] = 3.91720e+00;
	120
	121	// assume data if there is no mc handler
	122	if (!mgr->GetMCtruthEventHandler()){
	123	alephParameters[0] = 0.0283086/0.97;
	124	//alephParameters[0] = 0.0283086;
	125	alephParameters[1] = 2.63394e+01;
	126	alephParameters[2] = 5.04114e-11;
	127	alephParameters[3] = 2.12543e+00;
	128	alephParameters[4] = 4.88663e+00;
	129	//temporary solution
ee981ab3	130	//fESDpid->GetTPCResponse().SetMip(fMip);
ee981ab3	131	if (fDebugLevel>0) AliInfo(Form("Use Data parametrisation, Mip set to: %.3f\n",fMip));
e75408ba	132	//fESDpid->GetTPCResponse().SetMip(49.2);
e75408ba	133	} else {
f19a0c5e	134	//force no gain and P correction in MC
e75408ba	135	fGainCorrection=kFALSE;
f19a0c5e	136	fPcorrection=kFALSE;
ee981ab3	137	if (fDebugLevel>0) AliInfo("Use MC parametrisation\n");
e75408ba	138	}
	139
	140	fESDpid->GetTPCResponse().SetBetheBlochParameters(
	141	alephParameters[0],alephParameters[1],alephParameters[2],
	142	alephParameters[3],alephParameters[4]);
ee981ab3	143
e75408ba	144	//set detector resolution parametrisation
	145	fESDpid->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04);
	146	}
	147
	148	//_____________________________________________________
	149	void AliTPCTenderSupply::ProcessEvent()
	150	{
	151	//
	152	// Reapply pid information
	153	//
ee981ab3	154
e75408ba	155	AliESDEvent *event=fTender->GetEvent();
	156	if (!event) return;
	157
	158	//load gain correction if run has changed
	159	if (fTender->RunChanged()){
ee981ab3	160	if (fDebugLevel>0) AliInfo(Form("Run Changed (%d)\n",fTender->GetRun()));
ee981ab3	161	SetParametrisation();
e75408ba	162	if (fGainCorrection) SetSplines();
	163	}
	164
	165	//
	166	// get gain correction factor
	167	//
ee981ab3	168	Double_t corrFactor = 1;
ee981ab3	169	if (fGainCorrection) corrFactor=GetGainCorrection();
e75408ba	170
	171	//
	172	// - correct TPC signals
	173	// - recalculate PID probabilities for TPC
	174	//
	175	Int_t ntracks=event->GetNumberOfTracks();
	176	for(Int_t itrack = 0; itrack < ntracks; itrack++){
	177	AliESDtrack *track=event->GetTrack(itrack);
	178	if (fGainCorrection)
	179	track->SetTPCsignal(track->GetTPCsignal()*corrFactor,track->GetTPCsignalSigma(),track->GetTPCsignalN());
	180	fESDpid->MakeTPCPID(track);
	181	}
	182
	183	}
	184
	185	//_____________________________________________________
	186	void AliTPCTenderSupply::SetSplines()
	187	{
	188	//
	189	// Get Gain splines from OCDB
	190	//
ee981ab3	191
e75408ba	192	AliInfo("Update Gain splines");
ee981ab3	193
e75408ba	194	//
	195	// Get GPR info for pressure correction
	196	//
f19a0c5e	197	fPcorrection=kFALSE;
f19a0c5e	198
e75408ba	199	AliCDBEntry *entryGRP=fTender->GetCDBManager()->Get("GRP/GRP/Data",fTender->GetRun());
	200	if (!entryGRP) {
	201	AliError("No new GRP entry found");
	202	} else {
	203	fGRP = (AliGRPObject*)entryGRP->GetObject();
	204	}
ee981ab3	205	if (fDebugLevel>1) AliInfo(Form("GRP entry used: %s\n",entryGRP->GetId().ToString().Data()));
e75408ba	206
	207	fGainNew=0x0;
	208	fGainOld=0x0;
	209	//
	210	//find previous entry from the UserInfo
	211	//
0e373bac	212	// TTree tree=((TChain)fTender->GetInputData(0))->GetTree();
	213	AliAnalysisManager*mgr = AliAnalysisManager::GetAnalysisManager();
	214	AliAnalysisTaskSE task = (AliAnalysisTaskSE)mgr->GetTasks()->First();
	215	TTree tree=((TChain)task->GetInputData(0))->GetTree();
e75408ba	216	if (!tree) {
	217	AliError("Tree not found in ESDhandler");
	218	return;
	219	}
	220
	221	TList userInfo=(TList)tree->GetUserInfo();
	222	if (!userInfo) {
	223	AliError("No UserInfo found in tree");
	224	return;
	225	}
ee981ab3	226
e75408ba	227	TList cdbList=(TList)userInfo->FindObject("cdbList");
	228	if (!cdbList) {
	229	AliError("No cdbList found in UserInfo");
	230	if (AliLog::GetGlobalLogLevel()>=AliLog::kError) userInfo->Print();
	231	return;
	232	}
ee981ab3	233
e75408ba	234	TIter nextCDB(cdbList);
	235	TObjString *os=0x0;
	236	while ( (os=(TObjString*)nextCDB()) ){
	237	if (!(os->GetString().Contains("TPC/Calib/TimeGain"))) continue;
	238	AliCDBId *id=AliCDBId::MakeFromString(os->GetString());
	239
	240	AliCDBEntry entry=fTender->GetCDBManager()->Get(id);
	241	if (!entry) {
	242	AliError("No previous gain calibration entry found");
	243	return;
	244	}
	245
ee981ab3	246	if (fDebugLevel>1) AliInfo(Form("Used old Gain entry: %s\n",entry->GetId().ToString().Data()));
f19a0c5e	247
e75408ba	248	TObjArray arr=(TObjArray )entry->GetObject();
	249	if (!arr) {
	250	AliError("Gain Splines array not found in calibration entry");
	251	return;
	252	}
	253
	254	AliSplineFit fit=(AliSplineFit)arr->At(0);
	255	if (!fit) {
	256	AliError("Spline fit not found in array");
	257	return;
	258	}
	259
	260	fGainOld = fit;
	261	delete id;
	262	break;
	263	}
ee981ab3	264
e75408ba	265	//
	266	//new gain correction
	267	//
	268	AliCDBEntry *entryNew=fTender->GetCDBManager()->Get("TPC/Calib/TimeGain",fTender->GetRun());
	269	if (!entryNew) {
	270	AliError("No new gain calibration entry found");
	271	return;
	272	}
ee981ab3	273	if (fDebugLevel>1) AliInfo(Form("Used new Gain entry: %s\n",entryNew->GetId().ToString().Data()));
f19a0c5e	274
f19a0c5e	275	if (entryNew->GetId().GetLastRun()==AliCDBRunRange::Infinity()){
ee981ab3	276	if (fDebugLevel>0) AliInfo("Use P correction\n");
f19a0c5e	277	fPcorrection=kTRUE;
f19a0c5e	278	}
e75408ba	279
	280	TObjArray arrSplines=(TObjArray )entryNew->GetObject();
	281	if (!arrSplines) {
	282	AliError("Gain Splines array not found in new calibration entry");
	283	return;
	284	}
	285
	286	fGainNew = (AliSplineFit*)arrSplines->At(0);
	287
	288	if (!fGainNew) AliError("No recent spline fit object found");
	289	}
	290
	291	//_____________________________________________________
	292	Double_t AliTPCTenderSupply::GetGainCorrection()
	293	{
	294	//
	295	// Calculate gain correction factor
	296	//
	297	AliESDEvent *event=fTender->GetEvent();
	298	UInt_t time=event->GetTimeStamp();
ee981ab3	299
ee981ab3	300	Double_t gain=1;
e75408ba	301
ee981ab3	302
	303	if (fGainOld){
	304	//TODO TODO TODO
	305	//first correction for the eval
	306	// needs to be removed when the fix is in AliROOT and
	307	// the production was done with EvalGraphConst
	308	//TODO TODO TODO
	309	if (fTender->GetRun()<=138154\|\|fTender->GetRun()==138197){
	310	Double_t valDefault = fGainOld->Eval(time);
	311	Double_t valConst = AliTPCcalibDButil::EvalGraphConst(fGainOld, time);
	312	gain = valDefault/valConst;
	313	}
	314	if (fGainNew){
	315	gain *= AliTPCcalibDButil::EvalGraphConst(fGainOld,time)/AliTPCcalibDButil::EvalGraphConst(fGainNew,time);
	316	}
	317	}
e75408ba	318
f19a0c5e	319	//If there is only the default calibration, at least apply correction for pressure
f19a0c5e	320	if (fPcorrection){
e75408ba	321	if (fGRP) {
e75408ba	322	Double_t pressure=fGRP->GetCavernAtmosPressure()->GetValue(time);
ee981ab3	323	// gain=fGainOld->Eval(time)/(7.03814-0.00459798pressure)/49.53fMip;
	324	// gain=fGainOld->Eval(time)/(7.03814-0.00459798pressure)/51.51fMip;
	325	gain=AliTPCcalibDButil::EvalGraphConst(fGainOld,time)/(7.03814-0.00459798pressure)/51.51fMip;
e75408ba	326	}
	327	}
	328
	329	return gain;
	330	}
	331
ee981ab3	332	//_____________________________________________________
	333	void AliTPCTenderSupply::SetParametrisation()
	334	{
	335	//
	336	// Change BB parametrisation for current run
	337	//
	338
	339	//Get CDB Entry with pid response parametrisations
	340	AliCDBEntry *pidCDB=fTender->GetCDBManager()->Get("TPC/Calib/PidResponse",fTender->GetRun());
	341	if (pidCDB){
	342	fArrPidResponseMaster=(TObjArray*)pidCDB->GetObject();
	343	}
	344
	345	//Get the current file to check the reconstruction pass (UGLY, but not stored in ESD... )
0e373bac	346	AliESDInputHandler esdIH = dynamic_cast<AliESDInputHandler> (fTender->GetESDhandler());
	347	if (!esdIH) return;
	348	TTree tree= (TTree)esdIH->GetTree();
	349	TFile file= (TFile) tree->GetCurrentFile();
ee981ab3	350	if (!file) {
	351	AliError("File not found, not changing parametrisation");
	352	return;
	353	}
	354
	355	Int_t run=fTender->GetRun();
	356	AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
	357	Double_t alephParameters[5]={0,0,0,0,0};
	358
	359	//find the period by run number (UGLY, but not stored in ESD... )
	360	TString period;
	361	if (run>=114737&&run<=117223) period="LHC10B";
	362	else if (run>=118503&&run<=121040) period="LHC10C";
	363	else if (run>=122195&&run<=126437) period="LHC10D";
	364	else if (run>=127719&&run<=130850) period="LHC10E";
	365	else if (run>=133004&&run<=135029) period="LHC10F";
	366	else if (run>=135654&&run<=136377) period="LHC10G";
	367	else if (run>=136851&&run<=139517) period="LHC10H";
	368	else if (run>=139699) period="LHC11A";
	369
	370	//find pass from file name (UGLY, but not stored in ESD... )
	371	TString fileName(file->GetName());
	372	Int_t pass=0;
	373	if (fileName.Contains("/pass1")) {
	374	pass=1;
	375	} else if (fileName.Contains("/pass2")) {
	376	pass=2;
	377	}
	378
	379	//beam type
	380	TString beamtype=fTender->GetEvent()->GetBeamType();
	381	if (beamtype.IsNull()\|\|beamtype.Contains("No Beam")) beamtype="p-p";
	382	beamtype.ToUpper();
	383	beamtype.ReplaceAll("-","");
	384
	385	//
	386	// Set default parametrisations for data and MC
	387	//
	388	Bool_t isMC=mgr->GetMCtruthEventHandler();
	389	if (isMC){
	390	//MC data
	391	alephParameters[0] = 2.15898e+00/50.;
	392	alephParameters[1] = 1.75295e+01;
	393	alephParameters[2] = 3.40030e-09;
	394	alephParameters[3] = 1.96178e+00;
	395	alephParameters[4] = 3.91720e+00;
	396
	397	} else {
	398	//data parametrisation
	399
	400
	401	//use defaut data parametrisation in case no other will be selected
	402	alephParameters[0] = 0.0283086/0.97;
	403	alephParameters[1] = 2.63394e+01;
	404	alephParameters[2] = 5.04114e-11;
	405	alephParameters[3] = 2.12543e+00;
	406	alephParameters[4] = 4.88663e+00;
	407
	408	fESDpid->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04);
	409
	410	if (pass==1){
	411
	412	} else if (pass==2){
	413	//find period
414	if (run>=114737&&run<=117223){
415	//LHC10b
416	} else if (run>=118503&&run<=121040) {
417	//LHC10c
418	} else if (run>=122195){
419	//LHC10d +
420	// last run in LHC10d: &&run<126437
421	alephParameters[0] = 1.63246/50.;
422	alephParameters[1] = 2.20028e+01;
423	alephParameters[2] = TMath::Exp(-2.48879e+01);
424	alephParameters[3] = 2.39804e+00;
425	alephParameters[4] = 5.12090e+00;
426
427	//
428	fESDpid->GetTPCResponse().SetSigma(2.30176e-02, 5.60422e+02);
429	}
430	}
431
432	if ( beamtype == "PBPB" ){
433	AliInfo("BETHE-BLOCH parametrization for PbPb !!!!!!!!!!!!!!!!!!!!!!\n");
434
435	alephParameters[0] = 1.25202/50.; //was 1.79571/55.;
436	alephParameters[1] = 2.74992e+01; //was 22.0028;
437	alephParameters[2] = TMath::Exp(-3.31517e+01); //was1.55354e-11;
438	alephParameters[3] = 2.46246; //was 2.39804;
439	alephParameters[4] = 6.78938; //was 5.1209;
440	}
441	}
442
443	fESDpid->GetTPCResponse().SetBetheBlochParameters(
444	alephParameters[0],alephParameters[1],alephParameters[2],
445	alephParameters[3],alephParameters[4]);
446
447	//data type
448	TString datatype="DATA";
449	//in case of mc pass is per default 1
450	if (isMC) {
451	datatype="MC";
452	pass=1;
453	}
454
455	//
456	//set the new parametrisation
457	//
458
459	if (fArrPidResponseMaster){
460	TObject *grAll=0x0;
461	//for MC don't use period information
462	if (isMC) period="[A-Z0-9]*";
463	//pattern for the default entry (valid for all particles)
464	TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),pass,beamtype.Data()));
465
466	//loop over entries and filter them
467	for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
468	TObject *responseFunction=fArrPidResponseMaster->At(iresp);
469	TString responseName=responseFunction->GetName();
470
471	if (!reg.MatchB(responseName)) continue;
472
473	TObjArray *arr=reg.MatchS(responseName);
474	TString particleName=arr->At(1)->GetName();
475	delete arr;
476	if (particleName.IsNull()) continue;
477	if (particleName=="ALL") grAll=responseFunction;
478	else {
479	//find particle id
480	for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
481	TString particle=AliPID::ParticleName(ispec);
482	particle.ToUpper();
483	if ( particle == particleName ){
484	//test if there is already a function set. If yes, cleanup
485	TObject old=const_cast<TObject>(fESDpid->GetTPCResponse().GetResponseFunction((AliPID::EParticleType)ispec));
486	if (old) delete old;
487	fESDpid->GetTPCResponse().SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
488	fESDpid->GetTPCResponse().SetUseDatabase(kTRUE);
489	AliInfo(Form("Adding graph: %d - %s\n",ispec,responseFunction->GetName()));
490	break;
491	}
492	}
493	}
494	}
495
496	//set default response function to all particles which don't have a specific one
497	if (grAll){
498	for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
499	if (!fESDpid->GetTPCResponse().GetResponseFunction((AliPID::EParticleType)ispec)){
500	fESDpid->GetTPCResponse().SetResponseFunction((AliPID::EParticleType)ispec,grAll);
501	AliInfo(Form("Adding graph: %d - %s\n",ispec,grAll->GetName()));
502	}
503	}
504	}
505	}
506	}
507