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