PWGCF/EBYE/PIDFluctuation/task/AliEbyEPidRatioTask.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: ALICE Offline.                                                 *
   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 //             AliEbyE Analysis for Particle Ratio Fluctuation             //
  18 //                   Deepika Rathee  | Satyajit Jena                       //
  19 //                   drathee@cern.ch | sjena@cern.ch                       //
  20 //                  Date: Wed Jul  9 18:38:30 CEST 2014                    //
  21 //          New approch to find particle ratio to reduce memory            //
  22 //                             (Test Only)                                 //
  23 //=========================================================================//
  24
  25 #include "TFile.h"
  26 #include "TChain.h"
  27 #include "TTree.h"
  28 #include "TH1F.h"
  29 #include "TH2F.h"
  30 #include "TH3F.h"
  31 #include "TCanvas.h"
  32 #include "TStopwatch.h"
  33 #include "TMath.h"
  34 #include "THashList.h"
  35
  36 #include "AliAnalysisTask.h"
  37 #include "AliAnalysisManager.h"
  38 #include "AliTracker.h"
  39 #include "AliESDEvent.h"
  40 #include "AliESDInputHandler.h"
  41 #include "AliESDpid.h"
  42 #include "AliStack.h"
  43 #include "AliMCEvent.h"
  44 #include "AliMCEventHandler.h"
  45 #include "AliESDtrackCuts.h"
  46 #include "AliKineTrackCuts.h"
  47 #include "AliMCParticle.h"
  48 #include "AliESDVZERO.h"
  49 #include "AliEbyEPidRatioTask.h"
  50 #include "AliGenEventHeader.h"
  51 #include "AliCentrality.h"
  52 #include "AliAODEvent.h"
  53 #include "AliAODInputHandler.h"
  54
  55 using namespace std;
  56 ClassImp(AliEbyEPidRatioTask)
  57 //________________________________________________________________________
  58 AliEbyEPidRatioTask::AliEbyEPidRatioTask(const char *name) :
  59   AliAnalysisTaskSE(name),
  60   fHelper(NULL),
  61   fEffCont(NULL),
  62   fDCA(NULL),
  63   fDist(NULL),
  64   fQA(NULL),
  65
  66   fOutList(NULL),
  67   fOutListEff(NULL),
  68   fOutListCont(NULL),
  69   fOutListDCA(NULL),
  70   fOutListQA(NULL),
  71
  72   fESD(NULL),
  73   fESDHandler(NULL),
  74
  75   fESDTrackCutsBase(NULL),
  76   fESDTrackCuts(NULL),
  77   fESDTrackCutsBkg(NULL),
  78   fESDTrackCutsEff(NULL),
  79
  80   fAOD(NULL),
  81   fAODHandler(NULL),
  82
  83   fIsMC(kFALSE),
  84   fIsRatio(kFALSE),
  85   fIsPtBin(kFALSE),
  86   fIsAOD(kFALSE),
  87   fESDTrackCutMode(0),
  88   fModeEffCreation(0),
  89   fModeDCACreation(0),
  90   fModeDistCreation(0),
  91   fModeQACreation(0),
  92
  93   fMCEvent(NULL),
  94   fMCStack(NULL),
  95
  96   fEtaMax(0.8),
  97   fEtaMaxEff(0.9),
  98   fPtRange(),
  99   fPtRangeEff(),
 100
 101   fAODtrackCutBit(1024) {
 102
 103   AliLog::SetClassDebugLevel("AliEbyEPidRatioTask",10);
 104
 105   fPtRange[0]    = 0.4;
 106   fPtRange[1]    = 0.8;
 107   fPtRangeEff[0] = 0.2;
 108   fPtRangeEff[1] = 1.6;
 109
 110   DefineOutput(1, TList::Class());
 111   DefineOutput(2, TList::Class());
 112   DefineOutput(3, TList::Class());
 113   DefineOutput(4, TList::Class());
 114   DefineOutput(5, TList::Class());
 115 }
 116
 117 //________________________________________________________________________
 118 AliEbyEPidRatioTask::~AliEbyEPidRatioTask() {
 119   // Destructor
 120
 121   if (fESDTrackCutsBase) delete fESDTrackCutsBase;
 122   if (fESDTrackCuts)     delete fESDTrackCuts;
 123   if (fESDTrackCutsBkg)  delete fESDTrackCutsBkg;
 124   if (fESDTrackCutsEff)  delete fESDTrackCutsEff;
 125
 126   if (fEffCont)          delete fEffCont;
 127   if (fDCA)              delete fDCA;
 128   if (fDist)             delete fDist;
 129   if (fQA)               delete fQA;
 130   if (fHelper)           delete fHelper;
 131 }
 132
 133
 134
 135 //________________________________________________________________________
 136 void AliEbyEPidRatioTask::UserCreateOutputObjects() {
 137   Bool_t oldStatus = TH1::AddDirectoryStatus();
 138   TH1::AddDirectory(kFALSE);
 139
 140   fOutList = new TList;
 141   fOutList->SetName(GetName()) ;
 142   fOutList->SetOwner(kTRUE);
 143
 144   fOutListEff = new TList;
 145   fOutListEff->SetName(Form("%s_eff",GetName()));
 146   fOutListEff->SetOwner(kTRUE) ;
 147
 148   fOutListCont = new TList;
 149   fOutListCont->SetName(Form("%s_cont",GetName()));
 150   fOutListCont->SetOwner(kTRUE) ;
 151
 152   fOutListDCA = new TList;
 153   fOutListDCA->SetName(Form("%s_dca",GetName()));
 154   fOutListDCA->SetOwner(kTRUE) ;
 155
 156   fOutListQA = new TList;
 157   fOutListQA->SetName(Form("%s_qa",GetName()));
 158   fOutListQA->SetOwner(kTRUE) ;
 159
 160   Initialize();
 161
 162   fOutList->Add(new TList);
 163   TList *list = static_cast<TList*>(fOutList->Last());
 164   list->SetName(Form("fStat"));
 165   list->SetOwner(kTRUE);
 166
 167   list->Add(fHelper->GetHEventStat0());
 168   list->Add(fHelper->GetHEventStat1());
 169   list->Add(fHelper->GetHTriggerStat());
 170   list->Add(fHelper->GetHCentralityPercentile());
 171   list->Add(fHelper->GetHCentralityPercentileAll());
 172
 173
 174
 175
 176   if ((fIsAOD||fIsMC) && fModeEffCreation == 1) {
 177     fOutListEff->Add(fEffCont->GetHnEffMc());
 178     fOutListEff->Add(fEffCont->GetHnEffRec());
 179
 180     fOutListCont->Add(fEffCont->GetHnContMc());
 181     fOutListCont->Add(fEffCont->GetHnContRec());
 182   }
 183
 184   if (fModeDCACreation == 1)
 185     fOutListDCA->Add(fDCA->GetHnDCA());
 186
 187   if (fModeQACreation == 1) {
 188     fOutListQA->Add(fQA->GetHnQAPid());
 189     fOutListQA->Add(fQA->GetHnQADca());
 190   }
 191
 192   TH1::AddDirectory(oldStatus);
 193
 194   PostData(1,fOutList);
 195   PostData(2,fOutListEff);
 196   PostData(3,fOutListCont);
 197   PostData(4,fOutListDCA);
 198   PostData(5,fOutListQA);
 199
 200   return;
 201 }
 202
 203 //________________________________________________________________________
 204 void AliEbyEPidRatioTask::UserExec(Option_t *) {
 205
 206   if (SetupEvent() < 0) {
 207     PostData(1,fOutList);
 208     PostData(2,fOutListEff);
 209     PostData(3,fOutListCont);
 210     PostData(4,fOutListDCA);
 211     PostData(5,fOutListQA);
 212     return;
 213   }
 214
 215
 216   if ((fIsMC||fIsAOD) && fModeEffCreation == 1)
 217     fEffCont->Process();
 218
 219   if (fModeDCACreation == 1)
 220     fDCA->Process();
 221
 222   if (fModeDistCreation == 1)
 223     fDist->Process();
 224
 225   if (fModeQACreation == 1)
 226     fQA->Process();
 227
 228   PostData(1,fOutList);
 229   PostData(2,fOutListEff);
 230   PostData(3,fOutListCont);
 231   PostData(4,fOutListDCA);
 232   PostData(5,fOutListQA);
 233
 234   return;
 235 }
 236
 237 //________________________________________________________________________
 238 void AliEbyEPidRatioTask::Terminate(Option_t *){
 239   // Terminate
 240 }
 241
 242 Int_t AliEbyEPidRatioTask::Initialize() {
 243   // Initialize event
 244
 245   // ------------------------------------------------------------------
 246   // -- ESD TrackCuts
 247   // ------------------------------------------------------------------
 248   TString sModeName("");
 249
 250   // -- Create ESD track cuts
 251   // --------------------------
 252   fESDTrackCutsBase = new AliESDtrackCuts;
 253
 254   if (fESDTrackCutMode == 0) {
 255     fESDTrackCutsBase->SetMinNCrossedRowsTPC(70);                                             // TPC
 256     fESDTrackCutsBase->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);                    // TPC
 257   }
 258   else if (fESDTrackCutMode == 1) {
 259     fESDTrackCutsBase->SetMinNClustersTPC(70);                                                // TPC  2010
 260   }
 261
 262   fESDTrackCutsBase->SetMaxChi2PerClusterTPC(4);                                              // TPC  2010
 263   fESDTrackCutsBase->SetAcceptKinkDaughters(kFALSE);                                          // TPC  2010
 264   fESDTrackCutsBase->SetRequireTPCRefit(kTRUE);                                               // TPC  2010
 265
 266   if (fESDTrackCutMode == 0) {
 267     fESDTrackCutsBase->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kOff); // ITS
 268     fESDTrackCutsBase->SetClusterRequirementITS(AliESDtrackCuts::kSDD,AliESDtrackCuts::kOff); // ITS
 269     fESDTrackCutsBase->SetClusterRequirementITS(AliESDtrackCuts::kSSD,AliESDtrackCuts::kOff); // ITS
 270   }
 271   else if (fESDTrackCutMode == 1) {
 272     fESDTrackCutsBase->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny); // ITS 2010
 273   //  fESDTrackCutsBase->SetMinNClustersITS(4);
 274   }
 275
 276   fESDTrackCutsBase->SetRequireITSRefit(kTRUE);                                               // ITS 2010
 277   fESDTrackCutsBase->SetMaxChi2PerClusterITS(36);                                             // ITS 2010
 278
 279   fESDTrackCutsBase->SetDCAToVertex2D(kFALSE);                                                // VertexConstrained  2010
 280   fESDTrackCutsBase->SetRequireSigmaToVertex(kFALSE);                                         // VertexConstrained  2010
 281   fESDTrackCutsBase->SetMaxDCAToVertexZ(2);                                                   // VertexConstrained  2010
 282
 283   fESDTrackCutsBase->SetEtaRange(-1.*fEtaMax, fEtaMax);                                       // Acceptance
 284   fESDTrackCutsBase->SetPtRange(fPtRange[0],fPtRange[1]);                                     // Acceptance
 285
 286   // -- Mode : standard cuts
 287   if (fESDTrackCutMode == 0)
 288     sModeName = "Std";
 289   // -- Mode : for comparison to LF
 290   else if (fESDTrackCutMode == 1)
 291     sModeName = "LF";
 292   // -- Mode : Default
 293   else
 294     sModeName = "Base";
 295
 296   fESDTrackCutsBase->SetName(Form("NetParticleCuts2010_%s",sModeName.Data()));
 297
 298   // -- Create ESD track cuts -> Base + DCA
 299   // ------------------------------
 300   fESDTrackCuts = static_cast<AliESDtrackCuts*>(fESDTrackCutsBase->Clone());
 301   fESDTrackCuts->SetName(Form("NetParticleCuts2010_%s",sModeName.Data()));
 302   if (fESDTrackCutMode == 0)
 303     fESDTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");                       // 2010 VertexConstrained  ->  7*(0.0026+0.0050/pt^1.01)
 304     //    fESDTrackCuts->SetMaxDCAToVertexXY(0.3);
 305   else if (fESDTrackCutMode == 1)
 306     fESDTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");                       // 2010 VertexConstrained  ->  7*(0.0026+0.0050/pt^1.01)
 307
 308   //  fESDTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);                                    // golden cut off
 309
 310   // -- Create ESD BKG track cuts -> Base + Acceptance(Eff)
 311   // ------------------------------
 312   fESDTrackCutsBkg = static_cast<AliESDtrackCuts*>(fESDTrackCutsBase->Clone());
 313   fESDTrackCutsBkg->SetName(Form("NetParticleCuts2010_%s_Bkg",sModeName.Data()));
 314   fESDTrackCutsBkg->SetPtRange(fPtRangeEff[0],fPtRangeEff[1]);                              // Acceptance
 315   fESDTrackCutsBkg->SetEtaRange(-1.*fEtaMaxEff, fEtaMaxEff);                                // Acceptance
 316
 317   // -- Create ESD Eff track cuts -> Base + DCA + Acceptance(Eff)
 318   // ------------------------------
 319   fESDTrackCutsEff = static_cast<AliESDtrackCuts*>(fESDTrackCuts->Clone());
 320   fESDTrackCutsEff->SetName(Form("NetParticleCuts2010_%s_Eff",sModeName.Data()));
 321   fESDTrackCutsEff->SetPtRange(fPtRangeEff[0],fPtRangeEff[1]);                              // Acceptance
 322   fESDTrackCutsEff->SetEtaRange(-1.*fEtaMaxEff, fEtaMaxEff);                                // Acceptance
 323
 324   // ------------------------------------------------------------------
 325   // -- Initialize Helper
 326   // ------------------------------------------------------------------
 327
 328
 329   if (fHelper->Initialize(fESDTrackCutsEff, fIsMC,fIsRatio,fIsPtBin, fAODtrackCutBit, fModeDistCreation))
 330     return -1;
 331
 332   // fHelper->SetIsRatio(fIsRatio);
 333   // fHelper->SetIsPtBin(fIsPtBin);
 334
 335   // ------------------------------------------------------------------
 336   // -- Create / Initialize Efficiency/Contamination
 337   // ------------------------------------------------------------------
 338   if ((fIsMC||fIsAOD) && fModeEffCreation == 1) {
 339     fEffCont = new AliEbyEPidRatioEffCont;
 340     fEffCont->Initialize(fHelper, fESDTrackCuts);
 341   }
 342
 343   // ------------------------------------------------------------------
 344   // -- Create / Initialize DCA Determination
 345   // ------------------------------------------------------------------
 346   if (fModeDCACreation == 1) {
 347     fDCA = new AliEbyEPidRatioDCA;
 348     fDCA->SetESDTrackCutsBkg(fESDTrackCutsBkg);
 349     fDCA->Initialize(fHelper, fESDTrackCuts);
 350   }
 351
 352   // ------------------------------------------------------------------
 353   // -- Create / Initialize Phy Determination
 354   // ------------------------------------------------------------------
 355   if (fModeDistCreation == 1) {
 356     fDist = new AliEbyEPidRatioPhy;
 357     fDist->SetOutList(fOutList);
 358     fDist->Initialize(fHelper, fESDTrackCuts);
 359   }
 360
 361   // ------------------------------------------------------------------
 362   // -- Create / Initialize QA Determination
 363   // ------------------------------------------------------------------
 364   if (fModeQACreation == 1) {
 365     fQA = new AliEbyEPidRatioQA();
 366     fQA->Initialize(fHelper, fESDTrackCuts);
 367   }
 368
 369   // ------------------------------------------------------------------
 370   // -- Reset Event
 371   // ------------------------------------------------------------------
 372   ResetEvent();
 373
 374   return 0;
 375 }
 376
 377 //________________________________________________________________________
 378 Int_t AliEbyEPidRatioTask::SetupEvent() {
 379
 380   ResetEvent();
 381
 382   // -- ESD Event
 383   // ------------------------------------------------------------------
 384   if (!fIsAOD && SetupESDEvent() < 0) {
 385     AliError("Setup ESD Event failed");
 386     return -1;
 387   }
 388
 389   // -- AOD Event
 390   // ------------------------------------------------------------------
 391   if (fIsAOD && SetupAODEvent() < 0) {
 392     AliError("Setup AOD Event failed");
 393     return -1;
 394   }
 395
 396   // -- Setup MC Event
 397   // ------------------------------------------------------------------
 398   if (fIsMC && SetupMCEvent() < 0) {
 399     AliError("Setup MC Event failed");
 400     return -1;
 401   }
 402
 403   // -- Setup Event for Helper / EffCont  / DCA / Dist / QA classes
 404   // ------------------------------------------------------------------
 405   fHelper->SetupEvent(fESDHandler, fAODHandler, fMCEvent);
 406
 407
 408   if (fModeEffCreation && (fIsMC || fIsAOD) )
 409     fEffCont->SetupEvent();
 410
 411   if (fModeDCACreation == 1)
 412     fDCA->SetupEvent();
 413
 414   if (fModeDistCreation == 1)
 415     fDist->SetupEvent();
 416
 417   if (fModeQACreation == 1)
 418     fQA->SetupEvent();
 419
 420
 421   // -- Evaluate Event cuts
 422   // ------------------------------------------------------------------
 423   return fHelper->IsEventRejected() ? -2 : 0;
 424 }
 425
 426 //________________________________________________________________________
 427 Int_t AliEbyEPidRatioTask::SetupESDEvent() {
 428   // -- Setup ESD Event
 429   // > return 0 for success
 430   // > return -1 for failed setup
 431
 432
 433
 434
 435   fESDHandler= dynamic_cast<AliESDInputHandler*>
 436     (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 437   if (!fESDHandler) {
 438     AliError("Could not get ESD input handler");
 439     return -1;
 440   }
 441
 442   fESD = fESDHandler->GetEvent();
 443   if (!fESD) {
 444     AliError("Could not get ESD event");
 445     return -1;
 446   }
 447
 448   // -- Check PID response
 449   // ------------------------------------------------------------------
 450   if (!fESDHandler->GetPIDResponse()) {
 451     AliError("Could not get PID response");
 452     return -1;
 453   }
 454
 455   // -- Check Vertex
 456   // ------------------------------------------------------------------
 457   if (!fESD->GetPrimaryVertexTracks()) {
 458     AliError("Could not get vertex from tracks");
 459     return -1;
 460   }
 461
 462   // -- Check Centrality
 463   // ------------------------------------------------------------------
 464   if (!fESD->GetCentrality()) {
 465     AliError("Could not get centrality");
 466     return -1;
 467   }
 468
 469   return 0;
 470 }
 471
 472 //________________________________________________________________________
 473 Int_t AliEbyEPidRatioTask::SetupAODEvent() {
 474   // -- Setup AOD Event
 475   // > return 0 for success
 476   // > return -1 for failed setup
 477
 478   fAODHandler= dynamic_cast<AliAODInputHandler*>
 479     (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 480   if (!fAODHandler) {
 481     AliError("Could not get AOD input handler");
 482     return -1;
 483   }
 484
 485   fAOD = fAODHandler->GetEvent();
 486   if (!fAOD) {
 487     AliError("Could not get AOD event");
 488     return -1;
 489   }
 490
 491   // -- Check PID response
 492   // ------------------------------------------------------------------
 493   if (!fAODHandler->GetPIDResponse()) {
 494     AliError("Could not get PID response");
 495     return -1;
 496   }
 497
 498   // -- Check Vertex
 499   // ------------------------------------------------------------------
 500   if (!fAOD->GetPrimaryVertex()) {
 501     AliError("Could not get primary vertex");
 502     return -1;
 503   }
 504
 505   // -- Check Centrality
 506   // ------------------------------------------------------------------
 507   if (!fAOD->GetHeader()->GetCentralityP()) {
 508     AliError("Could not get centrality");
 509     return -1;
 510   }
 511
 512   return 0;
 513 }
 514
 515 //________________________________________________________________________
 516 Int_t AliEbyEPidRatioTask::SetupMCEvent() {
 517   AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*>
 518     (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 519
 520   if (!mcH) {
 521     AliError("MC event handler not available");
 522     return -1;
 523   }
 524
 525   fMCEvent = mcH->MCEvent();
 526   if (!fMCEvent) {
 527     AliError("MC event not available");
 528     return -1;
 529   }
 530
 531   // -- Get MC header
 532   // ------------------------------------------------------------------
 533   AliHeader* header = fMCEvent->Header();
 534   if (!header) {
 535     AliError("MC header not available");
 536     return -1;
 537   }
 538
 539   // -- Check Stack
 540   // ------------------------------------------------------------------
 541   fMCStack = fMCEvent->Stack();
 542   if (!fMCStack) {
 543     AliError("MC stack not available");
 544     return -1;
 545   }
 546
 547   // -- Check GenHeader
 548   // ------------------------------------------------------------------
 549   if (!header->GenEventHeader()) {
 550     AliError("Could not retrieve genHeader from header");
 551     return -1;
 552   }
 553
 554   // -- Check primary vertex
 555   // ------------------------------------------------------------------
 556   if (!fMCEvent->GetPrimaryVertex()){
 557     AliError("Could not get MC vertex");
 558     return -1;
 559   }
 560
 561   return 0;
 562 }
 563
 564 //________________________________________________________________________
 565 void AliEbyEPidRatioTask::ResetEvent() {
 566   // -- Reset event
 567
 568   // -- Reset ESD Event
 569   fESD       = NULL;
 570
 571   // -- Reset AOD Event
 572   fAOD       = NULL;
 573
 574   // -- Reset MC Event
 575   if (fIsMC)
 576     fMCEvent = NULL;
 577
 578   // -- Reset Dist Creation
 579   if (fModeDistCreation == 1)
 580     fDist->ResetEvent();
 581
 582   return;
 583 }
 584
 585