ANALYSIS/TenderSupplies/AliTPCTenderSupply.cxx

   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>
  24 #include <TObjArray.h>
  25 #include <TObjString.h>
  26 #include <TChain.h>
  27 #include <TFile.h>
  28 #include <TString.h>
  29 #include <TPRegexp.h>
  30 #include <TGraphErrors.h>
  31
  32 #include <AliDCSSensor.h>
  33 #include <AliGRPObject.h>
  34 #include <AliESDpid.h>
  35 #include <AliLog.h>
  36 #include <AliESDEvent.h>
  37 #include <AliExternalTrackParam.h>
  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>
  45 #include <AliCDBRunRange.h>
  46 #include <AliTender.h>
  47 #include <AliTPCcalibDButil.h>
  48 #include <AliPID.h>
  49
  50 #include "AliTPCTenderSupply.h"
  51
  52 ClassImp(AliTPCTenderSupply)
  53
  54 AliTPCTenderSupply::AliTPCTenderSupply() :
  55 AliTenderSupply(),
  56 fESDpid(0x0),
  57 fGainNew(0x0),
  58 fGainOld(0x0),
  59 fGainAttachment(0x0),
  60 fIsMC(kFALSE),
  61 fGainCorrection(kTRUE),
  62 fAttachmentCorrection(kFALSE),
  63 fPcorrection(kFALSE),
  64 fMultiCorrection(kFALSE),
  65 fArrPidResponseMaster(0x0),
  66 fMultiCorrMean(0x0),
  67 fMultiCorrSigma(0x0),
  68 fSpecificStorages(0x0),
  69 fDebugLevel(0),
  70 fMip(50),
  71 fGRP(0x0),
  72 fBeamType("PP"),
  73 fLHCperiod(),
  74 fMCperiod(),
  75 fRecoPass(0)
  76 {
  77   //
  78   // default ctor
  79   //
  80 }
  81
  82 //_____________________________________________________
  83 AliTPCTenderSupply::AliTPCTenderSupply(const char *name, const AliTender *tender) :
  84 AliTenderSupply(name,tender),
  85 fESDpid(0x0),
  86 fGainNew(0x0),
  87 fGainOld(0x0),
  88 fGainAttachment(0x0),
  89 fIsMC(kFALSE),
  90 fGainCorrection(kTRUE),
  91 fAttachmentCorrection(kFALSE),
  92 fPcorrection(kFALSE),
  93 fMultiCorrection(kFALSE),
  94 fArrPidResponseMaster(0x0),
  95 fMultiCorrMean(0x0),
  96 fMultiCorrSigma(0x0),
  97 fSpecificStorages(0x0),
  98 fDebugLevel(0),
  99 fMip(50),
 100 fGRP(0x0),
 101 fBeamType("PP"),
 102 fLHCperiod(),
 103 fMCperiod(),
 104 fRecoPass(0)
 105 {
 106   //
 107   // named ctor
 108   //
 109 }
 110
 111 //_____________________________________________________
 112 void AliTPCTenderSupply::Init()
 113 {
 114   //
 115   // Initialise TPC tender
 116   //
 117
 118   AliLog::SetClassDebugLevel("AliTPCTenderSupply",10);
 119   AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
 120   //
 121   // Setup PID object
 122   //
 123   fIsMC=mgr->GetMCtruthEventHandler();
 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) {
 129     fESDpid=new AliESDpid(fIsMC);
 130     fTender->GetESDhandler()->SetESDpid(fESDpid);
 131   }
 132
 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){
 144     //force no gain and P correction in MC
 145     fGainCorrection=kFALSE;
 146     fPcorrection=kFALSE;
 147     fAttachmentCorrection=kFALSE;
 148   }
 149 }
 150
 151 //_____________________________________________________
 152 void AliTPCTenderSupply::ProcessEvent()
 153 {
 154   //
 155   // Reapply pid information
 156   //
 157
 158   AliESDEvent *event=fTender->GetEvent();
 159   if (!event) return;
 160
 161   //load gain correction if run has changed
 162   if (fTender->RunChanged()){
 163     SetBeamType();
 164     SetRecoInfo();
 165     if ( fBeamType == "PBPB" ) fMultiCorrection=kTRUE;
 166
 167     if (fDebugLevel>0) AliInfo(Form("Run Changed (%d)",fTender->GetRun()));
 168     SetParametrisation();
 169     if (fGainCorrection) SetSplines();
 170   }
 171
 172   //
 173   // get gain correction factor
 174   //
 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());
 180
 181   //
 182   // - correct TPC signals
 183   // - recalculate PID probabilities for TPC
 184   // - correct TPC signal multiplicity dependence
 185
 186   Int_t ntracks=event->GetNumberOfTracks();
 187   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 188     AliESDtrack *track=event->GetTrack(itrack);
 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.5*TMath::Abs(2*inner->GetZ() + (247-83)*inner->GetTgl());
 199     Double_t corrGainTotal=corrFactor*(1 + corrAttachSlope*180.)/(1 + corrAttachSlope*meanDrift)/corrGainMultiplicityPbPb;
 200
 201     // apply gain correction
 202     track->SetTPCsignal(track->GetTPCsignal()*corrGainTotal ,track->GetTPCsignalSigma(), track->GetTPCsignalN());
 203
 204     // recalculate pid probabilities
 205     fESDpid->MakeTPCPID(track);
 206   }
 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;
 267
 268 }
 269
 270 //_____________________________________________________
 271 void AliTPCTenderSupply::SetSplines()
 272 {
 273   //
 274   // Get Gain splines from OCDB
 275   //
 276
 277   AliInfo("Update Gain splines");
 278
 279   //
 280   // Get GPR info for pressure correction
 281   //
 282   fPcorrection=kFALSE;
 283
 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   }
 290   if (fDebugLevel>1) AliInfo(Form("GRP entry used: %s",entryGRP->GetId().ToString().Data()));
 291
 292   fGainNew=0x0;
 293   fGainOld=0x0;
 294   //
 295   //find previous entry from the UserInfo
 296   //
 297 //   TTree *tree=((TChain*)fTender->GetInputData(0))->GetTree();
 298   AliAnalysisManager*mgr = AliAnalysisManager::GetAnalysisManager();
 299   AliAnalysisTaskSE *task = (AliAnalysisTaskSE*)mgr->GetTasks()->First();
 300   if (!task) return;
 301   TTree *tree=((TChain*)task->GetInputData(0))->GetTree();
 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   }
 312
 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   }
 319
 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
 332     if (fDebugLevel>1) AliInfo(Form("Used old Gain entry: %s",entry->GetId().ToString().Data()));
 333
 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   }
 350
 351   //
 352   //new gain correction
 353   //
 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   }
 366   if (!entryNew) {
 367     AliError("No new gain calibration entry found");
 368     return;
 369   }
 370   if (fDebugLevel>1) AliInfo(Form("Used new Gain entry: %s",entryNew->GetId().ToString().Data()));
 371
 372   if (entryNew->GetId().GetLastRun()==AliCDBRunRange::Infinity()){
 373     if (fDebugLevel>0) AliInfo("Use P correction");
 374     fPcorrection=kTRUE;
 375   }
 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
 383   fGainNew        = (AliSplineFit*)arrSplines->At(0);
 384   fGainAttachment = (TGraphErrors*)arrSplines->FindObject("TGRAPHERRORS_MEAN_ATTACHMENT_BEAM_ALL");
 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();
 397
 398   Double_t gain=1;
 399
 400
 401   if (fGainOld){
 402     //TODO TODO TODO
 403     //first correction for the eval const problem
 404     // needs to be removed when the fix is in AliROOT and
 405     // the production was done with EvalGraphConst
 406     // This should be the case from V4-20-Rev11 on
 407     //TODO TODO TODO
 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          ) {
 413       Double_t valDefault = fGainOld->Eval(time);
 414       Double_t valConst   = AliTPCcalibDButil::EvalGraphConst(fGainOld, time);
 415       gain = valDefault/valConst;
 416     }
 417
 418     if (fGainNew){
 419       gain *= AliTPCcalibDButil::EvalGraphConst(fGainOld,time)/AliTPCcalibDButil::EvalGraphConst(fGainNew,time);
 420     }
 421   }
 422
 423   //If there is only the default calibration, at least apply correction for pressure
 424   if (fPcorrection){
 425     if (fGRP) {
 426       Double_t pressure=fGRP->GetCavernAtmosPressure()->GetValue(time);
 427 //       gain=fGainOld->Eval(time)/(7.03814-0.00459798*pressure)/49.53*fMip;
 428 //       gain=fGainOld->Eval(time)/(7.03814-0.00459798*pressure)/51.51*fMip;
 429       gain=AliTPCcalibDButil::EvalGraphConst(fGainOld,time)/(7.03814-0.00459798*pressure)/51.51*fMip;
 430     }
 431   }
 432
 433   return gain;
 434 }
 435
 436 //_____________________________________________________
 437 void AliTPCTenderSupply::SetBeamType()
 438 {
 439   //
 440   // Set the beam type
 441   //
 442
 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   //
 455
 456   //reset information
 457   fRecoPass=0;
 458   fLHCperiod="";
 459   fMCperiod="";
 460
 461   //Get the current file to check the reconstruction pass (UGLY, but not stored in ESD... )
 462   AliESDInputHandler *esdIH = dynamic_cast<AliESDInputHandler*> (fTender->GetESDhandler());
 463   if (!esdIH) return;
 464   TTree *tree= (TTree*)esdIH->GetTree();
 465   TFile *file= (TFile*)tree->GetCurrentFile();
 466   if (!file) {
 467     AliError("Current file not found, cannot set reconstruction information");
 468     return;
 469   }
 470
 471   TString fileName(file->GetName());
 472
 473   Int_t run=fTender->GetRun();
 474
 475
 476   TPRegexp reg(".*(LHC11[a-z]+[0-9]+[a-z]*)/.*");
 477   //find the period by run number (UGLY, but not stored in ESD... )
 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";
 499
 500   //find pass from file name (UGLY, but not stored in ESD... )
 501   if (fileName.Contains("/pass1")) {
 502     fRecoPass=1;
 503   } else if (fileName.Contains("/pass2")) {
 504     fRecoPass=2;
 505   } else if (fileName.Contains("/pass3")) {
 506     fRecoPass=3;
 507   }
 508
 509 }
 510
 511 //_____________________________________________________
 512 void AliTPCTenderSupply::SetParametrisation()
 513 {
 514   //
 515   // Change PbPb multiplicity gain correction factor
 516   //
 517
 518   if (fLHCperiod.IsNull()) {
 519     AliError("No period set, not changing parametrisation");
 520     return;
 521   }
 522
 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   }
 529
 530   if (!fArrPidResponseMaster){
 531     AliError("No valid PidResponse master found in OCDB");
 532     return;
 533   }
 534   //data type
 535   TString datatype="DATA";
 536   TString period=fLHCperiod;
 537   //in case of mc fRecoPass is per default 1
 538   if (fIsMC) {
 539     datatype="MC";
 540     fRecoPass=1;
 541     period=fMCperiod;
 542   }
 543
 544   // Set PbPb correction
 545   fMultiCorrMean=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
 546
 547   if (fMultiCorrMean) AliInfo(Form("Setting multiplicity correction function: %s",fMultiCorrMean->GetName()));
 548
 549 }
 550
 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));
 559 }
 560