PWG4/totEt/AliAnalysisHadEtReconstructed.cxx

   1 //_________________________________________________________________________
   2 //  Utility Class for transverse energy studies, charged hadrons
   3 //  Base class for ESD analysis
   4 //  - reconstruction output
   5 // implementation file
   6 //
   7 //Created by Christine Nattrass, Rebecca Scott, Irakli Martashvili
   8 //University of Tennessee at Knoxville
   9 //_________________________________________________________________________
  10
  11 #include <TROOT.h>
  12 #include <TSystem.h>
  13 #include <TInterpreter.h>
  14 #include "AliAnalysisHadEtReconstructed.h"
  15 #include "AliAnalysisEtCuts.h"
  16 #include "AliESDtrack.h"
  17 #include "AliESDCaloCluster.h"
  18 #include "AliVEvent.h"
  19 #include "AliESDEvent.h"
  20 #include "AliESDtrackCuts.h"
  21 #include "AliESDpid.h"
  22 #include "AliVParticle.h"
  23 #include <iostream>
  24 #include "AliAnalysisHadEtCorrections.h"
  25 #include "TFile.h"
  26 #include "TString.h"
  27 #include "AliAnalysisEtCommon.h"
  28 #include "AliAnalysisHadEt.h"
  29
  30 using namespace std;
  31
  32 ClassImp(AliAnalysisHadEtReconstructed);
  33
  34
  35 AliAnalysisHadEtReconstructed::AliAnalysisHadEtReconstructed() :
  36         AliAnalysisHadEt()
  37         ,fCorrections(0)
  38         ,fConfigFile("ConfigHadEtAnalysis.C")
  39         ,fCorrTotEtFullAcceptanceTPC(0)
  40         ,fCorrTotEtFullAcceptanceITS(0)
  41         ,fCorrHadEtFullAcceptanceTPC(0)
  42         ,fCorrHadEtFullAcceptanceITS(0)
  43         ,fCorrTotEtEMCALAcceptanceTPC(0)
  44         ,fCorrTotEtEMCALAcceptanceITS(0)
  45         ,fCorrHadEtEMCALAcceptanceTPC(0)
  46         ,fCorrHadEtEMCALAcceptanceITS(0)
  47         ,fCorrTotEtPHOSAcceptanceTPC(0)
  48         ,fCorrTotEtPHOSAcceptanceITS(0)
  49         ,fCorrHadEtPHOSAcceptanceTPC(0)
  50         ,fCorrHadEtPHOSAcceptanceITS(0)
  51         ,fCorrectedHadEtFullAcceptanceTPCNoPID(0)
  52         ,fCorrectedHadEtFullAcceptanceITSNoPID(0)
  53         ,fCorrectedHadEtEMCALAcceptanceTPCNoPID(0)
  54         ,fCorrectedHadEtEMCALAcceptanceITSNoPID(0)
  55         ,fCorrectedHadEtPHOSAcceptanceTPCNoPID(0)
  56         ,fCorrectedHadEtPHOSAcceptanceITSNoPID(0)
  57         ,fCorrectedHadEtFullAcceptanceTPC(0)
  58         ,fCorrectedHadEtFullAcceptanceITS(0)
  59         ,fCorrectedHadEtEMCALAcceptanceTPC(0)
  60         ,fCorrectedHadEtEMCALAcceptanceITS(0)
  61         ,fCorrectedHadEtPHOSAcceptanceTPC(0)
  62         ,fCorrectedHadEtPHOSAcceptanceITS(0)
  63         ,fRawEtFullAcceptanceTPC(0)
  64         ,fRawEtFullAcceptanceITS(0)
  65         ,fRawEtEMCALAcceptanceTPC(0)
  66         ,fRawEtEMCALAcceptanceITS(0)
  67         ,fRawEtPHOSAcceptanceTPC(0)
  68         ,fRawEtPHOSAcceptanceITS(0)
  69         ,fRawEtFullAcceptanceTPCNoPID(0)
  70         ,fRawEtFullAcceptanceITSNoPID(0)
  71         ,fRawEtEMCALAcceptanceTPCNoPID(0)
  72         ,fRawEtEMCALAcceptanceITSNoPID(0)
  73         ,fRawEtPHOSAcceptanceTPCNoPID(0)
  74         ,fRawEtPHOSAcceptanceITSNoPID(0)
  75 {
  76 }
  77
  78 AliAnalysisHadEtReconstructed::~AliAnalysisHadEtReconstructed()
  79 {
  80   delete fCorrections;
  81 }
  82
  83 Int_t AliAnalysisHadEtReconstructed::AnalyseEvent(AliVEvent* ev)
  84 { // analyse ESD event
  85     ResetEventValues();
  86
  87     AliESDEvent *realEvent = dynamic_cast<AliESDEvent*>(ev);
  88     //for PID
  89     AliESDpid *pID = new AliESDpid();
  90     pID->MakePID(realEvent);
  91     TString *strTPC = new TString("TPC");
  92     TString *strITS = new TString("ITS");
  93     TString *strTPCITS = new TString("TPCITS");
  94     for(Int_t cutset=0;cutset<2;cutset++){
  95       bool isTPC = false;
  96       TString *cutName;
  97       TObjArray* list;
  98       switch(cutset){
  99       case 0:
 100         cutName = strTPCITS;
 101         list = fEsdtrackCutsITSTPC->GetAcceptedTracks(realEvent);
 102         isTPC = true;
 103         break;
 104       case 1:
 105         cutName = strITS;
 106         list = fEsdtrackCutsITS->GetAcceptedTracks(realEvent);
 107         break;
 108       case 2:
 109         cutName = strTPC;
 110         list = fEsdtrackCutsTPC->GetAcceptedTracks(realEvent);
 111         break;
 112       default:
 113         cerr<<"Error:  cannot fill histograms!"<<endl;
 114         return -1;
 115       }
 116       Int_t nGoodTracks = list->GetEntries();
 117       for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++)
 118         {
 119
 120
 121           AliESDtrack *track = dynamic_cast<AliESDtrack*> (list->At(iTrack));
 122           if (!track)
 123             {
 124               Printf("ERROR: Could not get track %d", iTrack);
 125               continue;
 126             }
 127           else{
 128             if(TMath::Abs(track->Eta())>fCorrections->GetEtaCut()) continue;
 129             Float_t nSigmaPion,nSigmaProton,nSigmaKaon,nSigmaElectron;
 130             if(cutset!=1){
 131               nSigmaPion = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kPion));
 132               nSigmaProton = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kProton));
 133               nSigmaKaon = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kKaon));
 134               nSigmaElectron = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kElectron));
 135             }
 136             else{
 137               nSigmaPion = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kPion));
 138               nSigmaProton = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kProton));
 139               nSigmaKaon = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kKaon));
 140               nSigmaElectron = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kElectron));
 141             }
 142             bool isPion = (nSigmaPion<3.0 && nSigmaProton>2.0 && nSigmaKaon>2.0);
 143             bool isElectron = (nSigmaElectron<2.0 && nSigmaPion>4.0 && nSigmaProton>3.0 && nSigmaKaon>3.0);
 144             bool isKaon = (nSigmaPion>3.0 && nSigmaProton>2.0 && nSigmaKaon<2.0);
 145             bool isProton = (nSigmaPion>3.0 && nSigmaProton<2.0 && nSigmaKaon>2.0);
 146
 147             bool unidentified = (!isProton && !isKaon && !isElectron && !isPion);
 148             Float_t dEdx = track->GetTPCsignal();
 149             if(cutset==1) dEdx = track->GetITSsignal();
 150             FillHisto2D(Form("dEdxDataAll%s",cutName->Data()),track->P(),dEdx,1.0);
 151
 152             bool inPHOS = IsInPHOS(track);
 153             bool inEMCAL = IsInEMCAL(track);
 154
 155             Float_t corrBkgd=0.0;
 156             Float_t corrNotID=0.0;
 157             Float_t corrNoID=0.0;// = fCorrections->GetNotIDCorrectionNoPID(track->Pt());
 158             Float_t corrEff = 0.0;
 159             Float_t corrEffNoID = 0.0;
 160             if(cutset==0){//TPC
 161               corrBkgd = fCorrections->GetBackgroundCorrectionTPC(track->Pt());
 162               corrEffNoID = fCorrections->GetTPCEfficiencyCorrectionHadron(track->Pt());
 163               corrNotID = fCorrections->GetNotIDConstCorrectionTPC();
 164               corrNoID = fCorrections->GetNotIDConstCorrectionTPCNoID();
 165             }
 166             if(cutset==1){//ITS
 167               corrBkgd = fCorrections->GetBackgroundCorrectionITS(track->Pt());
 168               corrEffNoID = fCorrections->GetITSEfficiencyCorrectionHadron(track->Pt());
 169               corrNotID = fCorrections->GetNotIDConstCorrectionITS();
 170               corrNoID = fCorrections->GetNotIDConstCorrectionITSNoID();
 171             }
 172             Float_t et = 0.0;
 173             Float_t etNoID = Et(track->P(),track->Theta(),fgPiPlusCode,track->Charge());
 174             Float_t etpartialcorrected = 0.0;
 175             Float_t etpartialcorrectedNoID = corrNoID*corrBkgd*corrEffNoID*etNoID;
 176             FillHisto2D(Form("EtDataRaw%sNoID",cutName->Data()),track->Pt(),track->Eta(),etpartialcorrectedNoID);
 177
 178             if(isPion){
 179               FillHisto2D(Form("dEdxDataPion%s",cutName->Data()),track->P(),dEdx,1.0);
 180               et = Et(track->P(),track->Theta(),fgPiPlusCode,track->Charge());
 181               if(cutset==0){corrEff = fCorrections->GetTPCEfficiencyCorrectionPion(track->Pt());}
 182               //else{corrEff = fCorrections->GetITSEfficiencyCorrectionPion(track->Pt());}
 183               etpartialcorrected = et*corrBkgd*corrEff*corrNotID;
 184
 185               if(track->Charge()>0.0){
 186                 FillHisto2D(Form("EtDataRaw%sPiPlus",cutName->Data()),track->Pt(),track->Eta(),et);
 187                 FillHisto2D(Form("EtDataCorrected%sPiPlus",cutName->Data()),track->Pt(),track->Eta(),etpartialcorrected);
 188               }
 189               else{
 190                 FillHisto2D(Form("EtDataRaw%sPiMinus",cutName->Data()),track->Pt(),track->Eta(),et);
 191                 FillHisto2D(Form("EtDataCorrected%sPiMinus",cutName->Data()),track->Pt(),track->Eta(),etpartialcorrected);
 192               }
 193             }
 194             if(isKaon){
 195               FillHisto2D(Form("dEdxDataKaon%s",cutName->Data()),track->P(),dEdx,1.0);
 196               et = Et(track->P(),track->Theta(),fgKPlusCode,track->Charge());
 197               if(cutset==0){corrEff = fCorrections->GetTPCEfficiencyCorrectionKaon(track->Pt());}
 198               //else{corrEff = fCorrections->GetITSEfficiencyCorrectionKaon(track->Pt());}
 199               etpartialcorrected = et*corrBkgd*corrEff*corrNotID;
 200
 201               if(track->Charge()>0.0){
 202                 FillHisto2D(Form("EtDataRaw%sKPlus",cutName->Data()),track->Pt(),track->Eta(),et);
 203                 FillHisto2D(Form("EtDataCorrected%sKPlus",cutName->Data()),track->Pt(),track->Eta(),etpartialcorrected);
 204               }
 205               else{
 206                 FillHisto2D(Form("EtDataRaw%sKMinus",cutName->Data()),track->Pt(),track->Eta(),et);
 207                 FillHisto2D(Form("EtDataCorrected%sKMinus",cutName->Data()),track->Pt(),track->Eta(),etpartialcorrected);
 208               }
 209             }
 210             if(isProton){
 211               FillHisto2D(Form("dEdxDataProton%s",cutName->Data()),track->P(),dEdx,1.0);
 212               et = Et(track->P(),track->Theta(),fgProtonCode,track->Charge());
 213               if(cutset==0){corrEff = fCorrections->GetTPCEfficiencyCorrectionProton(track->Pt());}
 214               //else{corrEff = fCorrections->GetITSEfficiencyCorrectionProton(track->Pt());}
 215               etpartialcorrected = et*corrBkgd*corrEff*corrNotID;
 216
 217               if(track->Charge()>0.0){
 218                 FillHisto2D(Form("EtDataRaw%sProton",cutName->Data()),track->Pt(),track->Eta(),et);
 219                 FillHisto2D(Form("EtDataCorrected%sProton",cutName->Data()),track->Pt(),track->Eta(),etpartialcorrected);
 220               }
 221               else{
 222                 FillHisto2D(Form("EtDataRaw%sAntiProton",cutName->Data()),track->Pt(),track->Eta(),et);
 223                 FillHisto2D(Form("EtDataCorrected%sAntiProton",cutName->Data()),track->Pt(),track->Eta(),etpartialcorrected);
 224               }
 225             }
 226             if(isElectron){
 227               FillHisto2D(Form("dEdxDataProton%s",cutName->Data()),track->P(),dEdx,1.0);
 228               //et = Et(track->P(),track->Theta(),fgPiPlusCode,track->Charge());
 229             }
 230             if(unidentified){
 231               //if(!isPion)
 232               FillHisto2D(Form("dEdxDataUnidentified%s",cutName->Data()),track->P(),dEdx,1.0);
 233               et = Et(track->P(),track->Theta(),fgPiPlusCode,track->Charge());
 234               etpartialcorrected = et*corrBkgd*corrEffNoID*corrNotID;
 235               //if(!isPion)
 236               FillHisto2D(Form("EtDataCorrected%sUnidentified",cutName->Data()),track->Pt(),track->Eta(),etpartialcorrected);
 237             }
 238             if(!isTPC){
 239               etpartialcorrected = etpartialcorrectedNoID;//Not using PID for ITS
 240             }
 241             AddEt(et,etNoID,etpartialcorrected,etpartialcorrectedNoID,track->Pt(),isTPC,inPHOS,inEMCAL);
 242           }
 243         }
 244     }
 245     if(GetCorrectedHadEtFullAcceptanceTPC()>0.0)FillHisto1D("RecoHadEtFullAcceptanceTPC",GetCorrectedHadEtFullAcceptanceTPC(),1.0);
 246     if(GetCorrectedTotEtFullAcceptanceTPC()>0.0)FillHisto1D("RecoTotEtFullAcceptanceTPC",GetCorrectedTotEtFullAcceptanceTPC(),1.0);
 247     if(GetCorrectedHadEtEMCALAcceptanceTPC()>0.0)FillHisto1D("RecoHadEtEMCALAcceptanceTPC",GetCorrectedHadEtEMCALAcceptanceTPC(),1.0);
 248     if(GetCorrectedTotEtEMCALAcceptanceTPC()>0.0)FillHisto1D("RecoTotEtEMCALAcceptanceTPC",GetCorrectedTotEtEMCALAcceptanceTPC(),1.0);
 249     if(GetCorrectedHadEtPHOSAcceptanceTPC()>0.0)FillHisto1D("RecoHadEtPHOSAcceptanceTPC",GetCorrectedHadEtPHOSAcceptanceTPC(),1.0);
 250     if(GetCorrectedTotEtPHOSAcceptanceTPC()>0.0)FillHisto1D("RecoTotEtPHOSAcceptanceTPC",GetCorrectedTotEtPHOSAcceptanceTPC(),1.0);
 251     if(GetCorrectedHadEtFullAcceptanceTPCNoPID()>0.0)FillHisto1D("RecoHadEtFullAcceptanceTPCNoPID",GetCorrectedHadEtFullAcceptanceTPCNoPID(),1.0);
 252     if(GetCorrectedTotEtFullAcceptanceTPCNoPID()>0.0)FillHisto1D("RecoTotEtFullAcceptanceTPCNoPID",GetCorrectedTotEtFullAcceptanceTPCNoPID(),1.0);
 253     if(GetCorrectedHadEtEMCALAcceptanceTPCNoPID()>0.0)FillHisto1D("RecoHadEtEMCALAcceptanceTPCNoPID",GetCorrectedHadEtEMCALAcceptanceTPCNoPID(),1.0);
 254     if(GetCorrectedTotEtEMCALAcceptanceTPCNoPID()>0.0)FillHisto1D("RecoTotEtEMCALAcceptanceTPCNoPID",GetCorrectedTotEtEMCALAcceptanceTPCNoPID(),1.0);
 255     if(GetCorrectedHadEtPHOSAcceptanceTPCNoPID()>0.0)FillHisto1D("RecoHadEtPHOSAcceptanceTPCNoPID",GetCorrectedHadEtPHOSAcceptanceTPCNoPID(),1.0);
 256     if(GetCorrectedTotEtPHOSAcceptanceTPCNoPID()>0.0)FillHisto1D("RecoTotEtPHOSAcceptanceTPCNoPID",GetCorrectedTotEtPHOSAcceptanceTPCNoPID(),1.0);
 257     if(GetCorrectedHadEtFullAcceptanceITS()>0.0)FillHisto1D("RecoHadEtFullAcceptanceITS",GetCorrectedHadEtFullAcceptanceITS(),1.0);
 258     if(GetCorrectedTotEtFullAcceptanceITS()>0.0)FillHisto1D("RecoTotEtFullAcceptanceITS",GetCorrectedTotEtFullAcceptanceITS(),1.0);
 259     if(GetCorrectedHadEtEMCALAcceptanceITS()>0.0)FillHisto1D("RecoHadEtEMCALAcceptanceITS",GetCorrectedHadEtEMCALAcceptanceITS(),1.0);
 260     if(GetCorrectedTotEtEMCALAcceptanceITS()>0.0)FillHisto1D("RecoTotEtEMCALAcceptanceITS",GetCorrectedTotEtEMCALAcceptanceITS(),1.0);
 261     if(GetCorrectedHadEtPHOSAcceptanceITS()>0.0)FillHisto1D("RecoHadEtPHOSAcceptanceITS",GetCorrectedHadEtPHOSAcceptanceITS(),1.0);
 262     if(GetCorrectedTotEtPHOSAcceptanceITS()>0.0)FillHisto1D("RecoTotEtPHOSAcceptanceITS",GetCorrectedTotEtPHOSAcceptanceITS(),1.0);
 263     if(GetCorrectedHadEtFullAcceptanceITSNoPID()>0.0)FillHisto1D("RecoHadEtFullAcceptanceITSNoPID",GetCorrectedHadEtFullAcceptanceITSNoPID(),1.0);
 264     if(GetCorrectedTotEtFullAcceptanceITSNoPID()>0.0)FillHisto1D("RecoTotEtFullAcceptanceITSNoPID",GetCorrectedTotEtFullAcceptanceITSNoPID(),1.0);
 265     if(GetCorrectedHadEtEMCALAcceptanceITSNoPID()>0.0)FillHisto1D("RecoHadEtEMCALAcceptanceITSNoPID",GetCorrectedHadEtEMCALAcceptanceITSNoPID(),1.0);
 266     if(GetCorrectedTotEtEMCALAcceptanceITSNoPID()>0.0)FillHisto1D("RecoTotEtEMCALAcceptanceITSNoPID",GetCorrectedTotEtEMCALAcceptanceITSNoPID(),1.0);
 267     if(GetCorrectedHadEtPHOSAcceptanceITSNoPID()>0.0)FillHisto1D("RecoHadEtPHOSAcceptanceITSNoPID",GetCorrectedHadEtPHOSAcceptanceITSNoPID(),1.0);
 268     if(GetCorrectedTotEtPHOSAcceptanceITSNoPID()>0.0)FillHisto1D("RecoTotEtPHOSAcceptanceITSNoPID",GetCorrectedTotEtPHOSAcceptanceITSNoPID(),1.0);
 269
 270     if(GetRawEtFullAcceptanceTPC()>0.0)FillHisto1D("RecoRawEtFullAcceptanceTPC",GetRawEtFullAcceptanceTPC(),1.0);
 271     if(GetRawEtEMCALAcceptanceTPC()>0.0)FillHisto1D("RecoRawEtEMCALAcceptanceTPC",GetRawEtEMCALAcceptanceTPC(),1.0);
 272     if(GetRawEtPHOSAcceptanceTPC()>0.0)FillHisto1D("RecoRawEtPHOSAcceptanceTPC",GetRawEtPHOSAcceptanceTPC(),1.0);
 273     if(GetRawEtFullAcceptanceTPCNoPID()>0.0)FillHisto1D("RecoRawEtFullAcceptanceTPCNoPID",GetRawEtFullAcceptanceTPCNoPID(),1.0);
 274     if(GetRawEtEMCALAcceptanceTPCNoPID()>0.0)FillHisto1D("RecoRawEtEMCALAcceptanceTPCNoPID",GetRawEtEMCALAcceptanceTPCNoPID(),1.0);
 275     if(GetRawEtPHOSAcceptanceTPCNoPID()>0.0)FillHisto1D("RecoRawEtPHOSAcceptanceTPCNoPID",GetRawEtPHOSAcceptanceTPCNoPID(),1.0);
 276     if(GetRawEtFullAcceptanceITS()>0.0)FillHisto1D("RecoRawEtFullAcceptanceITS",GetRawEtFullAcceptanceITS(),1.0);
 277     if(GetRawEtEMCALAcceptanceITS()>0.0)FillHisto1D("RecoRawEtEMCALAcceptanceITS",GetRawEtEMCALAcceptanceITS(),1.0);
 278     if(GetRawEtPHOSAcceptanceITS()>0.0)FillHisto1D("RecoRawEtPHOSAcceptanceITS",GetRawEtPHOSAcceptanceITS(),1.0);
 279     if(GetRawEtFullAcceptanceITSNoPID()>0.0)FillHisto1D("RecoRawEtFullAcceptanceITSNoPID",GetRawEtFullAcceptanceITSNoPID(),1.0);
 280     if(GetRawEtEMCALAcceptanceITSNoPID()>0.0)FillHisto1D("RecoRawEtEMCALAcceptanceITSNoPID",GetRawEtEMCALAcceptanceITSNoPID(),1.0);
 281     if(GetRawEtPHOSAcceptanceITSNoPID()>0.0)FillHisto1D("RecoRawEtPHOSAcceptanceITSNoPID",GetRawEtPHOSAcceptanceITSNoPID(),1.0);
 282     delete pID;
 283     delete strTPC;
 284     delete strITS;
 285     delete strTPCITS;
 286     return 1;
 287 }
 288 void AliAnalysisHadEtReconstructed::AddEt(Float_t rawEt, Float_t rawEtNoPID, Float_t corrEt, Float_t corrEtNoPID, Float_t pt, Bool_t IsTPC, Bool_t InPHOS, Bool_t InEMCAL) {//Adding Et to each of the variables that tracks et event by event
 289   if(pt>=AliAnalysisHadEt::fgPtTPCCutOff && IsTPC){//TPC tracks
 290     //adding to the raw Et
 291     fRawEtFullAcceptanceTPC += rawEt;
 292     if(InPHOS)fRawEtPHOSAcceptanceTPC += rawEt;
 293     if(InEMCAL)fRawEtEMCALAcceptanceTPC += rawEt;
 294     fRawEtFullAcceptanceTPCNoPID += rawEtNoPID;
 295     if(InPHOS)fRawEtPHOSAcceptanceTPCNoPID += rawEtNoPID;
 296     if(InEMCAL)fRawEtEMCALAcceptanceTPCNoPID += rawEtNoPID;
 297     //adding to the corrected Et
 298     fCorrectedHadEtFullAcceptanceTPC += corrEt;
 299     if(InPHOS)fCorrectedHadEtPHOSAcceptanceTPC += corrEt;
 300     if(InEMCAL)fCorrectedHadEtEMCALAcceptanceTPC += corrEt;
 301     fCorrectedHadEtFullAcceptanceTPCNoPID += corrEtNoPID;
 302     if(InPHOS)fCorrectedHadEtPHOSAcceptanceTPCNoPID += corrEtNoPID;
 303     if(InEMCAL)fCorrectedHadEtEMCALAcceptanceTPCNoPID += corrEtNoPID;
 304   }
 305   if(pt<AliAnalysisHadEt::fgPtTPCCutOff &&pt>=AliAnalysisHadEt::fgPtITSCutOff && !IsTPC){//ITS tracks
 306     //adding to the raw Et
 307     fRawEtFullAcceptanceITS += rawEt;
 308     if(InPHOS)fRawEtPHOSAcceptanceITS += rawEt;
 309     if(InEMCAL)fRawEtEMCALAcceptanceITS += rawEt;
 310     fRawEtFullAcceptanceITSNoPID += rawEtNoPID;
 311     if(InPHOS)fRawEtPHOSAcceptanceITSNoPID += rawEtNoPID;
 312     if(InEMCAL)fRawEtEMCALAcceptanceITSNoPID += rawEtNoPID;
 313     //adding to the corrected Et
 314     fCorrectedHadEtFullAcceptanceITS += corrEt;
 315     if(InPHOS)fCorrectedHadEtPHOSAcceptanceITS += corrEt;
 316     if(InEMCAL)fCorrectedHadEtEMCALAcceptanceITS += corrEt;
 317     fCorrectedHadEtFullAcceptanceITSNoPID += corrEtNoPID;
 318     if(InPHOS)fCorrectedHadEtPHOSAcceptanceITSNoPID += corrEtNoPID;
 319     if(InEMCAL)fCorrectedHadEtEMCALAcceptanceITSNoPID += corrEtNoPID;
 320   }
 321 }
 322
 323 Bool_t AliAnalysisHadEtReconstructed::IsInPHOS(AliESDtrack *track){//This function will need to be elaborated on later to include PHOS dead channels
 324   return   TMath::Abs(track->Eta()) < fCuts->GetGeometryPhosEtaAccCut()//in eta acceptance
 325     && track->Phi()*180.0/TMath::Pi() > fCuts->GetGeometryPhosPhiAccMinCut()//greater than the minimum phi
 326     && track->Phi()*180.0/TMath::Pi() < fCuts->GetGeometryPhosPhiAccMaxCut();//less than the maximum phi
 327 }
 328 Bool_t AliAnalysisHadEtReconstructed::IsInEMCAL(AliESDtrack *track){//This function will need to be elaborated on later to include EMCAL dead channels
 329   return   TMath::Abs(track->Eta()) < fCuts->GetGeometryEmcalEtaAccCut()//in eta acceptance
 330     && track->Phi()*180.0/TMath::Pi() > fCuts->GetGeometryEmcalPhiAccMinCut()//greater than the minimum phi
 331     && track->Phi()*180.0/TMath::Pi() < fCuts->GetGeometryEmcalPhiAccMaxCut();//less than the maximum phi
 332 }
 333 Bool_t AliAnalysisHadEtReconstructed::CheckGoodVertex(AliVParticle* track)
 334 { // check vertex
 335
 336     Float_t bxy = 999.;
 337     Float_t bz = 999.;
 338     dynamic_cast<AliESDtrack*>(track)->GetImpactParametersTPC(bxy,bz);
 339
 340     bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) &&
 341       (TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) &&
 342       (TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) &&
 343       (TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) &&
 344       (TMath::Abs(bz) < fCuts->GetReconstructedIPzCut());
 345
 346     return status;
 347 }
 348
 349 void AliAnalysisHadEtReconstructed::Init()
 350 { // Init
 351   AliAnalysisHadEt::Init();
 352   if(fCorrections){
 353     fCorrTotEtFullAcceptanceTPC = fCorrections->GetConstantCorrections(kTRUE,fgPtTPCCutOff,"Full");
 354     fCorrTotEtFullAcceptanceITS = fCorrections->GetConstantCorrections(kTRUE,fgPtITSCutOff,"Full");
 355     fCorrHadEtFullAcceptanceTPC = fCorrections->GetConstantCorrections(kFALSE,fgPtTPCCutOff,"Full");
 356     fCorrHadEtFullAcceptanceITS = fCorrections->GetConstantCorrections(kFALSE,fgPtITSCutOff,"Full");
 357     fCorrTotEtEMCALAcceptanceTPC = fCorrections->GetConstantCorrections(kTRUE,fgPtTPCCutOff,"EMCAL");
 358     fCorrTotEtEMCALAcceptanceITS = fCorrections->GetConstantCorrections(kTRUE,fgPtITSCutOff,"EMCAL");
 359     fCorrHadEtEMCALAcceptanceTPC = fCorrections->GetConstantCorrections(kFALSE,fgPtTPCCutOff,"EMCAL");
 360     fCorrHadEtEMCALAcceptanceITS = fCorrections->GetConstantCorrections(kFALSE,fgPtITSCutOff,"EMCAL");
 361     fCorrTotEtPHOSAcceptanceTPC = fCorrections->GetConstantCorrections(kTRUE,fgPtTPCCutOff,"PHOS");
 362     fCorrTotEtPHOSAcceptanceITS = fCorrections->GetConstantCorrections(kTRUE,fgPtITSCutOff,"PHOS");
 363     fCorrHadEtPHOSAcceptanceTPC = fCorrections->GetConstantCorrections(kFALSE,fgPtTPCCutOff,"PHOS");
 364     fCorrHadEtPHOSAcceptanceITS = fCorrections->GetConstantCorrections(kFALSE,fgPtITSCutOff,"PHOS");
 365   }
 366   else{
 367     cout<<"Warning!  You have not set corrections.  Your code will crash.  You have to set the corrections."<<endl;
 368   }
 369 }
 370
 371 void AliAnalysisHadEtReconstructed::ResetEventValues(){//resetting event by event et's
 372   AliAnalysisHadEt::ResetEventValues();
 373      fCorrectedHadEtFullAcceptanceTPCNoPID=0.0;
 374      fCorrectedHadEtFullAcceptanceITSNoPID=0.0;
 375      fCorrectedHadEtEMCALAcceptanceTPCNoPID=0.0;
 376      fCorrectedHadEtEMCALAcceptanceITSNoPID=0.0;
 377      fCorrectedHadEtPHOSAcceptanceTPCNoPID=0.0;
 378      fCorrectedHadEtPHOSAcceptanceITSNoPID=0.0;
 379      fCorrectedHadEtFullAcceptanceTPC=0.0;
 380      fCorrectedHadEtFullAcceptanceITS=0.0;
 381      fCorrectedHadEtEMCALAcceptanceTPC=0.0;
 382      fCorrectedHadEtEMCALAcceptanceITS=0.0;
 383      fCorrectedHadEtPHOSAcceptanceTPC=0.0;
 384      fCorrectedHadEtPHOSAcceptanceITS=0.0;
 385      fRawEtFullAcceptanceTPC=0.0;
 386      fRawEtFullAcceptanceITS=0.0;
 387      fRawEtEMCALAcceptanceTPC=0.0;
 388      fRawEtEMCALAcceptanceITS=0.0;
 389      fRawEtPHOSAcceptanceTPC=0.0;
 390      fRawEtPHOSAcceptanceITS=0.0;
 391      fRawEtFullAcceptanceTPCNoPID=0.0;
 392      fRawEtFullAcceptanceITSNoPID=0.0;
 393      fRawEtEMCALAcceptanceTPCNoPID=0.0;
 394      fRawEtEMCALAcceptanceITSNoPID=0.0;
 395      fRawEtPHOSAcceptanceTPCNoPID=0.0;
 396      fRawEtPHOSAcceptanceITSNoPID=0.0;
 397
 398      if(TMath::Abs(fCorrTotEtFullAcceptanceTPC)<1e-3){
 399        if (fConfigFile.Length()) {
 400          cout<<"Warning: Rereading fCorrections file..."<<endl;
 401          gROOT->LoadMacro(fConfigFile);
 402          fCorrections = (AliAnalysisHadEtCorrections *) gInterpreter->ProcessLine("ConfigHadEtAnalysis()");
 403          fCorrTotEtFullAcceptanceTPC = fCorrections->GetConstantCorrections(kTRUE,fgPtTPCCutOff,"Full");
 404          fCorrTotEtFullAcceptanceITS = fCorrections->GetConstantCorrections(kTRUE,fgPtITSCutOff,"Full");
 405          fCorrHadEtFullAcceptanceTPC = fCorrections->GetConstantCorrections(kFALSE,fgPtTPCCutOff,"Full");
 406          fCorrHadEtFullAcceptanceITS = fCorrections->GetConstantCorrections(kFALSE,fgPtITSCutOff,"Full");
 407          fCorrTotEtEMCALAcceptanceTPC = fCorrections->GetConstantCorrections(kTRUE,fgPtTPCCutOff,"EMCAL");
 408          fCorrTotEtEMCALAcceptanceITS = fCorrections->GetConstantCorrections(kTRUE,fgPtITSCutOff,"EMCAL");
 409          fCorrHadEtEMCALAcceptanceTPC = fCorrections->GetConstantCorrections(kFALSE,fgPtTPCCutOff,"EMCAL");
 410          fCorrHadEtEMCALAcceptanceITS = fCorrections->GetConstantCorrections(kFALSE,fgPtITSCutOff,"EMCAL");
 411          fCorrTotEtPHOSAcceptanceTPC = fCorrections->GetConstantCorrections(kTRUE,fgPtTPCCutOff,"PHOS");
 412          fCorrTotEtPHOSAcceptanceITS = fCorrections->GetConstantCorrections(kTRUE,fgPtITSCutOff,"PHOS");
 413          fCorrHadEtPHOSAcceptanceTPC = fCorrections->GetConstantCorrections(kFALSE,fgPtTPCCutOff,"PHOS");
 414          fCorrHadEtPHOSAcceptanceITS = fCorrections->GetConstantCorrections(kFALSE,fgPtITSCutOff,"PHOS");
 415        }
 416        else{cerr<<"Uh-oh!  Unable to open configuration file!"<<endl;}
 417      }
 418 }
 419 void AliAnalysisHadEtReconstructed::CreateHistograms(){//Creating histograms and adding them to the output TList
 420
 421   TString *strTPC = new TString("TPC");
 422   TString *strITS = new TString("ITS");
 423   TString *strTPCITS = new TString("TPCITS");
 424   for(Int_t i=0;i<2;i++){
 425     TString *cutName;
 426     Float_t maxPtdEdx = 10;
 427     Float_t mindEdx = 35;
 428     Float_t maxdEdx = 150.0;
 429     switch(i){
 430     case 0:
 431       cutName = strTPCITS;
 432       break;
 433     case 1:
 434       cutName = strITS;
 435       maxPtdEdx = 5;
 436       maxdEdx = 500.0;
 437       break;
 438     case 2:
 439       cutName = strTPC;
 440       maxPtdEdx = 5;
 441       maxdEdx = 500.0;
 442       break;
 443     default:
 444       cerr<<"Error:  cannot make histograms!"<<endl;
 445       return;
 446     }
 447
 448     CreateEtaPtHisto2D(Form("EtDataRaw%sPiPlus",cutName->Data()),"Raw reconstructed E_{T} from identified #pi^{+}");
 449     CreateEtaPtHisto2D(Form("EtDataRaw%sPiMinus",cutName->Data()),"Raw reconstructed E_{T} from identified #pi^{-}");
 450     CreateEtaPtHisto2D(Form("EtDataRaw%sKPlus",cutName->Data()),"Raw reconstructed E_{T} from identified K^{+}");
 451 //     CreateEtaPtHisto2D(Form("EtDataRaw%sEMinus",cutName->Data()),"Raw reconstructed E_{T} from identified e^{-}");
 452 //     CreateEtaPtHisto2D(Form("EtDataRaw%sEPlus",cutName->Data()),"Raw reconstructed E_{T} from identified e^{+}");
 453     CreateEtaPtHisto2D(Form("EtDataRaw%sKMinus",cutName->Data()),"Raw reconstructed E_{T} from identified K^{-}");
 454     CreateEtaPtHisto2D(Form("EtDataRaw%sProton",cutName->Data()),"Raw reconstructed E_{T} from identified p");
 455     CreateEtaPtHisto2D(Form("EtDataRaw%sAntiProton",cutName->Data()),"Raw reconstructed E_{T} from identified #bar{p}");
 456     CreateEtaPtHisto2D(Form("EtDataRaw%sUnidentified",cutName->Data()),"Raw reconstructed E_{T} from unidentified particles using real mass");
 457     CreateEtaPtHisto2D(Form("EtDataRaw%sNoID",cutName->Data()),"Raw reconstructed E_{T} from unidentified particles using real mass");
 458
 459     CreateEtaPtHisto2D(Form("EtDataCorrected%sPiPlus",cutName->Data()),"Corrected reconstructed E_{T} from identified #pi^{+}");
 460     CreateEtaPtHisto2D(Form("EtDataCorrected%sPiMinus",cutName->Data()),"Corrected reconstructed E_{T} from identified #pi^{-}");
 461     CreateEtaPtHisto2D(Form("EtDataCorrected%sKPlus",cutName->Data()),"Corrected reconstructed E_{T} from identified K^{+}");
 462 //     CreateEtaPtHisto2D(Form("EtDataCorrected%sEMinus",cutName->Data()),"Corrected reconstructed E_{T} from identified e^{-}");
 463 //     CreateEtaPtHisto2D(Form("EtDataCorrected%sEPlus",cutName->Data()),"Corrected reconstructed E_{T} from identified e^{+}");
 464     CreateEtaPtHisto2D(Form("EtDataCorrected%sKMinus",cutName->Data()),"Corrected reconstructed E_{T} from identified K^{-}");
 465     CreateEtaPtHisto2D(Form("EtDataCorrected%sProton",cutName->Data()),"Corrected reconstructed E_{T} from identified p");
 466     CreateEtaPtHisto2D(Form("EtDataCorrected%sAntiProton",cutName->Data()),"Corrected reconstructed E_{T} from identified #bar{p}");
 467     CreateEtaPtHisto2D(Form("EtDataCorrected%sUnidentified",cutName->Data()),"Corrected reconstructed E_{T} from unidentified particles using real mass");
 468     CreateEtaPtHisto2D(Form("EtDataCorrected%sNoID",cutName->Data()),"Corrected reconstructed E_{T} from unidentified particles using real mass");
 469
 470
 471     CreateEtaPtHisto2D(Form("EtNData%sPiPlus",cutName->Data()),"Number of reconstructed #pi^{+}");
 472     CreateEtaPtHisto2D(Form("EtNData%sPiMinus",cutName->Data()),"Number of reconstructed #pi^{-}");
 473     CreateEtaPtHisto2D(Form("EtNData%sKPlus",cutName->Data()),"Number of reconstructed K^{+}");
 474     CreateEtaPtHisto2D(Form("EtNData%sKMinus",cutName->Data()),"Number of reconstructed K^{-}");
 475     CreateEtaPtHisto2D(Form("EtNData%sProton",cutName->Data()),"Number of reconstructed p");
 476     CreateEtaPtHisto2D(Form("EtNData%sAntiProton",cutName->Data()),"Number of reconstructed #bar{p}");
 477     CreateEtaPtHisto2D(Form("EtNData%sUnidentified",cutName->Data()),"Number of Reconstructed unidentified particles");
 478
 479     CreateHisto2D(Form("dEdxDataAll%s",cutName->Data()),"dE/dx for all particles","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
 480     CreateHisto2D(Form("dEdxDataPion%s",cutName->Data()),"dE/dx for #pi^{#pm}","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
 481     CreateHisto2D(Form("dEdxDataKaon%s",cutName->Data()),"dE/dx for K^{#pm}","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
 482     CreateHisto2D(Form("dEdxDataProton%s",cutName->Data()),"dE/dx for p(#bar{p})","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
 483     CreateHisto2D(Form("dEdxDataElectron%s",cutName->Data()),"dE/dx for e^{#pm}","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
 484     CreateHisto2D(Form("dEdxDataUnidentified%s",cutName->Data()),"dE/dx for unidentified particles","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
 485   }
 486
 487   Float_t minEt = 0.0;
 488   Float_t maxEt = 100.0;
 489   Int_t nbinsEt = 200;
 490   char histoname[200];
 491   char histotitle[200];
 492   char xtitle[50];
 493   TString *ytitle = new TString("Number of events");
 494   TString *sTPC = new TString("TPC");
 495   TString *sITS = new TString("ITS");
 496   TString *sTPCpt = new TString("0.15");
 497   TString *sITSpt = new TString("0.10");
 498   TString *sPID = new TString("");
 499   TString *sNoPID = new TString("NoPID");
 500   TString *sNoPIDString = new TString(", No PID");
 501   TString *sHadEt = new TString("HadEt");
 502   TString *sRawEt = new TString("RawEt");
 503   TString *sTotEt = new TString("TotEt");
 504   TString *sTotEtString = new TString("total E_{T}");
 505   TString *sHadEtString = new TString("hadronic E_{T}");
 506   TString *sRawEtString = new TString("raw E_{T}");
 507   TString *sFull = new TString("Full");
 508   TString *sEMCAL = new TString("EMCAL");
 509   TString *sPHOS = new TString("PHOS");
 510
 511   for(int tpc = 0;tpc<2;tpc++){
 512     for(int hadet = 0;hadet<3;hadet++){
 513       for(int type = 0;type<3;type++){
 514         for(int pid = 0;pid<2;pid++){
 515           TString *detector;
 516           TString *partid;
 517           TString *et = sHadEt;
 518           TString *acceptance;
 519           TString *ptstring;
 520           TString *partidstring;
 521           TString *etstring = sHadEtString;
 522           if(tpc==1) {detector = sTPC; ptstring = sTPCpt;}
 523           else{detector = sITS; ptstring = sITSpt;}
 524           if(pid==1){partid = sPID; partidstring = sPID;}
 525           else{partid = sNoPID; partidstring = sNoPIDString;}
 526           if(hadet==1) {et = sHadEt; etstring = sHadEtString;}
 527           if(hadet==0){et = sTotEt; etstring = sTotEtString;}
 528           if(hadet==2){et = sRawEt; etstring = sRawEtString;}
 529           switch(type){
 530           case 0:
 531             acceptance = sFull;
 532             break;
 533           case 1:
 534             acceptance = sEMCAL;
 535             break;
 536           case 2:
 537             acceptance = sPHOS;
 538             break;
 539           default:
 540             acceptance = sFull;
 541           }
 542           sprintf(histoname,"Reco%s%sAcceptance%s%s",et->Data(),acceptance->Data(),detector->Data(),partid->Data());
 543           sprintf(histotitle,"Reconstructed %s with %s acceptance for p_{T}>%s GeV/c%s",etstring->Data(),acceptance->Data(),ptstring->Data(),partidstring->Data());
 544           sprintf(xtitle,"Reconstructed %s",etstring->Data());
 545           CreateHisto1D(histoname,histotitle,xtitle,ytitle->Data(),nbinsEt*2,minEt,maxEt);
 546         }
 547       }
 548     }
 549   }
 550
 551   //CreateHisto2D("Efficiency","Efficiency","pT","efficiency",
 552
 553   delete sTPC;
 554   delete sITS;
 555   delete sTPCpt;
 556   delete sITSpt;
 557   delete sPID;
 558   delete sNoPID;
 559   delete sNoPIDString;
 560   delete sHadEt;
 561   delete sTotEt;
 562   delete sTotEtString;
 563   delete sHadEtString;
 564   delete sFull;
 565   delete sEMCAL;
 566   delete sPHOS;
 567
 568 }