PWGLF/totEt/AliAnalysisTaskHadEt.cxx

   1 //_________________________________________________________________________
   2 //  Utility Class for transverse energy studies; charged hadrons
   3 //  Task for analysis
   4 //  - reconstruction and MC output
   5 // implementation file
   6 //
   7 //Created by Christine Nattrass, Rebecca Scott, Irakli Martashvili
   8 //University of Tennessee at Knoxville
   9 //_________________________________________________________________________
  10 //Necessary to read config macros
  11 #include <TROOT.h>
  12 #include <TSystem.h>
  13 #include <TInterpreter.h>
  14
  15 #include "TChain.h"
  16 #include "TList.h"
  17 #include "TH2F.h"
  18
  19 #include "AliESDEvent.h"
  20 #include "AliMCEvent.h"
  21 #include "AliESDtrackCuts.h"
  22
  23 #include "AliAnalysisTaskHadEt.h"
  24 #include "AliAnalysisHadEtReconstructed.h"
  25 #include "AliAnalysisHadEtMonteCarlo.h"
  26 #include "AliPWG0Helper.h"
  27 #include "AliTriggerAnalysis.h"
  28 #include "AliAnalysisManager.h"
  29 #include "AliPIDResponse.h"
  30 #include "AliTPCPIDResponse.h"
  31 #include "AliInputEventHandler.h"
  32
  33 #include <iostream>
  34 #include "AliLog.h"
  35
  36 using namespace std;
  37
  38 ClassImp(AliAnalysisTaskHadEt)
  39
  40
  41
  42 //________________________________________________________________________
  43   AliAnalysisTaskHadEt::AliAnalysisTaskHadEt(const char *name, Bool_t isMc, TString recoConfigFile, TString mcConfigFile) :
  44         AliAnalysisTaskTransverseEnergy(name, isMc)
  45         ,fPIDResponse(0)
  46         ,fRecAnalysis(0)
  47         ,fMCAnalysis(0)
  48         ,fIsSim(isMc)
  49         ,kIsOfflineV0AND(0)
  50         ,kIsOfflineMB(0)
  51 {
  52     // Constructor
  53   //input hander
  54   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
  55   if (!man) {
  56     AliFatal("Analysis manager needed");
  57     return;
  58   }
  59
  60   AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
  61   if (!inputHandler) {
  62     AliFatal("Input handler needed");
  63     return;
  64   }
  65   inputHandler->SetNeedField();
  66
  67   //pid response object
  68   fPIDResponse=inputHandler->GetPIDResponse();
  69   if (!fPIDResponse) AliError("PIDResponse object was not created");
  70   else{cout<<"PIDResponse was created!"<<endl;}
  71
  72
  73   fMCConfigFile = mcConfigFile;
  74   fRecoConfigFile = recoConfigFile;
  75
  76   if(fMCAnalysis) delete fMCAnalysis;
  77   if(fRecAnalysis) delete fRecAnalysis;
  78
  79   if (fRecoConfigFile.Length()) {
  80     cout<<"Rereading AliAnalysisHadEtReconstructed configuration file "<<fRecoConfigFile<<endl;
  81     gROOT->LoadMacro(fRecoConfigFile);
  82     fRecAnalysis = (AliAnalysisHadEtReconstructed *) gInterpreter->ProcessLine("ConfigHadEtReconstructed()");
  83   }
  84
  85   if (fMCConfigFile.Length()) {
  86     cout<<"Rereading AliAnalysisHadEtMonteCarlo configuration file "<<fMCConfigFile<<endl;
  87     gROOT->LoadMacro(fMCConfigFile);
  88     fMCAnalysis = (AliAnalysisHadEtMonteCarlo *) gInterpreter->ProcessLine("ConfigHadEtMonteCarlo()");
  89     fMCAnalysis->SetHadEtReconstructed(fRecAnalysis);
  90   }
  91
  92     // Define input and output slots here
  93     // Input slot #0 works with a TChain
  94     DefineInput(0, TChain::Class());
  95     // Output slot #1 writes into a TH1 container
  96
  97     DefineOutput(1, TList::Class());
  98 }
  99 AliAnalysisTaskHadEt::~AliAnalysisTaskHadEt(){//Destructor
 100   delete fRecAnalysis;
 101   delete fMCAnalysis;
 102   delete fPIDResponse;
 103 }
 104
 105
 106 //________________________________________________________________________
 107 void AliAnalysisTaskHadEt::UserCreateOutputObjects()
 108 {
 109     // Create histograms
 110
 111     // Called once
 112
 113
 114   //input hander
 115   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 116   if (!man) {
 117     AliFatal("Analysis manager needed");
 118     return;
 119   }
 120
 121   AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
 122   if (!inputHandler) {
 123     AliFatal("Input handler needed");
 124     return;
 125   }
 126
 127   //pid response object
 128   fPIDResponse=inputHandler->GetPIDResponse();
 129   if (!fPIDResponse) AliError("PIDResponse object was not created");
 130   else{cout<<"PIDResponse was created!"<<endl;}
 131
 132
 133   fOutputList = new TList;
 134   fOutputList->SetOwner();
 135   fMCAnalysis->SetHistoList(fOutputList);
 136   fRecAnalysis->SetHistoList(fOutputList);
 137   if(fIsSim) fMCAnalysis->CreateHistograms();
 138   fRecAnalysis->CreateHistograms();
 139
 140
 141   if(fRecAnalysis->DataSet() != fMCAnalysis->DataSet()){
 142     cout<<"Warning: Reconstruction data set and Monte Carlo data set are not the same!  Setting data set to "<<fRecAnalysis->DataSet()<<endl;
 143   }
 144
 145   Bool_t selectPrimaries=kTRUE;
 146   if(fEsdtrackCutsITSTPC) delete fEsdtrackCutsITSTPC;
 147   if(fEsdtrackCutsITS) delete fEsdtrackCutsITS;
 148   if(fEsdtrackCutsTPC) delete fEsdtrackCutsTPC;
 149   //We do not use these because we are using the 2010 900 GeV data
 150 //   if(fRecAnalysis->DataSet()==2009){
 151 //     cout<<"Setting track cuts for the 2009 p+p collisions at 900 GeV"<<endl;
 152 //     fEsdtrackCutsITSTPC = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selectPrimaries);
 153 //     fEsdtrackCutsITSTPC->SetName("fEsdTrackCuts");
 154 //     fEsdtrackCutsTPC = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 155 //     fEsdtrackCutsTPC->SetName("fEsdTrackCutsTPCOnly");
 156 //     //ITS stand alone cuts - similar to 2009 cuts but with only ITS hits required
 157 //     fEsdtrackCutsITS =  AliESDtrackCuts::GetStandardITSPureSATrackCuts2009(kTRUE,kFALSE);//we do want primaries but we do not want to require PID info
 158 //     fEsdtrackCutsITS->SetName("fEsdTrackCutsITS");
 159 //   }
 160   if(fRecAnalysis->DataSet()==2010 || fRecAnalysis->DataSet()==20111||fRecAnalysis->DataSet()==2009 || fRecAnalysis->DataSet()==2012 || fRecAnalysis->DataSet()==2013){
 161     // AliAnalysisTaskSE::      SelectCollisionCandidates(AliVEvent::kINT7 ) ;
 162     if(fRecAnalysis->DataSet()==2010)cout<<"Setting track cuts for the 2010 p+p collisions at 7 TeV"<<endl;
 163     else{
 164       if(fRecAnalysis->DataSet()==2012)cout<<"Setting track cuts for the 2012 p+p collisions at 8 TeV"<<endl;
 165       else{
 166         if(fRecAnalysis->DataSet()==2013)cout<<"Setting track cuts for the 2013 p+Pb collisions at 5 TeV"<<endl;
 167         else{
 168           if(fRecAnalysis->DataSet()==2009){cout<<"Setting track cuts for the 2010 p+p collisions at 900 GeV"<<endl;}
 169           else{cout<<"Setting track cuts for the 2011 p+p collisions at 2.76 TeV"<<endl;}
 170         }
 171       }
 172     }
 173     //cout<<"Warning:  Have not set 2010 track cuts yet!!"<<endl;
 174     fEsdtrackCutsITSTPC = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
 175     fEsdtrackCutsITSTPC->SetName("fEsdTrackCuts");
 176     fEsdtrackCutsTPC = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 177     fEsdtrackCutsTPC->SetName("fEsdTrackCutsTPCOnly");
 178     //ITS stand alone cuts - similar to 2009 cuts but with only ITS hits required
 179     fEsdtrackCutsITS =  AliESDtrackCuts::GetStandardITSSATrackCuts2010(kTRUE,kFALSE);//we do want primaries but we do not want to require PID info
 180     fEsdtrackCutsITS->SetName("fEsdTrackCutsITS");
 181   }
 182   if(fRecAnalysis->DataSet()==20100 || fRecAnalysis->DataSet()==2011){
 183     cout<<"Setting track cuts for the 2010 Pb+Pb collisions at 2.76 TeV"<<endl;
 184     //cout<<"Warning:  Have not set 2010 track cuts yet!!"<<endl;
 185     fEsdtrackCutsITSTPC = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
 186     fEsdtrackCutsITSTPC->SetName("fEsdTrackCuts");
 187     fEsdtrackCutsTPC = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 188     fEsdtrackCutsTPC->SetName("fEsdTrackCutsTPCOnly");
 189     //ITS stand alone cuts - similar to 2009 cuts but with only ITS hits required
 190     fEsdtrackCutsITS =  AliESDtrackCuts::GetStandardITSSATrackCutsPbPb2010(kTRUE,kFALSE);//we do want primaries but we do not want to require PID info
 191     // fEsdtrackCutsITS =  AliESDtrackCuts::GetStandardITSPureSATrackCuts2010(kTRUE,kFALSE);//we do want primaries but we do not want to require PID info
 192    fEsdtrackCutsITS->SetName("fEsdTrackCutsITS");
 193   }
 194
 195   fOutputList->Add(fEsdtrackCutsITSTPC);
 196   fOutputList->Add(fEsdtrackCutsTPC);
 197   fOutputList->Add(fEsdtrackCutsITS);
 198   if(fEsdtrackCutsITSTPC && fEsdtrackCutsTPC){
 199     fRecAnalysis->SetITSTrackCuts( GetITSTrackCuts());
 200     fMCAnalysis->SetITSTrackCuts( GetITSTrackCuts());
 201     fRecAnalysis->SetTPCITSTrackCuts( GetTPCITSTrackCuts());
 202     fMCAnalysis->SetTPCITSTrackCuts( GetTPCITSTrackCuts());
 203     fRecAnalysis->SetTPCOnlyTrackCuts( GetTPCOnlyTrackCuts());
 204     fMCAnalysis->SetTPCOnlyTrackCuts( GetTPCOnlyTrackCuts());
 205     //add ITS stuff!
 206   }
 207   else{
 208     Printf("Error: no track cuts!");
 209   }
 210
 211
 212
 213  PostData(1, fOutputList);
 214 }
 215
 216 //________________________________________________________________________
 217 void AliAnalysisTaskHadEt::UserExec(Option_t *)
 218 { // execute method
 219   fESDEvent = dynamic_cast<AliESDEvent*>(InputEvent());
 220 if (!fESDEvent) {
 221   Printf("ERROR: Could not retrieve event");
 222   return;
 223  }
 224 //cout<<"AliAnalysisTaskHadEt 165"<<endl;
 225
 226 //Int_t res = CheckPhysicsSelection(fESDEvent->GetRunNumber()); // Check if the physics selection is valid for this run
 227
 228 //AliCentrality *cent = GetCentralityObject();
 229
 230 //if(res == 0 && cent){
 231 //if(cent){
 232  AliTriggerAnalysis *fTriggerAnalysis = new AliTriggerAnalysis();
 233
 234   kIsOfflineV0AND = fTriggerAnalysis->IsOfflineTriggerFired(fESDEvent, AliTriggerAnalysis::kV0AND);
 235   kIsOfflineMB = fTriggerAnalysis->IsOfflineTriggerFired(fESDEvent, AliTriggerAnalysis::kMB1);
 236   fRecAnalysis->SetIsOfflineV0AND(kIsOfflineV0AND);
 237   fMCAnalysis->SetIsOfflineV0AND(kIsOfflineV0AND);
 238   fMCAnalysis->SetIsOfflineMB(kIsOfflineMB);
 239
 240   Int_t eventtype =     AliPWG0Helper::kInvalidProcess;
 241   if(fIsSim &&( fRecAnalysis->DataSet()!=20100 || fRecAnalysis->DataSet()!=2011) ) eventtype = (Int_t) AliPWG0Helper::GetEventProcessType(MCEvent()->Header());
 242   //only do the analysis if it meets the offline trigger cut
 243   if(kIsOfflineV0AND) fRecAnalysis->AnalyseEvent(fESDEvent,eventtype);
 244   //else{cout<<"Not analyzing this event!  Does not meet trigger condition!"<<endl;}
 245   if(fIsSim){
 246     AliMCEvent* mcEvent = MCEvent();
 247     if(!mcEvent){
 248       AliFatal("ERROR: MC Event does not exist");
 249       return;
 250     }
 251     if (fESDEvent){
 252       ((AliAnalysisHadEtMonteCarlo*)fMCAnalysis)->AnalyseEvent((AliVEvent*)mcEvent,(AliVEvent*)fESDEvent);
 253     if( fRecAnalysis->DataSet()==20100 || AliPWG0Helper::GetEventProcessType(mcEvent->Header()) == AliPWG0Helper::kND || fRecAnalysis->DataSet()==2013|| fRecAnalysis->DataSet()==2011){//either non-diffractive or Pb+Pb or p+Pb
 254       if(fMCAnalysis->Full()){
 255         fMCAnalysis->FillSimTotEtMinusRecoTotEtFullAcceptanceTPC( fRecAnalysis->GetCorrectedTotEtFullAcceptanceTPC() );
 256         fMCAnalysis->FillSimTotEtMinusRecoTotEtFullAcceptanceITS( fRecAnalysis->GetCorrectedTotEtFullAcceptanceITS() );
 257         fMCAnalysis->FillSimTotEtMinusRecoTotEtFullAcceptanceTPCNoPID( fRecAnalysis->GetCorrectedTotEtFullAcceptanceTPCNoPID() );
 258         fMCAnalysis->FillSimTotEtMinusRecoTotEtFullAcceptanceITSNoPID( fRecAnalysis->GetCorrectedTotEtFullAcceptanceITSNoPID() );
 259         fMCAnalysis->FillSimHadEtMinusRecoHadEtFullAcceptanceTPC( fRecAnalysis->GetCorrectedHadEtFullAcceptanceTPC() );
 260         fMCAnalysis->FillSimHadEtMinusRecoHadEtFullAcceptanceITS( fRecAnalysis->GetCorrectedHadEtFullAcceptanceITS() );
 261         fMCAnalysis->FillSimHadEtMinusRecoHadEtFullAcceptanceTPCNoPID( fRecAnalysis->GetCorrectedHadEtFullAcceptanceTPCNoPID() );
 262         fMCAnalysis->FillSimHadEtMinusRecoHadEtFullAcceptanceITSNoPID( fRecAnalysis->GetCorrectedHadEtFullAcceptanceITSNoPID() );
 263
 264         fMCAnalysis->FillSimTotEtVsRecoTotEtFullAcceptanceTPC( fRecAnalysis->GetCorrectedTotEtFullAcceptanceTPC() );
 265         fMCAnalysis->FillSimTotEtVsRecoTotEtFullAcceptanceITS( fRecAnalysis->GetCorrectedTotEtFullAcceptanceITS() );
 266         fMCAnalysis->FillSimTotEtVsRecoTotEtFullAcceptanceTPCNoPID( fRecAnalysis->GetCorrectedTotEtFullAcceptanceTPCNoPID() );
 267         fMCAnalysis->FillSimTotEtVsRecoTotEtFullAcceptanceITSNoPID( fRecAnalysis->GetCorrectedTotEtFullAcceptanceITSNoPID() );
 268         fMCAnalysis->FillSimHadEtVsRecoHadEtFullAcceptanceTPC( fRecAnalysis->GetCorrectedHadEtFullAcceptanceTPC() );
 269         fMCAnalysis->FillSimHadEtVsRecoHadEtFullAcceptanceITS( fRecAnalysis->GetCorrectedHadEtFullAcceptanceITS() );
 270         fMCAnalysis->FillSimHadEtVsRecoHadEtFullAcceptanceTPCNoPID( fRecAnalysis->GetCorrectedHadEtFullAcceptanceTPCNoPID() );
 271         fMCAnalysis->FillSimHadEtVsRecoHadEtFullAcceptanceITSNoPID( fRecAnalysis->GetCorrectedHadEtFullAcceptanceITSNoPID() );
 272
 273         fMCAnalysis->FillSimPiKPEtVsRecoPiKPEtFullAcceptanceTPC( fRecAnalysis->GetCorrectedPiKPEtFullAcceptanceTPC() );
 274         fMCAnalysis->FillSimPiKPEtVsRecoPiKPEtFullAcceptanceITS( fRecAnalysis->GetCorrectedPiKPEtFullAcceptanceITS() );
 275         fMCAnalysis->FillSimPiKPEtVsRecoPiKPEtFullAcceptanceTPCNoPID( fRecAnalysis->GetCorrectedPiKPEtFullAcceptanceTPCNoPID() );
 276         fMCAnalysis->FillSimPiKPEtVsRecoPiKPEtFullAcceptanceITSNoPID( fRecAnalysis->GetCorrectedPiKPEtFullAcceptanceITSNoPID() );//Had
 277
 278         fMCAnalysis->FillSimRawEtVsRecoRawEtFullAcceptanceTPC( fRecAnalysis->GetRawEtFullAcceptanceTPC() );
 279         fMCAnalysis->FillSimRawEtVsRecoRawEtFullAcceptanceITS( fRecAnalysis->GetRawEtFullAcceptanceITS() );
 280         fMCAnalysis->FillSimRawEtVsRecoRawEtFullAcceptanceTPCNoPID( fRecAnalysis->GetRawEtFullAcceptanceTPCNoPID() );
 281         fMCAnalysis->FillSimRawEtVsRecoRawEtFullAcceptanceITSNoPID( fRecAnalysis->GetRawEtFullAcceptanceITSNoPID() );//Had
 282
 283       }
 284       if(fMCAnalysis->EMCAL()){
 285         fMCAnalysis->FillSimTotEtMinusRecoTotEtEMCALAcceptanceTPC( fRecAnalysis->GetCorrectedTotEtEMCALAcceptanceTPC() );
 286         fMCAnalysis->FillSimTotEtMinusRecoTotEtEMCALAcceptanceITS( fRecAnalysis->GetCorrectedTotEtEMCALAcceptanceITS() );
 287         fMCAnalysis->FillSimTotEtMinusRecoTotEtEMCALAcceptanceTPCNoPID( fRecAnalysis->GetCorrectedTotEtEMCALAcceptanceTPCNoPID() );
 288         fMCAnalysis->FillSimTotEtMinusRecoTotEtEMCALAcceptanceITSNoPID( fRecAnalysis->GetCorrectedTotEtEMCALAcceptanceITSNoPID() );
 289         fMCAnalysis->FillSimHadEtMinusRecoHadEtEMCALAcceptanceTPC( fRecAnalysis->GetCorrectedHadEtEMCALAcceptanceTPC() );
 290         fMCAnalysis->FillSimHadEtMinusRecoHadEtEMCALAcceptanceITS( fRecAnalysis->GetCorrectedHadEtEMCALAcceptanceITS() );
 291         fMCAnalysis->FillSimHadEtMinusRecoHadEtEMCALAcceptanceTPCNoPID( fRecAnalysis->GetCorrectedHadEtEMCALAcceptanceTPCNoPID() );
 292         fMCAnalysis->FillSimHadEtMinusRecoHadEtEMCALAcceptanceITSNoPID( fRecAnalysis->GetCorrectedHadEtEMCALAcceptanceITSNoPID() );
 293       }
 294       if(fMCAnalysis->PHOS()){
 295         fMCAnalysis->FillSimTotEtMinusRecoTotEtPHOSAcceptanceTPC( fRecAnalysis->GetCorrectedTotEtPHOSAcceptanceTPC() );
 296         fMCAnalysis->FillSimTotEtMinusRecoTotEtPHOSAcceptanceITS( fRecAnalysis->GetCorrectedTotEtPHOSAcceptanceITS() );
 297         fMCAnalysis->FillSimTotEtMinusRecoTotEtPHOSAcceptanceTPCNoPID( fRecAnalysis->GetCorrectedTotEtPHOSAcceptanceTPCNoPID() );
 298         fMCAnalysis->FillSimTotEtMinusRecoTotEtPHOSAcceptanceITSNoPID( fRecAnalysis->GetCorrectedTotEtPHOSAcceptanceITSNoPID() );
 299         fMCAnalysis->FillSimHadEtMinusRecoHadEtPHOSAcceptanceTPC( fRecAnalysis->GetCorrectedHadEtPHOSAcceptanceTPC() );
 300         fMCAnalysis->FillSimHadEtMinusRecoHadEtPHOSAcceptanceITS( fRecAnalysis->GetCorrectedHadEtPHOSAcceptanceITS() );
 301         fMCAnalysis->FillSimHadEtMinusRecoHadEtPHOSAcceptanceTPCNoPID( fRecAnalysis->GetCorrectedHadEtPHOSAcceptanceTPCNoPID() );
 302         fMCAnalysis->FillSimHadEtMinusRecoHadEtPHOSAcceptanceITSNoPID( fRecAnalysis->GetCorrectedHadEtPHOSAcceptanceITSNoPID() );
 303       }
 304       if(fMCAnalysis->PiKP() && fMCAnalysis->Full()){
 305         fMCAnalysis->FillSimPiKPMinusRecoPiKPFullAcceptanceTPC(fRecAnalysis->GetCorrectedPiKPEtFullAcceptanceTPC());
 306         fMCAnalysis->FillSimPiKPMinusRecoPiKPFullAcceptanceITS(fRecAnalysis->GetCorrectedPiKPEtFullAcceptanceITS());
 307         fMCAnalysis->FillSimPiKPMinusRecoPiKPFullAcceptanceTPCNoPID(fRecAnalysis->GetCorrectedPiKPEtFullAcceptanceTPCNoPID());
 308         fMCAnalysis->FillSimPiKPMinusRecoPiKPFullAcceptanceITSNoPID(fRecAnalysis->GetCorrectedPiKPEtFullAcceptanceITSNoPID());
 309       }
 310     }
 311   }
 312   }
 313   delete fTriggerAnalysis;
 314   //}
 315 //cout<<"End Event"<<endl<<endl;
 316 // Post output data.
 317  PostData(1, fOutputList);
 318 }
 319
 320 //________________________________________________________________________
 321 void AliAnalysisTaskHadEt::Terminate(Option_t *)
 322 {
 323     // Draw result to the screen
 324     // Called once at the end of the query
 325
 326     fOutputList = dynamic_cast<TList*> (GetOutputData(1));
 327     if (!fOutputList) {
 328         printf("ERROR: Output list not available\n");
 329         return;
 330     }
 331 }
 332
 333