PWG4/totEt/AliAnalysisHadEtMonteCarlo.cxx

   1 //_________________________________________________________________________
   2 //  Utility Class for transverse energy studies, charged hadrons
   3 //  Base class for MC analysis
   4 //  - MC output
   5 // implementation file
   6 //
   7 //Created by Christine Nattrass, Rebecca Scott, Irakli Martashvili
   8 //University of Tennessee at Knoxville
   9 //_________________________________________________________________________
  10 #include "AliAnalysisHadEtMonteCarlo.h"
  11 #include "AliAnalysisEtCuts.h"
  12
  13 #include "AliStack.h"
  14 #include "AliMCEvent.h"
  15 #include "AliESDEvent.h"
  16 #include "AliESDtrackCuts.h"
  17 #include "AliESDpid.h"
  18 #include "AliPID.h"
  19 #include "AliESDtrack.h"
  20 #include "AliVParticle.h"
  21 #include "AliAnalysisTask.h"
  22 #include "AliAnalysisManager.h"
  23 #include "AliAnalysisHadEtReconstructed.h"
  24 #include "AliAnalysisEtCuts.h"
  25 #include <iostream>
  26 #include "TRandom.h"
  27 #include "AliAnalysisEtCommon.h"
  28 #include "AliCentrality.h"
  29 #include "AliLog.h"
  30
  31 using namespace std;
  32
  33 ClassImp(AliAnalysisHadEtMonteCarlo);
  34
  35
  36 Int_t AliAnalysisHadEtMonteCarlo::fgNumSmearWidths = 4;
  37 Float_t AliAnalysisHadEtMonteCarlo::fgSmearWidths[4] = {0.005,0.006,0.007,0.008};
  38
  39 AliAnalysisHadEtMonteCarlo::AliAnalysisHadEtMonteCarlo():AliAnalysisHadEt()
  40                                                         ,fSimPiKPEt(0)
  41                                                         ,fSimHadEt(0)
  42                                                         ,fSimTotEt(0)
  43                                                         ,fInvestigateSmearing(0)
  44                                                         ,fInvestigateFull(0)
  45                                                         ,fInvestigateEMCal(0)
  46                                                         ,fInvestigatePHOS(0)
  47                                                         ,fInvestigatePiKP(0)
  48                                                         ,fRequireITSHits(0)
  49                                                         ,fBaryonEnhancement(0)
  50                                                         ,fUseRecoPt(0)
  51                                                         ,fPtSmearer(0)
  52 {
  53 }
  54 AliAnalysisHadEtMonteCarlo::~AliAnalysisHadEtMonteCarlo(){//destructor
  55   if(fPtSmearer) delete fPtSmearer;
  56 }
  57
  58 void AliAnalysisHadEtMonteCarlo::ResetEventValues(){//resetting event variables
  59   AliAnalysisHadEt::ResetEventValues();
  60     fSimHadEt=0.0;
  61     fSimTotEt=0.0;
  62     fSimPiKPEt=0.0;
  63 }
  64 Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
  65 { // analyse MC and real event info
  66   FillHisto1D("NEvents",0.5,1);
  67   if(!ev || !ev2){
  68     AliFatal("ERROR: Event does not exist");
  69     return 0;
  70   }
  71   AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*>(ev);
  72   AliESDEvent *realEvent = dynamic_cast<AliESDEvent*>(ev2);
  73   if(!mcEvent || !realEvent){
  74     AliFatal("ERROR: mcEvent or realEvent does not exist");
  75     return 0;
  76   }
  77   AliStack *stack = mcEvent->Stack();
  78   fCentBin= -1;
  79   fGoodEvent = kTRUE;//for p+p collisions if we made it this far we have a good event
  80   if(fDataSet==20100){//If this is Pb+Pb
  81     AliCentrality *centrality = realEvent->GetCentrality();
  82     if(fNCentBins<21) fCentBin= centrality->GetCentralityClass10(fCentralityMethod);
  83     else{ fCentBin= centrality->GetCentralityClass5(fCentralityMethod);}
  84     if(fCentBin ==-1) fGoodEvent = kFALSE;//but for Pb+Pb events we don't want to count events where we did not find a centrality
  85   }
  86   AnalyseEvent(ev);
  87
  88   //for PID
  89   AliESDpid *pID = new AliESDpid();//This is identified as a memory leak in valgrind but I delete this object so I think it may be a problem with AliESDpid.
  90
  91   //=============================================
  92
  93   //Roughly following $ALICE_ROOT/PWG0/dNdEta/AlidNdEtaCorrectionTask
  94
  95   //=============================================TPC&&ITS=============================================
  96   TString *strTPC = new TString("TPC");
  97   TString *strITS = new TString("ITS");
  98   TString *strTPCITS = new TString("TPCITS");
  99   Int_t lastcutset = 1;
 100   if(fRequireITSHits) lastcutset = 2;
 101   for(Int_t cutset=0;cutset<=lastcutset;cutset++){
 102     TString *cutName = NULL;
 103     TObjArray* list = NULL;
 104     switch(cutset){
 105     case 0:
 106       cutName = strTPC;
 107       list = fEsdtrackCutsTPC->GetAcceptedTracks(realEvent);
 108       break;
 109     case 1:
 110       cutName = strITS;
 111       list = fEsdtrackCutsITS->GetAcceptedTracks(realEvent);
 112       break;
 113     case 2:
 114       cutName = strTPCITS;
 115       list = fEsdtrackCutsITSTPC->GetAcceptedTracks(realEvent);
 116       break;
 117     default:
 118       cerr<<"Error:  cannot fill histograms!"<<endl;
 119       return -1;
 120     }
 121     Int_t nGoodTracks = list->GetEntries();
 122     for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++)
 123       {
 124         AliESDtrack *track = dynamic_cast<AliESDtrack*> (list->At(iTrack));
 125         if (!track)
 126           {
 127             Printf("ERROR: Could not get track %d", iTrack);
 128             continue;
 129           }
 130         else{
 131           Float_t nSigmaPion,nSigmaProton,nSigmaKaon,nSigmaElectron;
 132           pID->MakeTPCPID(track);
 133           pID->MakeITSPID(track);
 134           if(cutset!=1){
 135             nSigmaPion = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kPion));
 136             nSigmaProton = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kProton));
 137             nSigmaKaon = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kKaon));
 138             nSigmaElectron = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kElectron));
 139           }
 140           else{
 141             nSigmaPion = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kPion));
 142             nSigmaProton = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kProton));
 143             nSigmaKaon = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kKaon));
 144             nSigmaElectron = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kElectron));
 145           }
 146 //        bool isPion = (nSigmaPion<3.0 && nSigmaProton>2.0 && nSigmaKaon>2.0);
 147 //        bool isElectron = (nSigmaElectron<2.0 && nSigmaPion>4.0 && nSigmaProton>3.0 && nSigmaKaon>3.0);
 148 //        bool isKaon = (nSigmaPion>3.0 && nSigmaProton>2.0 && nSigmaKaon<2.0);
 149 //        bool isProton = (nSigmaPion>3.0 && nSigmaProton<2.0 && nSigmaKaon>2.0);
 150           bool isPion = (nSigmaPion<3.0 && nSigmaProton>2.0 && nSigmaKaon>2.0);
 151           bool isElectron = (nSigmaElectron<2.0 && nSigmaPion>4.0 && nSigmaProton>3.0 && nSigmaKaon>3.0);
 152           bool isKaon = (nSigmaPion>3.0 && nSigmaProton>3.0 && nSigmaKaon<3.0 && track->Pt()<0.45);
 153           bool isProton = (nSigmaPion>3.0 && nSigmaProton<3.0 && nSigmaKaon>3.0 && track->Pt()<0.9);
 154
 155           bool unidentified = (!isProton && !isKaon && !isElectron);
 156           Float_t dEdx = track->GetTPCsignal();
 157           if(cutset==1) dEdx = track->GetITSsignal();
 158
 159           FillHisto2D(Form("dEdxAll%s",cutName->Data()),track->P(),dEdx,1.0);
 160
 161           UInt_t label = (UInt_t)TMath::Abs(track->GetLabel());
 162           TParticle  *simPart  = stack->Particle(label);
 163           if(!simPart) {
 164             Printf("no MC particle\n");
 165             continue;
 166           }
 167           else{//analysis
 168             if(stack->IsPhysicalPrimary(label)){
 169               if (TMath::Abs(simPart->Eta()) < fCuts->GetCommonEtaCut()){
 170                 Int_t pdgCode =  simPart->GetPDG(0)->PdgCode();
 171                 Int_t mypid = 0;
 172                 if(pdgCode==AliAnalysisHadEt::fgPiPlusCode) mypid = 1;
 173                 if(pdgCode==fgProtonCode) mypid = 2;
 174                 if(pdgCode==fgKPlusCode) mypid = 3;
 175                 if(pdgCode==fgEPlusCode) mypid = 4;
 176                 if(pdgCode==fgPiMinusCode) mypid = 1;
 177                 if(pdgCode==fgAntiProtonCode) mypid = 2;
 178                 if(pdgCode==fgKMinusCode) mypid = 3;
 179                 if(pdgCode==fgEMinusCode) mypid = 4;
 180                 bool filled = false;
 181                 //============Charged hadrons===================================
 182                 //identified...
 183                 if(isPion){
 184                   if(pdgCode!=fgPiPlusCode && pdgCode!=fgPiMinusCode){
 185                     FillHisto2D(Form("MisidentifiedPIDs%s",cutName->Data()),1,mypid,1);
 186                     //if(mypid==0)cerr<<"I was misidentified! I'm not a pion! I am a "<<simPart->GetName()<<endl;
 187                   }
 188                   float myEt = Et(simPart);
 189                   if(track->Charge()>0){ FillHisto2D(Form("EtReconstructed%sIdentifiedPiPlus",cutName->Data()),track->Pt(),track->Eta(),myEt);}
 190                   else{ FillHisto2D(Form("EtReconstructed%sIdentifiedPiMinus",cutName->Data()),track->Pt(),track->Eta(),myEt);}
 191                   FillHisto2D(Form("dEdxPion%s",cutName->Data()),track->P(),dEdx,1.0);
 192                 }
 193                 if(isProton){
 194                   if(pdgCode!=fgProtonCode && pdgCode!=fgAntiProtonCode){
 195                     FillHisto2D(Form("MisidentifiedPIDs%s",cutName->Data()),2,mypid,1);
 196                   }
 197                   float myEt = Et(simPart);
 198                   if(track->Charge()>0){ FillHisto2D(Form("EtReconstructed%sIdentifiedProton",cutName->Data()),track->Pt(),track->Eta(),myEt);}
 199                   else{ FillHisto2D(Form("EtReconstructed%sIdentifiedAntiProton",cutName->Data()),track->Pt(),track->Eta(),myEt);}
 200                   if(fBaryonEnhancement){
 201                     myEt = myEt*ProtonBaryonEnhancement(track->Pt());
 202                     if(track->Charge()>0){ FillHisto2D(Form("EtReconstructed%sIdentifiedProtonEnhanced",cutName->Data()),track->Pt(),track->Eta(),myEt);}
 203                     else{ FillHisto2D(Form("EtReconstructed%sIdentifiedAntiProtonEnhanced",cutName->Data()),track->Pt(),track->Eta(),myEt);}
 204                   }
 205                   FillHisto2D(Form("dEdxProton%s",cutName->Data()),track->P(),dEdx,1.0);
 206                 }
 207                 if(isKaon){
 208                   if(pdgCode!=fgKMinusCode && pdgCode!=fgKPlusCode){
 209                     FillHisto2D(Form("MisidentifiedPIDs%s",cutName->Data()),3,mypid,1);
 210                   }
 211                   float myEt = Et(simPart);
 212                   if(track->Charge()>0){ FillHisto2D(Form("EtReconstructed%sIdentifiedKPlus",cutName->Data()),track->Pt(),track->Eta(),myEt);}
 213                   else{ FillHisto2D(Form("EtReconstructed%sIdentifiedKMinus",cutName->Data()),track->Pt(),track->Eta(),myEt);}
 214                   FillHisto2D(Form("dEdxKaon%s",cutName->Data()),track->P(),dEdx,1.0);
 215                 }
 216                 if(isElectron){
 217                   if(pdgCode!=fgEMinusCode && pdgCode!=fgEPlusCode){
 218                     FillHisto2D(Form("MisidentifiedPIDs%s",cutName->Data()),4,mypid,1);
 219                   }
 220                   float myEt = Et(simPart);
 221                   if(track->Charge()>0){ FillHisto2D(Form("EtReconstructed%sIdentifiedEPlus",cutName->Data()),track->Pt(),track->Eta(),myEt);}
 222                   else{ FillHisto2D(Form("EtReconstructed%sIdentifiedEMinus",cutName->Data()),track->Pt(),track->Eta(),myEt);}
 223                   FillHisto2D(Form("dEdxElectron%s",cutName->Data()),track->P(),dEdx,1.0);
 224                 }
 225                 if(unidentified){
 226                   if(pdgCode!=fgEMinusCode && pdgCode!=fgEPlusCode){
 227                     float myEtPi = Et(simPart,fgPionMass);
 228                     float myEtP = Et(simPart,fgProtonMass);
 229                     float myEtK = Et(simPart,fgKaonMass);
 230                     float myEt = Et(simPart);
 231                     FillHisto2D(Form("EtReconstructed%sUnidentifiedAssumingPion",cutName->Data()),track->Pt(),track->Eta(),myEtPi);
 232                     FillHisto2D(Form("EtReconstructed%sUnidentifiedAssumingProton",cutName->Data()),track->Pt(),track->Eta(),myEtP);
 233                     FillHisto2D(Form("EtReconstructed%sUnidentifiedAssumingKaon",cutName->Data()),track->Pt(),track->Eta(),myEtK);
 234                     FillHisto2D(Form("EtReconstructed%sUnidentified",cutName->Data()),track->Pt(),track->Eta(),myEt);
 235                     FillHisto2D(Form("EtNReconstructed%sUnidentified",cutName->Data()),track->Pt(),track->Eta(),1.0);
 236                     if(pdgCode == fgPiPlusCode||pdgCode == fgPiMinusCode){
 237                       FillHisto2D(Form("EtReconstructed%sUnidentifiedPionAssumingPion",cutName->Data()),track->Pt(),track->Eta(),myEtPi);
 238                       FillHisto2D(Form("EtReconstructed%sUnidentifiedPionAssumingProton",cutName->Data()),track->Pt(),track->Eta(),myEtP);
 239                       FillHisto2D(Form("EtReconstructed%sUnidentifiedPionAssumingKaon",cutName->Data()),track->Pt(),track->Eta(),myEtK);
 240                       FillHisto2D(Form("EtReconstructed%sUnidentifiedPion",cutName->Data()),track->Pt(),track->Eta(),myEt);
 241                       FillHisto2D(Form("EtNReconstructed%sUnidentifiedPion",cutName->Data()),track->Pt(),track->Eta(),1.0);
 242                     }
 243                     if(pdgCode == fgKPlusCode||pdgCode == fgKMinusCode){
 244                       FillHisto2D(Form("EtReconstructed%sUnidentifiedKaonAssumingPion",cutName->Data()),track->Pt(),track->Eta(),myEtPi);
 245                       FillHisto2D(Form("EtReconstructed%sUnidentifiedKaonAssumingProton",cutName->Data()),track->Pt(),track->Eta(),myEtP);
 246                       FillHisto2D(Form("EtReconstructed%sUnidentifiedKaonAssumingKaon",cutName->Data()),track->Pt(),track->Eta(),myEtK);
 247                       FillHisto2D(Form("EtReconstructed%sUnidentifiedKaon",cutName->Data()),track->Pt(),track->Eta(),myEt);
 248                       FillHisto2D(Form("EtNReconstructed%sUnidentifiedKaon",cutName->Data()),track->Pt(),track->Eta(),1.0);
 249                     }
 250                     if(pdgCode == fgProtonCode||pdgCode == fgAntiProtonCode){
 251                       FillHisto2D(Form("EtReconstructed%sUnidentifiedProtonAssumingPion",cutName->Data()),track->Pt(),track->Eta(),myEtPi);
 252                       FillHisto2D(Form("EtReconstructed%sUnidentifiedProtonAssumingProton",cutName->Data()),track->Pt(),track->Eta(),myEtP);
 253                       FillHisto2D(Form("EtReconstructed%sUnidentifiedProtonAssumingKaon",cutName->Data()),track->Pt(),track->Eta(),myEtK);
 254                       FillHisto2D(Form("EtReconstructed%sUnidentifiedProton",cutName->Data()),track->Pt(),track->Eta(),myEt);
 255                       FillHisto2D(Form("EtNReconstructed%sUnidentifiedProton",cutName->Data()),track->Pt(),track->Eta(),1.0);
 256                       if(fBaryonEnhancement){
 257                         myEt = myEt*ProtonBaryonEnhancement(track->Pt());
 258                         FillHisto2D(Form("EtReconstructed%sUnidentifiedProtonAssumingPionEnhanced",cutName->Data()),track->Pt(),track->Eta(),myEtPi);
 259                         FillHisto2D(Form("EtReconstructed%sUnidentifiedProtonEnhanced",cutName->Data()),track->Pt(),track->Eta(),myEt);
 260                         FillHisto2D(Form("EtNReconstructed%sUnidentifiedProtonEnhanced",cutName->Data()),track->Pt(),track->Eta(),1.0);
 261                       }
 262                     }
 263                   }
 264                   FillHisto2D(Form("dEdxUnidentified%s",cutName->Data()),track->P(),dEdx,1.0);
 265                   FillHisto1D(Form("UnidentifiedPIDs%s",cutName->Data()),mypid,1);
 266                 }
 267                 //...simulated
 268                 float myEtSim = Et(simPart);
 269                 float myEtReco = 0.0;
 270                 if(pdgCode == fgPiPlusCode){
 271                   float myEt = Et(simPart);
 272                   float myEtP = Et(simPart,fgProtonMass);
 273                   float myEtK = Et(simPart,fgKaonMass);
 274                   myEtReco = Et(track->P(),track->Theta(),fgPiPlusCode,track->Charge());
 275                   float pT = simPart->Pt();
 276                   float eta = simPart->Eta();
 277                   if(fUseRecoPt){//Then we switch the pT and the Et
 278                     myEt = myEtReco;
 279                     pT = track->Pt();
 280                     eta = track->Eta();
 281                   }
 282                   FillHisto2D(Form("EtReconstructed%sPiPlus",cutName->Data()),pT,eta,myEt);
 283                   FillHisto2D(Form("EtReconstructed%sChargedHadron",cutName->Data()),pT,eta,myEt);
 284                   FillHisto2D(Form("EtNReconstructed%sPiPlus",cutName->Data()),pT,eta,myEt);
 285                   FillHisto2D(Form("EtNReconstructed%sChargedHadron",cutName->Data()),pT,eta,myEt);
 286                   if(fCentBin>=0){//if a centrality bin was defined
 287                     FillHisto2D(Form("EtNReconstructed%sPiPlusCB%i",cutName->Data(),fCentBin),pT,eta,myEt);
 288                     FillHisto2D(Form("EtNReconstructed%sChargedHadronCB%i",cutName->Data(),fCentBin),pT,eta,myEt);
 289                   }
 290                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingPion",cutName->Data()),pT,eta,myEt);
 291                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingProton",cutName->Data()),pT,eta,myEtP);
 292                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingKaon",cutName->Data()),pT,eta,myEtK);
 293                   FillHisto2D(Form("EtReconstructed%sPiPlusAssumingKaon",cutName->Data()),pT,eta,myEtK);
 294                   FillHisto2D(Form("EtReconstructed%sPiPlusAssumingProton",cutName->Data()),pT,eta,myEtP);
 295                   filled = true;
 296                 }
 297                 if(pdgCode == fgPiMinusCode){
 298                   float myEt = Et(simPart);
 299                   float myEtP = Et(simPart,fgProtonMass);
 300                   float myEtK = Et(simPart,fgKaonMass);
 301                   myEtReco = Et(track->P(),track->Theta(),fgPiMinusCode,track->Charge());
 302                   float pT = simPart->Pt();
 303                   float eta = simPart->Eta();
 304                   if(fUseRecoPt){//Then we switch the pT and the Et
 305                     myEt = myEtReco;
 306                     pT = track->Pt();
 307                     eta = track->Eta();
 308                   }
 309                   FillHisto2D(Form("EtReconstructed%sPiMinus",cutName->Data()),pT,eta,myEt);
 310                   FillHisto2D(Form("EtReconstructed%sChargedHadron",cutName->Data()),pT,eta,myEt);
 311                   FillHisto2D(Form("EtNReconstructed%sPiMinus",cutName->Data()),pT,eta,myEt);
 312                   FillHisto2D(Form("EtNReconstructed%sChargedHadron",cutName->Data()),pT,eta,myEt);
 313                   if(fCentBin>=0){//if a centrality bin was defined
 314                     FillHisto2D(Form("EtNReconstructed%sPiMinusCB%i",cutName->Data(),fCentBin),pT,eta,myEt);
 315                     FillHisto2D(Form("EtNReconstructed%sChargedHadronCB%i",cutName->Data(),fCentBin),pT,eta,myEt);
 316                   }
 317                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingPion",cutName->Data()),pT,eta,myEt);
 318                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingProton",cutName->Data()),pT,eta,myEtP);
 319                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingKaon",cutName->Data()),pT,eta,myEtK);
 320                   FillHisto2D(Form("EtReconstructed%sPiMinusAssumingKaon",cutName->Data()),pT,eta,myEtK);
 321                   FillHisto2D(Form("EtReconstructed%sPiMinusAssumingProton",cutName->Data()),pT,eta,myEtP);
 322                   filled = true;
 323                 }
 324                 if(pdgCode == fgKPlusCode){
 325                   float myEt = Et(simPart);
 326                   float myEtPi = Et(simPart,fgPionMass);
 327                   float myEtP = Et(simPart,fgProtonMass);
 328                   myEtReco = Et(track->P(),track->Theta(),fgKPlusCode,track->Charge());
 329                   float pT = simPart->Pt();
 330                   float eta = simPart->Eta();
 331                   if(fUseRecoPt){//Then we switch the pT and the Et
 332                     myEt = myEtReco;
 333                     pT = track->Pt();
 334                     eta = track->Eta();
 335                   }
 336                   FillHisto2D(Form("EtReconstructed%sKPlus",cutName->Data()),pT,eta,myEt);
 337                   FillHisto2D(Form("EtReconstructed%sChargedHadron",cutName->Data()),pT,eta,myEt);
 338                   FillHisto2D(Form("EtNReconstructed%sKPlus",cutName->Data()),pT,eta,myEt);
 339                   FillHisto2D(Form("EtNReconstructed%sChargedHadron",cutName->Data()),pT,eta,myEt);
 340                   if(fCentBin>=0){//if a centrality bin was defined
 341                     FillHisto2D(Form("EtNReconstructed%sKPlusCB%i",cutName->Data(),fCentBin),pT,eta,myEt);
 342                     FillHisto2D(Form("EtNReconstructed%sChargedHadronCB%i",cutName->Data(),fCentBin),pT,eta,myEt);
 343                   }
 344                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingPion",cutName->Data()),pT,eta,myEtPi);
 345                   FillHisto2D(Form("EtReconstructed%sKPlusAssumingPion",cutName->Data()),pT,eta,myEtPi);
 346                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingKaon",cutName->Data()),pT,eta,myEt);
 347                   FillHisto2D(Form("EtReconstructed%sKPlusAssumingKaon",cutName->Data()),pT,eta,myEt);
 348                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingProton",cutName->Data()),pT,eta,myEtP);
 349                   FillHisto2D(Form("EtReconstructed%sKPlusAssumingProton",cutName->Data()),pT,eta,myEtP);
 350                   filled = true;
 351                 }
 352                 if(pdgCode == fgKMinusCode){
 353                   float myEt = Et(simPart);
 354                   float myEtPi = Et(simPart,fgPionMass);
 355                   float myEtP = Et(simPart,fgProtonMass);
 356                   myEtReco = Et(track->P(),track->Theta(),fgKMinusCode,track->Charge());
 357                   float pT = simPart->Pt();
 358                   float eta = simPart->Eta();
 359                   if(fUseRecoPt){//Then we switch the pT and the Et
 360                     myEt = myEtReco;
 361                     pT = track->Pt();
 362                     eta = track->Eta();
 363                   }
 364                   FillHisto2D(Form("EtReconstructed%sKMinus",cutName->Data()),pT,eta,myEt);
 365                   FillHisto2D(Form("EtReconstructed%sChargedHadron",cutName->Data()),pT,eta,myEt);
 366                   FillHisto2D(Form("EtNReconstructed%sKMinus",cutName->Data()),pT,eta,myEt);
 367                   FillHisto2D(Form("EtNReconstructed%sChargedHadron",cutName->Data()),pT,eta,myEt);
 368                   if(fCentBin>=0){//if a centrality bin was defined
 369                     FillHisto2D(Form("EtNReconstructed%sKMinusCB%i",cutName->Data(),fCentBin),pT,eta,myEt);
 370                     FillHisto2D(Form("EtNReconstructed%sChargedHadronCB%i",cutName->Data(),fCentBin),pT,eta,myEt);
 371                   }
 372                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingPion",cutName->Data()),pT,eta,myEtPi);
 373                   FillHisto2D(Form("EtReconstructed%sKMinusAssumingPion",cutName->Data()),pT,eta,myEtPi);
 374                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingKaon",cutName->Data()),pT,eta,myEt);
 375                   FillHisto2D(Form("EtReconstructed%sKMinusAssumingKaon",cutName->Data()),pT,eta,myEt);
 376                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingProton",cutName->Data()),pT,eta,myEtP);
 377                   FillHisto2D(Form("EtReconstructed%sKMinusAssumingProton",cutName->Data()),pT,eta,myEtP);
 378                   filled = true;
 379                 }
 380                 if(pdgCode == fgProtonCode){
 381                   float myEt = Et(simPart);
 382                   float myEtPi = Et(simPart,fgPionMass);
 383                   float myEtK = Et(simPart,fgKaonMass);
 384                   myEtReco = Et(track->P(),track->Theta(),fgProtonCode,track->Charge());
 385                   float pT = simPart->Pt();
 386                   float eta = simPart->Eta();
 387                   if(fUseRecoPt){//Then we switch the pT and the Et
 388                     myEt = myEtReco;
 389                     pT = track->Pt();
 390                     eta = track->Eta();
 391                   }
 392                   FillHisto2D(Form("EtReconstructed%sProton",cutName->Data()),pT,eta,myEt);
 393                   FillHisto2D(Form("EtReconstructed%sChargedHadron",cutName->Data()),pT,eta,myEt);
 394                   FillHisto2D(Form("EtNReconstructed%sProton",cutName->Data()),pT,eta,myEt);
 395                   FillHisto2D(Form("EtNReconstructed%sChargedHadron",cutName->Data()),pT,eta,myEt);
 396                   if(fCentBin>=0){//if a centrality bin was defined
 397                     FillHisto2D(Form("EtNReconstructed%sProtonCB%i",cutName->Data(),fCentBin),pT,eta,myEt);
 398                     FillHisto2D(Form("EtNReconstructed%sChargedHadronCB%i",cutName->Data(),fCentBin),pT,eta,myEt);
 399                   }
 400                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingPion",cutName->Data()),pT,eta,myEtPi);
 401                   FillHisto2D(Form("EtReconstructed%sProtonAssumingPion",cutName->Data()),pT,eta,myEtPi);
 402                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingKaon",cutName->Data()),pT,eta,myEtK);
 403                   FillHisto2D(Form("EtReconstructed%sProtonAssumingKaon",cutName->Data()),pT,eta,myEtK);
 404                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingProton",cutName->Data()),pT,eta,myEt);
 405                   FillHisto2D(Form("EtReconstructed%sProtonAssumingProton",cutName->Data()),pT,eta,myEt);
 406                   filled = true;
 407
 408                   if(fBaryonEnhancement){
 409                     float enhancement = ProtonBaryonEnhancement(track->Pt());
 410                     FillHisto2D(Form("EtReconstructed%sProtonEnhanced",cutName->Data()),pT,eta,myEt*enhancement);
 411                     FillHisto2D(Form("EtNReconstructed%sProtonEnhanced",cutName->Data()),pT,eta,myEt*enhancement);
 412                     FillHisto2D(Form("EtReconstructed%sProtonAssumingPionEnhanced",cutName->Data()),pT,eta,myEtPi*enhancement);
 413                   }
 414
 415                 }
 416                 if(pdgCode == fgAntiProtonCode){
 417                   float myEt = Et(simPart);
 418                   float myEtPi = Et(simPart,fgPionMass);
 419                   float myEtK = Et(simPart,fgKaonMass);
 420                   myEtReco = Et(track->P(),track->Theta(),fgAntiProtonCode,track->Charge());
 421                   float pT = simPart->Pt();
 422                   float eta = simPart->Eta();
 423                   if(fUseRecoPt){//Then we switch the pT and the Et
 424                     myEt = myEtReco;
 425                     pT = track->Pt();
 426                     eta = track->Eta();
 427                   }
 428                   FillHisto2D(Form("EtReconstructed%sAntiProton",cutName->Data()),pT,eta,myEt);
 429                   FillHisto2D(Form("EtReconstructed%sChargedHadron",cutName->Data()),pT,eta,myEt);
 430                   FillHisto2D(Form("EtNReconstructed%sAntiProton",cutName->Data()),pT,eta,myEt);
 431                   FillHisto2D(Form("EtNReconstructed%sChargedHadron",cutName->Data()),pT,eta,myEt);
 432                   if(fCentBin>=0){//if a centrality bin was defined
 433                     FillHisto2D(Form("EtNReconstructed%sAntiProtonCB%i",cutName->Data(),fCentBin),pT,eta,myEt);
 434                     FillHisto2D(Form("EtNReconstructed%sChargedHadronCB%i",cutName->Data(),fCentBin),pT,eta,myEt);
 435                   }
 436                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingPion",cutName->Data()),pT,eta,myEtPi);
 437                   FillHisto2D(Form("EtReconstructed%sAntiProtonAssumingPion",cutName->Data()),pT,eta,myEtPi);
 438                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingKaon",cutName->Data()),pT,eta,myEtK);
 439                   FillHisto2D(Form("EtReconstructed%sAntiProtonAssumingKaon",cutName->Data()),pT,eta,myEtK);
 440                   FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingProton",cutName->Data()),pT,eta,myEt);
 441                   FillHisto2D(Form("EtReconstructed%sAntiProtonAssumingProton",cutName->Data()),pT,eta,myEt);
 442                   filled = true;
 443                   if(fBaryonEnhancement){
 444                         float enhancement = ProtonBaryonEnhancement(track->Pt());
 445                         FillHisto2D(Form("EtReconstructed%sAntiProtonEnhanced",cutName->Data()),pT,eta,myEt*enhancement);
 446                     FillHisto2D(Form("EtNReconstructed%sAntiProtonEnhanced",cutName->Data()),pT,eta,myEt*enhancement);
 447                     FillHisto2D(Form("EtReconstructed%sAntiProtonAssumingPionEnhanced",cutName->Data()),pT,eta,myEtPi*enhancement);
 448                   }
 449                 }
 450                 if(pdgCode == fgEPlusCode){
 451                   float myEt = Et(simPart);
 452                   FillHisto2D(Form("EtReconstructed%sEPlus",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
 453                   if(!isElectron || unidentified){
 454                     float myEtPi = Et(simPart,fgPionMass);
 455                     FillHisto2D(Form("EtReconstructed%sMisidentifiedElectrons",cutName->Data()),track->Pt(),track->Eta(),myEtPi);
 456                   }
 457                   filled = true;
 458                 }
 459                 if(pdgCode == fgEMinusCode){
 460                   if(!isElectron || unidentified){
 461                     float myEtPi = Et(simPart,fgPionMass);
 462                     FillHisto2D(Form("EtReconstructed%sMisidentifiedElectrons",cutName->Data()),track->Pt(),track->Eta(),myEtPi);
 463                   }
 464                   float myEt = Et(simPart);
 465                     FillHisto2D(Form("EtReconstructed%sEMinus",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
 466                     filled = true;
 467                 }
 468                 if(myEtReco>0.0){FillHisto2D(Form("ETresolution%s",cutName->Data()),myEtReco,(myEtSim-myEtReco)/myEtReco,1.0);}
 469                 if(track->Pt()>0.0){FillHisto2D(Form("pTresolution%s",cutName->Data()),track->Pt(),(simPart->Pt() - track->Pt())/track->Pt(),1.0);}
 470                 if(track->P()>0.0){FillHisto2D(Form("presolution%s",cutName->Data()),track->P(),(simPart->P() - track->P())/track->P(),1.0);}
 471                 FillHisto1D(Form("pTsim%s",cutName->Data()),simPart->Pt(),1.0);
 472                 FillHisto1D(Form("pTrec%s",cutName->Data()),track->Pt(),1.0);
 473                 if(fCentBin!=-1){
 474                   FillHisto1D(Form("pTsim%sCB%i",cutName->Data(),fCentBin),simPart->Pt(),1.0);
 475                   FillHisto1D(Form("pTrec%sCB%i",cutName->Data(),fCentBin),track->Pt(),1.0);
 476                 }
 477
 478               }
 479
 480             }
 481             else{//not a primary - we're after V0 daughters!
 482               if (TMath::Abs(simPart->Eta()) < fCuts->GetCommonEtaCut()){
 483                 TParticle *mom = stack->Particle(simPart->GetFirstMother());
 484                 if(mom){
 485                   TParticlePDG *pc = mom->GetPDG(0);
 486                   if(pc){
 487                     Int_t pdgCode =  mom->GetPDG(0)->PdgCode();
 488                     if(pdgCode == fgLambdaCode){
 489                       float myEt = Et(simPart);
 490                       float pT = simPart->Pt();
 491                       float eta = simPart->Eta();
 492                       if(fUseRecoPt){//Then we switch the pT and the Et
 493                         myEt = Et(track->P(),track->Theta(),simPart->GetPDG(0)->PdgCode(),track->Charge());
 494                         pT = track->Pt();
 495                         eta = track->Eta();
 496                       }
 497                       FillHisto2D(Form("EtReconstructed%sLambdaDaughters",cutName->Data()),pT,eta,myEt);
 498                       Float_t weight = LambdaWeight(mom->Pt());
 499                       if(fBaryonEnhancement){
 500                         float enhancement = ProtonBaryonEnhancement(track->Pt());
 501                         weight = weight*enhancement;
 502                       }
 503                       FillHisto2D(Form("EtReconstructed%sLambdaDaughtersReweighted",cutName->Data()),pT,eta,myEt*weight);
 504                     }
 505                     if(pdgCode == fgAntiLambdaCode){
 506                       float myEt = Et(simPart);
 507                       float pT = simPart->Pt();
 508                       float eta = simPart->Eta();
 509                       if(fUseRecoPt){//Then we switch the pT and the Et
 510                         myEt = Et(track->P(),track->Theta(),simPart->GetPDG(0)->PdgCode(),track->Charge());
 511                         pT = track->Pt();
 512                         eta = track->Eta();
 513                       }
 514                       FillHisto2D(Form("EtReconstructed%sAntiLambdaDaughters",cutName->Data()),pT,eta,myEt);
 515                       Float_t weight = AntiLambdaWeight(mom->Pt());
 516                       if(fBaryonEnhancement){
 517                         float enhancement = ProtonBaryonEnhancement(track->Pt());
 518                         weight = weight*enhancement;
 519                       }
 520                       FillHisto2D(Form("EtReconstructed%sAntiLambdaDaughtersReweighted",cutName->Data()),pT,eta,myEt*weight);
 521                     }
 522                     if(pdgCode == fgK0SCode){
 523                       float myEt = Et(simPart);
 524                       float pT = simPart->Pt();
 525                       float eta = simPart->Eta();
 526                       if(fUseRecoPt){//Then we switch the pT and the Et
 527                         myEt = Et(track->P(),track->Theta(),simPart->GetPDG(0)->PdgCode(),track->Charge());
 528                         pT = track->Pt();
 529                         eta = track->Eta();
 530                       }
 531                       FillHisto2D(Form("EtReconstructed%sK0SDaughters",cutName->Data()),pT,eta,myEt);
 532                       Float_t weight = K0Weight(mom->Pt());
 533                       FillHisto2D(Form("EtReconstructed%sK0SDaughtersReweighted",cutName->Data()),pT,eta,myEt*weight);
 534                     }
 535                     if(pdgCode == fgXiCode){
 536                       float myEt = Et(simPart);
 537                       float pT = simPart->Pt();
 538                       float eta = simPart->Eta();
 539                       if(fUseRecoPt){//Then we switch the pT and the Et
 540                         myEt = Et(track->P(),track->Theta(),simPart->GetPDG(0)->PdgCode(),track->Charge());
 541                         pT = track->Pt();
 542                         eta = track->Eta();
 543                       }
 544                       FillHisto2D(Form("EtReconstructed%sXiDaughters",cutName->Data()),pT,eta,myEt);
 545                     }
 546                     if(pdgCode == fgAntiXiCode){
 547                       float myEt = Et(simPart);
 548                       float pT = simPart->Pt();
 549                       float eta = simPart->Eta();
 550                       if(fUseRecoPt){//Then we switch the pT and the Et
 551                         myEt = Et(track->P(),track->Theta(),simPart->GetPDG(0)->PdgCode(),track->Charge());
 552                         pT = track->Pt();
 553                         eta = track->Eta();
 554                       }
 555                       FillHisto2D(Form("EtReconstructed%sAntiXiDaughters",cutName->Data()),pT,eta,myEt);
 556                     }
 557                     if(pdgCode == fgOmegaCode){
 558                       float myEt = Et(simPart);
 559                       float pT = simPart->Pt();
 560                       float eta = simPart->Eta();
 561                       if(fUseRecoPt){//Then we switch the pT and the Et
 562                         myEt = Et(track->P(),track->Theta(),simPart->GetPDG(0)->PdgCode(),track->Charge());
 563                         pT = track->Pt();
 564                         eta = track->Eta();
 565                       }
 566                       FillHisto2D(Form("EtReconstructed%sOmegaDaughters",cutName->Data()),pT,eta,myEt);
 567                     }
 568                     if(pdgCode == fgXiCode){
 569                       float myEt = Et(simPart);
 570                       float pT = simPart->Pt();
 571                       float eta = simPart->Eta();
 572                       if(fUseRecoPt){//Then we switch the pT and the Et
 573                         myEt = Et(track->P(),track->Theta(),simPart->GetPDG(0)->PdgCode(),track->Charge());
 574                         pT = track->Pt();
 575                         eta = track->Eta();
 576                       }
 577                       FillHisto2D(Form("EtReconstructed%sAntiOmegaDaughters",cutName->Data()),pT,eta,myEt);
 578                     }
 579
 580                     if(mom->GetFirstMother()>0){
 581                       TParticle *grandma = stack->Particle(mom->GetFirstMother());
 582                       if(grandma){
 583                         Int_t pdgCodeMom =  mom->GetPDG(0)->PdgCode();
 584                         if(pdgCodeMom==fgPiPlusCode || pdgCodeMom==fgPiMinusCode || pdgCodeMom==fgProtonCode ||pdgCodeMom==fgAntiProtonCode || pdgCodeMom==fgKPlusCode || pdgCode==fgKMinusCode){
 585                           Int_t pdgCodeGrandma =  grandma->GetPDG(0)->PdgCode();
 586
 587                           if(pdgCodeGrandma == fgXiCode){
 588                             float myEt = Et(simPart);
 589                             float pT = simPart->Pt();
 590                             float eta = simPart->Eta();
 591                             if(fUseRecoPt){//Then we switch the pT and the Et
 592                               myEt = Et(track->P(),track->Theta(),simPart->GetPDG(0)->PdgCode(),track->Charge());
 593                               pT = track->Pt();
 594                               eta = track->Eta();
 595                             }
 596                             FillHisto2D(Form("EtReconstructed%sXiDaughters",cutName->Data()),pT,eta,myEt);
 597                           }
 598                           if(pdgCodeGrandma == fgAntiXiCode){
 599                             float myEt = Et(simPart);
 600                             float pT = simPart->Pt();
 601                             float eta = simPart->Eta();
 602                             if(fUseRecoPt){//Then we switch the pT and the Et
 603                               myEt = Et(track->P(),track->Theta(),simPart->GetPDG(0)->PdgCode(),track->Charge());
 604                               pT = track->Pt();
 605                               eta = track->Eta();
 606                             }
 607                             FillHisto2D(Form("EtReconstructed%sAntiXiDaughters",cutName->Data()),pT,eta,myEt);
 608                           }
 609                           if(pdgCodeGrandma == fgOmegaCode){
 610                             float myEt = Et(simPart);
 611                             float pT = simPart->Pt();
 612                             float eta = simPart->Eta();
 613                             if(fUseRecoPt){//Then we switch the pT and the Et
 614                               myEt = Et(track->P(),track->Theta(),simPart->GetPDG(0)->PdgCode(),track->Charge());
 615                               pT = track->Pt();
 616                               eta = track->Eta();
 617                             }
 618                             FillHisto2D(Form("EtReconstructed%sOmegaDaughters",cutName->Data()),pT,eta,myEt);
 619                           }
 620                           if(pdgCodeGrandma == fgXiCode){
 621                             float myEt = Et(simPart);
 622                             float pT = simPart->Pt();
 623                             float eta = simPart->Eta();
 624                             if(fUseRecoPt){//Then we switch the pT and the Et
 625                               myEt = Et(track->P(),track->Theta(),simPart->GetPDG(0)->PdgCode(),track->Charge());
 626                               pT = track->Pt();
 627                               eta = track->Eta();
 628                             }
 629                             FillHisto2D(Form("EtReconstructed%sAntiOmegaDaughters",cutName->Data()),pT,eta,myEt);
 630                           }
 631
 632                         }
 633                       }
 634                     }
 635                   }
 636                 }
 637               }
 638             }
 639           }
 640
 641         }
 642       }
 643     delete list;
 644   }
 645   delete pID;
 646   delete strTPC;
 647   delete strITS;
 648   delete strTPCITS;
 649   //delete AliESDpid;
 650   return 1;
 651 }
 652 Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
 653 { // analyse MC event
 654      ResetEventValues();
 655      if(!ev){
 656             AliFatal("ERROR: Event does not exist");
 657             return 0;
 658      }
 659
 660     // Get us an mc event
 661     AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*>(ev);
 662     if(!mcEvent){
 663       AliFatal("ERROR: MC Event does not exist");
 664       return 0;
 665     }
 666     // Let's play with the stack!
 667     AliStack *stack = mcEvent->Stack();
 668
 669     Int_t nPrim = stack->GetNtrack();
 670
 671     Float_t fSimPiKPEtPtSmeared = 0;
 672     Float_t fSimPiKPEtEfficiencySmeared = 0;
 673     Float_t fSimPiKPEtPtCutSmearedTPC = 0;
 674     Float_t fSimPiKPEtPtCutSmearedITS = 0;
 675     Float_t fSimPiKPEtPIDSmeared = 0;
 676     Float_t fSimPiKPEtPIDSmearedNoID = 0;
 677     //=================Tracks which may or may not have been reconstructed=================
 678
 679     for (Int_t iPart = 0; iPart < nPrim; iPart++)
 680     {
 681
 682       TParticle *part = stack->Particle(iPart);//This line is identified as a loss of memory by valgrind, however, the pointer still belongs to the stack, so it's the stack's problem
 683
 684         if (!part)
 685           {
 686             Printf("ERROR: Could not get particle %d", iPart);
 687             continue;
 688           }
 689
 690         //TParticlePDG *pc = part->GetPDG(0);
 691
 692         // Check if it is a primary particle
 693         if (stack->IsPhysicalPrimary(iPart)){//primaries
 694
 695           if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())           {
 696
 697             Int_t pdgCode =  part->GetPDG(0)->PdgCode();
 698             bool filled = false;
 699             //Investigating smearing...
 700             //Numbers are realistic correction factors from previous studies
 701             if(fInvestigateSmearing){
 702               if(pdgCode==fgPiPlusCode ||pdgCode==fgPiMinusCode ||pdgCode==fgKPlusCode ||pdgCode==fgKMinusCode ||pdgCode==fgProtonCode ||pdgCode==fgAntiProtonCode){
 703                 //To investigate Smearing...
 704                 Float_t myet = Et(part);
 705                 fSimPiKPEt += myet;
 706                 Float_t theta = part->Theta();
 707                 Short_t charge = 1;
 708                 Float_t momentum = part->P();
 709                 //pt smearing
 710                 Float_t pSmeared = momentum *  fPtSmearer->Gaus(1,0.005);//Gaussian centered around 1
 711                 fSimPiKPEtPtSmeared += Et(pSmeared,theta,pdgCode,charge);
 712                 //Efficiency smearing
 713                 //to mock up the difference between TPC only tracks in p+p (~90% efficiency at high pT) and TPC+ITS tracks in Pb+Pb (about 70% efficiency at high pT) a factor of 7/9 was added in front
 714                 float efficiency = 7.0/9.0*2.26545*TMath::Exp(-TMath::Power(9.99977e-01/part->Pt(),7.85488e-02));//simple rough efficiency from fitting curve
 715                 if(fPtSmearer->Binomial(1,efficiency) ==1){
 716                   fSimPiKPEtEfficiencySmeared += (1.0/efficiency)*myet;
 717                 }
 718                 //pT cut smeared
 719                 if(part->Pt()>0.10){fSimPiKPEtPtCutSmearedITS +=1.00988*myet;}
 720                 if(part->Pt()>0.15){fSimPiKPEtPtCutSmearedTPC +=1.02994*myet;}
 721                 //PID smearing
 722                 fSimPiKPEtPIDSmearedNoID += 1.03018015790601458*Et(momentum,theta,fgPiPlusCode,charge);
 723                 if(part->P()<1.0){//then the particle would have been ID'd
 724                   fSimPiKPEtPIDSmeared += 1.00918051514628582*myet;
 725                 }
 726                 else{//Then it would have been assumed to be a pion
 727                   fSimPiKPEtPIDSmeared += 1.00918051514628582*Et(momentum,theta,fgPiPlusCode,charge);
 728                 }
 729               }
 730             }
 731
 732             //============Charged hadrons===================================
 733             if(pdgCode == fgPiPlusCode){
 734               float myEt = Et(part);
 735               float myEtP = Et(part,fgProtonMass);
 736               float myEtK = Et(part,fgKaonMass);
 737
 738               fSimHadEt += myEt;
 739               fSimTotEt += myEt;
 740
 741               FillHisto2D("EtSimulatedPiPlus",part->Pt(),part->Eta(),myEt);
 742               FillHisto2D("EtNSimulatedPiPlus",part->Pt(),part->Eta(),1.0);
 743               FillHisto2D("EtSimulatedChargedHadron",part->Pt(),part->Eta(),myEt);
 744               FillHisto2D("EtNSimulatedChargedHadron",part->Pt(),part->Eta(),1.0);
 745               if(fCentBin>=0){//if a centrality bin was defined
 746                 FillHisto2D(Form("EtNSimulatedPiPlusCB%i",fCentBin),part->Pt(),part->Eta(),1.0);
 747                 FillHisto2D(Form("EtNSimulatedChargedHadronCB%i",fCentBin),part->Pt(),part->Eta(),1.0);
 748               }
 749               FillHisto2D("EtSimulatedChargedHadronAssumingPion",part->Pt(),part->Eta(),myEt);
 750               FillHisto2D("EtSimulatedChargedHadronAssumingProton",part->Pt(),part->Eta(),myEtP);
 751               FillHisto2D("EtSimulatedPiPlusAssumingProton",part->Pt(),part->Eta(),myEtP);
 752               FillHisto2D("EtSimulatedChargedHadronAssumingKaon",part->Pt(),part->Eta(),myEtK);
 753               FillHisto2D("EtSimulatedPiPlusAssumingKaon",part->Pt(),part->Eta(),myEtK);
 754               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
 755               Short_t charge = 1;
 756               Float_t myEtLow = Et(0.0,part->Theta(),pdgCode,charge);
 757               Float_t myEtITS = Et(0.10/TMath::Sin(part->Theta()),part->Theta(),pdgCode,charge);
 758               Float_t myEtTPC = Et(0.15/TMath::Sin(part->Theta()),part->Theta(),pdgCode,charge);
 759               FillHisto2D("EtSimulatedChargedHadronAssumingNoPt",part->Pt(),part->Eta(),myEtLow);
 760               FillHisto2D("EtSimulatedChargedHadronAssumingPtTPCCut",part->Pt(),part->Eta(),myEtTPC);
 761               FillHisto2D("EtSimulatedChargedHadronAssumingPtITSCut",part->Pt(),part->Eta(),myEtITS);
 762               filled = true;
 763             }
 764             if(pdgCode == fgPiMinusCode){
 765               float myEt = Et(part);
 766               float myEtP = Et(part,fgProtonMass);
 767               float myEtK = Et(part,fgKaonMass);
 768               fSimHadEt += myEt;
 769               fSimTotEt += myEt;
 770               FillHisto2D("EtSimulatedPiMinus",part->Pt(),part->Eta(),myEt);
 771               FillHisto2D("EtNSimulatedPiMinus",part->Pt(),part->Eta(),1.0);
 772               FillHisto2D("EtSimulatedChargedHadron",part->Pt(),part->Eta(),myEt);
 773               FillHisto2D("EtNSimulatedChargedHadron",part->Pt(),part->Eta(),1.0);
 774               if(fCentBin>=0){//if a centrality bin was defined
 775                 FillHisto2D(Form("EtNSimulatedPiMinusCB%i",fCentBin),part->Pt(),part->Eta(),1.0);
 776                 FillHisto2D(Form("EtNSimulatedChargedHadronCB%i",fCentBin),part->Pt(),part->Eta(),1.0);
 777               }
 778               FillHisto2D("EtSimulatedChargedHadronAssumingPion",part->Pt(),part->Eta(),myEt);
 779               FillHisto2D("EtSimulatedChargedHadronAssumingProton",part->Pt(),part->Eta(),myEtP);
 780               FillHisto2D("EtSimulatedPiMinusAssumingProton",part->Pt(),part->Eta(),myEtP);
 781               FillHisto2D("EtSimulatedChargedHadronAssumingKaon",part->Pt(),part->Eta(),myEtK);
 782               FillHisto2D("EtSimulatedPiMinusAssumingKaon",part->Pt(),part->Eta(),myEtK);
 783               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
 784               Short_t charge = -1;
 785               Float_t myEtLow = Et(0.0,part->Theta(),pdgCode,charge);
 786               Float_t myEtITS = Et(0.10/TMath::Sin(part->Theta()),part->Theta(),pdgCode,charge);
 787               Float_t myEtTPC = Et(0.15/TMath::Sin(part->Theta()),part->Theta(),pdgCode,charge);
 788               FillHisto2D("EtSimulatedChargedHadronAssumingNoPt",part->Pt(),part->Eta(),myEtLow);
 789               FillHisto2D("EtSimulatedChargedHadronAssumingPtTPCCut",part->Pt(),part->Eta(),myEtTPC);
 790               FillHisto2D("EtSimulatedChargedHadronAssumingPtITSCut",part->Pt(),part->Eta(),myEtITS);
 791               filled = true;
 792             }
 793             if(pdgCode == fgKPlusCode){
 794               float myEt = Et(part);
 795               float myEtPi = Et(part,fgPionMass);
 796               float myEtP = Et(part,fgProtonMass);
 797               fSimHadEt += myEt;
 798               fSimTotEt += myEt;
 799               FillHisto2D("EtSimulatedKPlus",part->Pt(),part->Eta(),myEt);
 800               FillHisto2D("EtNSimulatedKPlus",part->Pt(),part->Eta(),1.0);
 801               FillHisto2D("EtSimulatedChargedHadron",part->Pt(),part->Eta(),myEt);
 802               FillHisto2D("EtNSimulatedChargedHadron",part->Pt(),part->Eta(),1.0);
 803               if(fCentBin>=0){//if a centrality bin was defined
 804                 FillHisto2D(Form("EtNSimulatedKPlusCB%i",fCentBin),part->Pt(),part->Eta(),1.0);
 805                 FillHisto2D(Form("EtNSimulatedChargedHadronCB%i",fCentBin),part->Pt(),part->Eta(),1.0);
 806               }
 807               FillHisto2D("EtSimulatedChargedHadronAssumingPion",part->Pt(),part->Eta(),myEtPi);
 808               FillHisto2D("EtSimulatedKPlusAssumingPion",part->Pt(),part->Eta(),myEtPi);
 809               FillHisto2D("EtSimulatedChargedHadronAssumingProton",part->Pt(),part->Eta(),myEtP);
 810               FillHisto2D("EtSimulatedKPlusAssumingProton",part->Pt(),part->Eta(),myEtP);
 811               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
 812               Short_t charge = 1;
 813               Float_t myEtLow = Et(0.0,part->Theta(),pdgCode,charge);
 814               Float_t myEtITS = Et(0.10/TMath::Sin(part->Theta()),part->Theta(),pdgCode,charge);
 815               Float_t myEtTPC = Et(0.15/TMath::Sin(part->Theta()),part->Theta(),pdgCode,charge);
 816               FillHisto2D("EtSimulatedChargedHadronAssumingNoPt",part->Pt(),part->Eta(),myEtLow);
 817               FillHisto2D("EtSimulatedChargedHadronAssumingPtTPCCut",part->Pt(),part->Eta(),myEtTPC);
 818               FillHisto2D("EtSimulatedChargedHadronAssumingPtITSCut",part->Pt(),part->Eta(),myEtITS);
 819               filled = true;
 820             }
 821             if(pdgCode == fgKMinusCode){
 822               float myEt = Et(part);
 823               float myEtPi = Et(part,fgPionMass);
 824               float myEtP = Et(part,fgProtonMass);
 825               fSimHadEt += myEt;
 826               fSimTotEt += myEt;
 827               FillHisto2D("EtSimulatedKMinus",part->Pt(),part->Eta(),myEt);
 828               FillHisto2D("EtNSimulatedKMinus",part->Pt(),part->Eta(),1.0);
 829               FillHisto2D("EtSimulatedChargedHadron",part->Pt(),part->Eta(),myEt);
 830               FillHisto2D("EtNSimulatedChargedHadron",part->Pt(),part->Eta(),1.0);
 831               if(fCentBin>=0){//if a centrality bin was defined
 832                 FillHisto2D(Form("EtNSimulatedKMinusCB%i",fCentBin),part->Pt(),part->Eta(),1.0);
 833                 FillHisto2D(Form("EtNSimulatedChargedHadronCB%i",fCentBin),part->Pt(),part->Eta(),1.0);
 834               }
 835               FillHisto2D("EtSimulatedChargedHadronAssumingPion",part->Pt(),part->Eta(),myEtPi);
 836               FillHisto2D("EtSimulatedKMinusAssumingPion",part->Pt(),part->Eta(),myEtPi);
 837               FillHisto2D("EtSimulatedChargedHadronAssumingProton",part->Pt(),part->Eta(),myEtP);
 838               FillHisto2D("EtSimulatedKMinusAssumingProton",part->Pt(),part->Eta(),myEtP);
 839               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
 840               Short_t charge = -1;
 841               Float_t myEtLow = Et(0.0,part->Theta(),pdgCode,charge);
 842               Float_t myEtITS = Et(0.10/TMath::Sin(part->Theta()),part->Theta(),pdgCode,charge);
 843               Float_t myEtTPC = Et(0.15/TMath::Sin(part->Theta()),part->Theta(),pdgCode,charge);
 844               FillHisto2D("EtSimulatedChargedHadronAssumingNoPt",part->Pt(),part->Eta(),myEtLow);
 845               FillHisto2D("EtSimulatedChargedHadronAssumingPtTPCCut",part->Pt(),part->Eta(),myEtTPC);
 846               FillHisto2D("EtSimulatedChargedHadronAssumingPtITSCut",part->Pt(),part->Eta(),myEtITS);
 847               filled = true;
 848             }
 849             if(pdgCode == fgProtonCode){
 850               float myEt = Et(part);
 851               float myEtPi = Et(part,fgPionMass);
 852               float myEtK = Et(part,fgKaonMass);
 853               fSimHadEt += myEt;
 854               fSimTotEt += myEt;
 855               FillHisto2D("EtSimulatedProton",part->Pt(),part->Eta(),myEt);
 856               FillHisto2D("EtNSimulatedProton",part->Pt(),part->Eta(),1.0);
 857               FillHisto2D("EtSimulatedChargedHadron",part->Pt(),part->Eta(),myEt);
 858               FillHisto2D("EtNSimulatedChargedHadron",part->Pt(),part->Eta(),1.0);
 859               if(fCentBin>=0){//if a centrality bin was defined
 860                 FillHisto2D(Form("EtNSimulatedProtonCB%i",fCentBin),part->Pt(),part->Eta(),1.0);
 861                 FillHisto2D(Form("EtNSimulatedChargedHadronCB%i",fCentBin),part->Pt(),part->Eta(),1.0);
 862               }
 863               FillHisto2D("EtSimulatedChargedHadronAssumingPion",part->Pt(),part->Eta(),myEtPi);
 864               FillHisto2D("EtSimulatedProtonAssumingPion",part->Pt(),part->Eta(),myEtPi);
 865               FillHisto2D("EtSimulatedChargedHadronAssumingKaon",part->Pt(),part->Eta(),myEtK);
 866               FillHisto2D("EtSimulatedProtonAssumingKaon",part->Pt(),part->Eta(),myEtK);
 867               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
 868               Short_t charge = 1;
 869               Float_t myEtLow = Et(0.0,part->Theta(),pdgCode,charge);
 870               Float_t myEtITS = Et(0.10/TMath::Sin(part->Theta()),part->Theta(),pdgCode,charge);
 871               Float_t myEtTPC = Et(0.15/TMath::Sin(part->Theta()),part->Theta(),pdgCode,charge);
 872               FillHisto2D("EtSimulatedChargedHadronAssumingNoPt",part->Pt(),part->Eta(),myEtLow);
 873               FillHisto2D("EtSimulatedChargedHadronAssumingPtTPCCut",part->Pt(),part->Eta(),myEtTPC);
 874               FillHisto2D("EtSimulatedChargedHadronAssumingPtITSCut",part->Pt(),part->Eta(),myEtITS);
 875               filled = true;
 876               if(fBaryonEnhancement){
 877                 float enhancement = ProtonBaryonEnhancement(part->Pt());
 878                 FillHisto2D("EtSimulatedProtonEnhanced",part->Pt(),part->Eta(),myEt*enhancement);
 879                 FillHisto2D("EtNSimulatedProtonEnhanced",part->Pt(),part->Eta(),1.0*enhancement);
 880                 FillHisto2D("EtSimulatedProtonAssumingPionEnhanced",part->Pt(),part->Eta(),myEtPi*enhancement);
 881               }
 882             }
 883             if(pdgCode == fgAntiProtonCode){
 884               float myEt = Et(part);
 885               float myEtPi = Et(part,fgPionMass);
 886               float myEtK = Et(part,fgKaonMass);
 887               fSimHadEt += myEt;
 888               fSimTotEt += myEt;
 889               FillHisto2D("EtSimulatedAntiProton",part->Pt(),part->Eta(),myEt);
 890               FillHisto2D("EtNSimulatedAntiProton",part->Pt(),part->Eta(),1.0);
 891               FillHisto2D("EtSimulatedChargedHadron",part->Pt(),part->Eta(),myEt);
 892               FillHisto2D("EtNSimulatedChargedHadron",part->Pt(),part->Eta(),1.0);
 893               if(fCentBin>=0){//if a centrality bin was defined
 894                 FillHisto2D(Form("EtNSimulatedAntiProtonCB%i",fCentBin),part->Pt(),part->Eta(),1.0);
 895                 FillHisto2D(Form("EtNSimulatedChargedHadronCB%i",fCentBin),part->Pt(),part->Eta(),1.0);
 896               }
 897               FillHisto2D("EtSimulatedChargedHadronAssumingPion",part->Pt(),part->Eta(),myEtPi);
 898               FillHisto2D("EtSimulatedAntiProtonAssumingPion",part->Pt(),part->Eta(),myEtPi);
 899               FillHisto2D("EtSimulatedChargedHadronAssumingKaon",part->Pt(),part->Eta(),myEtK);
 900               FillHisto2D("EtSimulatedAntiProtonAssumingKaon",part->Pt(),part->Eta(),myEtK);
 901               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
 902               Short_t charge = -1;
 903               Float_t myEtLow = Et(0.0,part->Theta(),pdgCode,charge);
 904               Float_t myEtITS = Et(0.10/TMath::Sin(part->Theta()),part->Theta(),pdgCode,charge);
 905               Float_t myEtTPC = Et(0.15/TMath::Sin(part->Theta()),part->Theta(),pdgCode,charge);
 906               FillHisto2D("EtSimulatedChargedHadronAssumingNoPt",part->Pt(),part->Eta(),myEtLow);
 907               FillHisto2D("EtSimulatedChargedHadronAssumingPtTPCCut",part->Pt(),part->Eta(),myEtTPC);
 908               FillHisto2D("EtSimulatedChargedHadronAssumingPtITSCut",part->Pt(),part->Eta(),myEtITS);
 909               filled = true;
 910               if(fBaryonEnhancement){
 911                 float enhancement = ProtonBaryonEnhancement(part->Pt());
 912                 FillHisto2D("EtSimulatedAntiProtonEnhanced",part->Pt(),part->Eta(),myEt*enhancement);
 913                 FillHisto2D("EtNSimulatedAntiProtonEnhanced",part->Pt(),part->Eta(),1.0*enhancement);
 914                 FillHisto2D("EtSimulatedAntiProtonAssumingPionEnhanced",part->Pt(),part->Eta(),myEtPi*enhancement);
 915               }
 916             }
 917             //============Other hadrons===================================
 918
 919             if(pdgCode == fgNeutronCode){
 920               float myEt = Et(part);
 921               fSimHadEt += myEt;
 922               fSimTotEt += myEt;
 923               FillHisto2D("EtSimulatedNeutron",part->Pt(),part->Eta(),myEt);
 924               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
 925               filled = true;
 926             }
 927             if(pdgCode == fgAntiNeutronCode){
 928               float myEt = Et(part);
 929               fSimHadEt += myEt;
 930               fSimTotEt += myEt;
 931               FillHisto2D("EtSimulatedAntiNeutron",part->Pt(),part->Eta(),myEt);
 932               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
 933               filled = true;
 934             }
 935             if(pdgCode == fgLambdaCode){
 936               float myEt = Et(part);
 937               fSimHadEt += myEt;
 938               fSimTotEt += myEt;
 939               //cout<<"I am a simulated lambda! pt "<<part->Pt()<<" eta "<<part->Eta()<<endl;
 940               FillHisto2D("EtSimulatedLambda",part->Pt(),part->Eta(),myEt);
 941               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
 942               Float_t weight = LambdaWeight(part->Pt());
 943               if(fBaryonEnhancement){
 944                 float enhancement = ProtonBaryonEnhancement(part->Pt());
 945                 weight = weight*enhancement;
 946               }
 947               FillHisto2D("EtSimulatedLambdaReweighted",part->Pt(),part->Eta(),myEt*weight);
 948               Int_t ndaughters = part->GetNDaughters();
 949               for(Int_t idaughter = 0;idaughter<ndaughters;idaughter++){
 950                 Int_t daughterindex = part->GetDaughter(idaughter);
 951                 if(daughterindex<0 || daughterindex>1e5) continue;
 952                 TParticle *daughter = stack->ParticleFromTreeK(daughterindex);
 953                 if(daughter){
 954                   if(daughter->GetPDG(0)){
 955                     Int_t daughtercode = daughter->GetPDG(0)->PdgCode();
 956                     if(daughtercode==fgPiMinusCode || daughtercode==fgProtonCode){
 957                       myEt = Et(daughter);
 958                       FillHisto2D("EtSimulatedLambdaDaughters",daughter->Pt(),daughter->Eta(),myEt);
 959                       FillHisto2D("EtSimulatedLambdaDaughtersReweighted",daughter->Pt(),daughter->Eta(),myEt*weight);
 960                       //cout<<"Lambda daughter is a "<<daughter->GetName()<<endl;
 961                     }
 962                   }
 963                   else{
 964                     //cout<<"Lambda daughter is a "<<daughter->GetName()<<endl;
 965                   }
 966                 }
 967               }
 968               filled = true;
 969             }
 970             if(pdgCode == fgAntiLambdaCode){
 971               float myEt = Et(part);
 972               fSimHadEt += myEt;
 973               fSimTotEt += myEt;
 974               FillHisto2D("EtSimulatedAntiLambda",part->Pt(),part->Eta(),myEt);
 975               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
 976               Float_t weight = AntiLambdaWeight(part->Pt());
 977               if(fBaryonEnhancement){
 978                 float enhancement = ProtonBaryonEnhancement(part->Pt());
 979                 weight = weight*enhancement;
 980               }
 981               FillHisto2D("EtSimulatedAntiLambdaReweighted",part->Pt(),part->Eta(),myEt*weight);
 982               Int_t ndaughters = part->GetNDaughters();
 983               for(Int_t idaughter = 0;idaughter<ndaughters;idaughter++){
 984                 Int_t daughterindex = part->GetDaughter(idaughter);
 985                 if(daughterindex<0 || daughterindex>1e5) continue;
 986                 TParticle *daughter = stack->ParticleFromTreeK(daughterindex);
 987                 if(daughter){
 988                   if(daughter->GetPDG(0)){
 989                     Int_t daughtercode = daughter->GetPDG(0)->PdgCode();
 990                     if(daughtercode==fgPiPlusCode || daughtercode==fgAntiProtonCode){
 991                       myEt = Et(daughter);
 992                       FillHisto2D("EtSimulatedAntiLambdaDaughters",daughter->Pt(),daughter->Eta(),myEt);
 993                       FillHisto2D("EtSimulatedAntiLambdaDaughtersReweighted",daughter->Pt(),daughter->Eta(),myEt*weight);
 994                       //cout<<"AntiLambda daughter is a "<<daughter->GetName()<<endl;
 995                     }
 996                   }
 997                   else{
 998                     //cout<<"AntiLambda daughter is a "<<daughter->GetName()<<endl;
 999                   }
1000                 }
1001               }
1002               filled = true;
1003             }
1004             if(pdgCode == fgK0SCode){
1005               float myEt = Et(part);
1006               fSimHadEt += myEt;
1007               fSimTotEt += myEt;
1008               FillHisto2D("EtSimulatedK0S",part->Pt(),part->Eta(),myEt);
1009               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
1010               Float_t weight = K0Weight(part->Pt());
1011               FillHisto2D("EtSimulatedK0SReweighted",part->Pt(),part->Eta(),myEt*weight);
1012               Int_t ndaughters = part->GetNDaughters();
1013               for(Int_t idaughter = 0;idaughter<ndaughters;idaughter++){
1014                 Int_t daughterindex = part->GetDaughter(idaughter);
1015                 if(daughterindex<0 || daughterindex>1e5) continue;
1016                 TParticle *daughter = stack->ParticleFromTreeK(daughterindex);
1017                 if(daughter){
1018                   if(daughter->GetPDG(0)){
1019
1020                     Int_t daughtercode = daughter->GetPDG(0)->PdgCode();
1021                     if(daughtercode==fgPiMinusCode || daughtercode==fgPiPlusCode){
1022                       myEt = Et(daughter);
1023                       FillHisto2D("EtSimulatedK0SDaughters",daughter->Pt(),daughter->Eta(),myEt);
1024                       FillHisto2D("EtSimulatedK0SDaughtersReweighted",daughter->Pt(),daughter->Eta(),myEt*weight);
1025                       //cout<<"K0S daughter is a "<<daughter->GetName()<<endl;
1026                     }
1027                   }
1028                   else{
1029                     //cout<<"K0S daughter is a "<<daughter->GetName()<<endl;
1030                   }
1031                 }
1032               }
1033               filled = true;
1034             }
1035             if(pdgCode == fgK0LCode){
1036               float myEt = Et(part);
1037               fSimHadEt += myEt;
1038               fSimTotEt += myEt;
1039               FillHisto2D("EtSimulatedK0L",part->Pt(),part->Eta(),myEt);
1040               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
1041               Float_t weight = K0Weight(part->Pt());
1042               FillHisto2D("EtSimulatedK0LReweighted",part->Pt(),part->Eta(),myEt*weight);
1043               filled = true;
1044             }
1045             if(pdgCode == fgOmegaCode){
1046               float myEt = Et(part);
1047               fSimHadEt += myEt;
1048               fSimTotEt += myEt;
1049               FillHisto2D("EtSimulatedOmega",part->Pt(),part->Eta(),myEt);
1050               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
1051               Int_t ndaughters = part->GetNDaughters();
1052               for(Int_t idaughter = 0;idaughter<ndaughters;idaughter++){
1053                 Int_t daughterindex = part->GetDaughter(idaughter);
1054                 if(daughterindex<0 || daughterindex>1e5) continue;
1055                 TParticle *daughter = stack->Particle(daughterindex);
1056                 if(daughter){
1057                   if(daughter->GetPDG(0)){
1058
1059                     Int_t daughtercode = daughter->GetPDG(0)->PdgCode();
1060                     if(daughtercode==fgPiPlusCode || daughtercode==fgProtonCode || daughtercode==fgKMinusCode){
1061                       myEt = Et(daughter);
1062                       FillHisto2D("EtSimulatedOmegaDaughters",daughter->Pt(),daughter->Eta(),myEt);
1063                     //cout<<"Omega daughter is a "<<daughter->GetName()<<endl;
1064                     }
1065                   }
1066                   else{
1067                     //cout<<"Omega daughter is a "<<daughter->GetName()<<endl;
1068                   }
1069                 }
1070               }
1071               filled = true;
1072             }
1073             if(pdgCode == fgAntiOmegaCode){
1074               float myEt = Et(part);
1075               fSimHadEt += myEt;
1076               fSimTotEt += myEt;
1077               FillHisto2D("EtSimulatedOmega",part->Pt(),part->Eta(),myEt);
1078               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
1079               Int_t ndaughters = part->GetNDaughters();
1080               for(Int_t idaughter = 0;idaughter<ndaughters;idaughter++){
1081                 Int_t daughterindex = part->GetDaughter(idaughter);
1082                 if(daughterindex<0 || daughterindex>1e5) continue;
1083                 TParticle *daughter = stack->ParticleFromTreeK(daughterindex);
1084                 if(daughter){
1085                   if(daughter->GetPDG(0)){
1086                     Int_t daughtercode = daughter->GetPDG(0)->PdgCode();
1087                     if(daughtercode==fgPiMinusCode || daughtercode==fgAntiProtonCode || daughtercode==fgKPlusCode){
1088                       myEt = Et(daughter);
1089                       FillHisto2D("EtSimulatedAntiOmegaDaughters",daughter->Pt(),daughter->Eta(),myEt);
1090                       //cout<<"AntiOmega daughter is a "<<daughter->GetName()<<endl;
1091                     }
1092                   }
1093                   else{
1094                     //cout<<"AntiOmega daughter is a "<<daughter->GetName()<<endl;
1095                   }
1096                 }
1097               }
1098               filled = true;
1099             }
1100             //There are two codes for Sigmas
1101             if(pdgCode == fgSigmaCode || pdgCode == -3222){
1102               float myEt = Et(part);
1103               fSimHadEt += myEt;
1104               fSimTotEt += myEt;
1105               FillHisto2D("EtSimulatedSigma",part->Pt(),part->Eta(),myEt);
1106               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
1107               filled = true;
1108             }
1109             if(pdgCode == fgAntiSigmaCode || pdgCode == 3222){
1110               float myEt = Et(part);
1111               fSimHadEt += myEt;
1112               fSimTotEt += myEt;
1113               FillHisto2D("EtSimulatedAntiSigma",part->Pt(),part->Eta(),myEt);
1114               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
1115               filled = true;
1116             }
1117             if(pdgCode == fgXiCode){
1118               float myEt = Et(part);
1119               fSimHadEt += myEt;
1120               fSimTotEt += myEt;
1121               FillHisto2D("EtSimulatedXi",part->Pt(),part->Eta(),myEt);
1122               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
1123               Int_t ndaughters = part->GetNDaughters();
1124               for(Int_t idaughter = 0;idaughter<ndaughters;idaughter++){
1125                 Int_t daughterindex = part->GetDaughter(idaughter);
1126                 if(daughterindex<0 || daughterindex>1e5 || daughterindex>1e5) continue;
1127                 //cerr<<"Daughter index "<<daughterindex<<" npart "<<nPrim<<endl;
1128                 TParticle *daughter = stack->ParticleFromTreeK(daughterindex);
1129                 if(daughter){
1130                   if(daughter->GetPDG(0)){
1131
1132                     Int_t daughtercode = daughter->GetPDG(0)->PdgCode();
1133                     if(daughtercode==fgPiPlusCode || daughtercode==fgProtonCode || daughtercode==fgPiMinusCode){
1134                       myEt = Et(daughter);
1135                       FillHisto2D("EtSimulatedXiDaughters",daughter->Pt(),daughter->Eta(),myEt);
1136                     //cout<<"Xi daughter is a "<<daughter->GetName()<<endl;
1137                     }
1138                   }
1139                   else{
1140                     //cout<<"Xi daughter is a "<<daughter->GetName()<<endl;
1141                   }
1142                 }
1143               }
1144               filled = true;
1145             }
1146             if(pdgCode == fgAntiXiCode){
1147               float myEt = Et(part);
1148               fSimHadEt += myEt;
1149               fSimTotEt += myEt;
1150               FillHisto2D("EtSimulatedAntiXi",part->Pt(),part->Eta(),myEt);
1151               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
1152               Int_t ndaughters = part->GetNDaughters();
1153               for(Int_t idaughter = 0;idaughter<ndaughters;idaughter++){
1154                 Int_t daughterindex = part->GetDaughter(idaughter);
1155                 if(daughterindex<0 || daughterindex>1e5) continue;
1156                 TParticle *daughter = stack->ParticleFromTreeK(daughterindex);
1157                 if(daughter){
1158                   if(daughter->GetPDG(0)){
1159                     Int_t daughtercode = daughter->GetPDG(0)->PdgCode();
1160                     if(daughtercode==fgPiPlusCode || daughtercode==fgAntiProtonCode || daughtercode==fgPiMinusCode){
1161                       myEt = Et(daughter);
1162                       FillHisto2D("EtSimulatedAntiXiDaughters",daughter->Pt(),daughter->Eta(),myEt);
1163                       //cout<<"AntiXi daughter is a "<<daughter->GetName()<<endl;
1164                     }
1165                   }
1166                   else{
1167                     //cout<<"AntiXi daughter is a "<<daughter->GetName()<<endl;
1168                   }
1169                 }
1170               }
1171               filled = true;
1172             }
1173             if(pdgCode == fgXi0Code){
1174               float myEt = Et(part);
1175               fSimHadEt += myEt;
1176               fSimTotEt += myEt;
1177               FillHisto2D("EtSimulatedXi0",part->Pt(),part->Eta(),myEt);
1178               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
1179               filled = true;
1180             }
1181             if(pdgCode == fgAntiXi0Code){
1182               float myEt = Et(part);
1183               fSimHadEt += myEt;
1184               fSimTotEt += myEt;
1185               FillHisto2D("EtSimulatedAntiXi0",part->Pt(),part->Eta(),myEt);
1186               FillHisto2D("EtSimulatedAllHadron",part->Pt(),part->Eta(),myEt);
1187               filled = true;
1188             }
1189             //============electrons===================================
1190
1191             if(pdgCode == fgEPlusCode){
1192               float myEt = Et(part);
1193               fSimTotEt += myEt;
1194               FillHisto2D("EtSimulatedEPlus",part->Pt(),part->Eta(),myEt);
1195               filled = true;
1196             }
1197             if(pdgCode == fgEMinusCode){
1198               float myEt = Et(part);
1199               fSimTotEt += myEt;
1200               FillHisto2D("EtSimulatedEMinus",part->Pt(),part->Eta(),myEt);
1201               filled = true;
1202             }
1203             //============neutrals===================================
1204             if(pdgCode == fgGammaCode){
1205               TParticle *mom = NULL;
1206               Int_t pdgCodeMom = -99999999;
1207               float momEta = -30;
1208               float mompT = -5;
1209               if(part->GetFirstMother()){
1210                 mom = stack->Particle(part->GetFirstMother());
1211                 pdgCodeMom =  mom->GetPDG(0)->PdgCode();
1212                 momEta = mom->Eta();
1213                 mompT = mom->Pt();
1214               }
1215               //cout<<"I am a gamma and my mom is "<<mom->GetName()<<endl;
1216               //We want to separate the gammas by pi0, eta, omega0 but we don't want to double count energy so we get the et from the gamma daughter
1217               if(pdgCodeMom == fgEtaCode){
1218                 float myEt = Et(part);
1219                 fSimTotEt += myEt;
1220                 FillHisto2D("EtSimulatedEta",mompT,momEta,myEt);
1221                 filled = true;
1222               }
1223               if(pdgCodeMom == fgPi0Code){
1224                 float myEt = Et(part);
1225                 fSimTotEt += myEt;
1226                 FillHisto2D("EtSimulatedPi0",mompT,momEta,myEt);
1227                 filled = true;
1228               }
1229               if(pdgCodeMom == fgOmega0Code){
1230                 float myEt = Et(part);
1231                 fSimTotEt += myEt;
1232                 FillHisto2D("EtSimulatedOmega0",mompT,momEta,myEt);
1233                 filled = true;
1234               }
1235               if(!filled){
1236                 float myEt = Et(part);
1237                 fSimTotEt += myEt;
1238                 FillHisto2D("EtSimulatedGamma",part->Pt(),part->Eta(),myEt);
1239                 filled = true;
1240               }
1241             }
1242             if(pdgCode == fgEtaCode){
1243               float myEt = Et(part);
1244               fSimTotEt += myEt;
1245               FillHisto2D("EtSimulatedEta",part->Pt(),part->Eta(),myEt);
1246               filled = true;
1247             }
1248             if(pdgCode == fgPi0Code){
1249               float myEt = Et(part);
1250               fSimTotEt += myEt;
1251               FillHisto2D("EtSimulatedPi0",part->Pt(),part->Eta(),myEt);
1252               filled = true;
1253             }
1254             if(pdgCode == fgOmega0Code){
1255               float myEt = Et(part);
1256               fSimTotEt += myEt;
1257               FillHisto2D("EtSimulatedOmega0",part->Pt(),part->Eta(),myEt);
1258               filled = true;
1259             }
1260             if(!filled){
1261               //if( strcmp(pc->ParticleClass(),"Baryon")==0 || strcmp(pc->ParticleClass(),"Meson")==0 ){
1262               //cout<<"Did not find a place for "<<part->GetName()<<" "<<pdgCode<<" which is a "<<pc->ParticleClass()<<endl;
1263                 //}
1264             }
1265           }
1266         }
1267     }
1268
1269     if(fSimTotEt>0.0)FillHisto1D("SimTotEt",fSimTotEt,1.0);
1270     if(fSimHadEt>0.0)FillHisto1D("SimHadEt",fSimHadEt,1.0);
1271     if(fSimPiKPEt>0.0)FillHisto1D("SimPiKPEt",fSimPiKPEt,1.0);
1272     if(fCentBin != -1){//if we have Pb+Pb and a centrality bin was found
1273       if(fSimTotEt>0.0) FillHisto1D(Form("SimTotEtCB%i",fCentBin),fSimTotEt,1.0);
1274       if(fSimHadEt>0.0) FillHisto1D(Form("SimHadEtCB%i",fCentBin),fSimHadEt,1.0);
1275       if(fSimPiKPEt>0.0)FillHisto1D(Form("SimPiKPEtCB%i",fCentBin),fSimPiKPEt,1.0);
1276     }
1277
1278     if(fInvestigateSmearing){
1279       //Smearing histograms
1280       if(fSimPiKPEt>0.0) FillHisto2D("SimPiKPEtMinusSimPtSmeared",fSimPiKPEt,(fSimPiKPEt-fSimPiKPEtPtSmeared)/fSimPiKPEt,1.0);
1281       FillHisto1D("SimPiKPEtPtSmeared",fSimPiKPEtPtSmeared,1.0);
1282       if(fSimPiKPEt>0.0) FillHisto2D("SimPiKPEtMinusSimEfficiencySmeared",fSimPiKPEt,(fSimPiKPEt-fSimPiKPEtEfficiencySmeared)/fSimPiKPEt,1.0);
1283       FillHisto1D("SimPiKPEtEfficiencySmeared",fSimPiKPEtEfficiencySmeared,1.0);
1284       if(fSimPiKPEt>0.0) FillHisto2D("SimPiKPEtMinusSimPtCutSmearedTPC",fSimPiKPEt,(fSimPiKPEt-fSimPiKPEtPtCutSmearedTPC)/fSimPiKPEt,1.0);
1285       FillHisto1D("SimPiKPEtPtCutSmearedTPC",fSimPiKPEtPtCutSmearedTPC,1.0);
1286       if(fSimPiKPEt>0.0) FillHisto2D("SimPiKPEtMinusSimPtCutSmearedITS",fSimPiKPEt,(fSimPiKPEt-fSimPiKPEtPtCutSmearedITS)/fSimPiKPEt,1.0);
1287       FillHisto1D("SimPiKPEtPtCutSmearedITS",fSimPiKPEtPtCutSmearedTPC,1.0);
1288       if(fSimPiKPEt>0.0) FillHisto2D("SimPiKPEtMinusSimPIDSmeared",fSimPiKPEt,(fSimPiKPEt-fSimPiKPEtPIDSmeared)/fSimPiKPEt,1.0);
1289       FillHisto1D("SimPiKPEtPIDSmeared",fSimPiKPEtPIDSmeared,1.0);
1290       if(fSimPiKPEt>0.0) FillHisto2D("SimPiKPEtMinusSimPIDSmearedNoID",fSimPiKPEt,(fSimPiKPEt-fSimPiKPEtPIDSmearedNoID)/fSimPiKPEt,1.0);
1291       FillHisto1D("SimPiKPEtPIDSmearedNoID",fSimPiKPEtPIDSmearedNoID,1.0);
1292     }
1293     return 1;
1294
1295 }
1296
1297 void AliAnalysisHadEtMonteCarlo::Init()
1298 { // Init
1299     AliAnalysisHadEt::Init();
1300     if(!fPtSmearer) fPtSmearer = new TRandom();
1301 }
1302 void AliAnalysisHadEtMonteCarlo::CreateHistograms(){
1303   //for simulated Et only (no reconstruction)
1304   CreateEtaPtHisto2D(TString("EtSimulatedPiPlus"),TString("Simulated E_{T} from #pi^{+}"));
1305   CreateEtaPtHisto2D("EtSimulatedPiMinus","Simulated E_{T} from #pi^{-}");
1306   CreateEtaPtHisto2D("EtSimulatedKPlus","Simulated E_{T} from K^{+}");
1307   CreateEtaPtHisto2D("EtSimulatedKMinus","Simulated E_{T} from K^{-}");
1308   CreateEtaPtHisto2D("EtSimulatedProton","Simulated E_{T} from p");
1309   CreateEtaPtHisto2D("EtSimulatedAntiProton","Simulated E_{T} from #bar{p}");//Both baryon enhancement and strangeness rescaling
1310   if(fBaryonEnhancement){
1311     CreateEtaPtHisto2D("EtSimulatedProtonEnhanced","Simulated E_{T} from p");
1312     CreateEtaPtHisto2D("EtSimulatedAntiProtonEnhanced","Simulated E_{T} from #bar{p}");
1313   }
1314   CreateEtaPtHisto2D("EtSimulatedChargedHadron","Simulated E_{T} from charged hadrons");
1315   CreateEtaPtHisto2D("EtNSimulatedPiPlus","Number of Simulated #pi^{+}");
1316   CreateEtaPtHisto2D("EtNSimulatedPiMinus","Number of simulated #pi^{-}");
1317   CreateEtaPtHisto2D("EtNSimulatedKPlus","Number of simulated K^{+}");
1318   CreateEtaPtHisto2D("EtNSimulatedKMinus","Number of simulated K^{-}");
1319   CreateEtaPtHisto2D("EtNSimulatedProton","Number of simulated p");
1320   CreateEtaPtHisto2D("EtNSimulatedAntiProton","Number of simulated #bar{p}");
1321   if(fBaryonEnhancement){
1322     CreateEtaPtHisto2D("EtNSimulatedProtonEnhanced","Number of simulated p");
1323     CreateEtaPtHisto2D("EtNSimulatedAntiProtonEnhanced","Number of simulated #bar{p}");
1324   }
1325   CreateEtaPtHisto2D("EtNSimulatedChargedHadron","Number of simulated charged hadrons");
1326   if(fDataSet==20100){//If this is Pb+Pb
1327     Int_t width = 5;
1328     if(fNCentBins<21) width = 10;
1329     for(Int_t i=0;i<fNCentBins;i++){
1330       CreateEtaPtHisto2D(Form("EtNSimulatedPiPlusCB%i",i),Form("Number of Simulated #pi^{+} for %i-%i central",i*width,(i+1)*width));
1331       CreateEtaPtHisto2D(Form("EtNSimulatedPiMinusCB%i",i),Form("Number of simulated #pi^{-} for %i-%i central",i*width,(i+1)*width));
1332       CreateEtaPtHisto2D(Form("EtNSimulatedKPlusCB%i",i),Form("Number of simulated K^{+} for %i-%i central",i*width,(i+1)*width));
1333       CreateEtaPtHisto2D(Form("EtNSimulatedKMinusCB%i",i),Form("Number of simulated K^{-} for %i-%i central",i*width,(i+1)*width));
1334       CreateEtaPtHisto2D(Form("EtNSimulatedProtonCB%i",i),Form("Number of simulated p for %i-%i central",i*width,(i+1)*width));
1335       CreateEtaPtHisto2D(Form("EtNSimulatedAntiProtonCB%i",i),Form("Number of simulated #bar{p} for %i-%i central",i*width,(i+1)*width));
1336       CreateEtaPtHisto2D(Form("EtNSimulatedChargedHadronCB%i",i),Form("Number of simulated charged hadrons for %i-%i central",i*width,(i+1)*width));
1337     }
1338   }
1339   CreateEtaPtHisto2D("EtSimulatedChargedHadronAssumingNoPt","Simulated E_{T} from charged hadrons assuming p_{T}=0");
1340   CreateEtaPtHisto2D("EtSimulatedChargedHadronAssumingPtTPCCut","Simulated E_{T} from charged hadrons assuming p_{T}=0.15");
1341   CreateEtaPtHisto2D("EtSimulatedChargedHadronAssumingPtITSCut","Simulated E_{T} from charged hadrons assuming p_{T}=0.10");
1342
1343   CreateEtaPtHisto2D("EtSimulatedChargedHadronAssumingPion","Simulated E_{T} from charged hadrons assuming they are all pions");
1344   CreateEtaPtHisto2D("EtSimulatedChargedHadronAssumingProton","Simulated E_{T} from charged hadrons assuming they are all pions");
1345   CreateEtaPtHisto2D("EtSimulatedChargedHadronAssumingKaon","Simulated E_{T} from charged hadrons assuming they are all pions");
1346   CreateEtaPtHisto2D("EtSimulatedKPlusAssumingPion","Simulated E_{T} from K^{+} assuming #pi mass");
1347   CreateEtaPtHisto2D("EtSimulatedKMinusAssumingPion","Simulated E_{T} from K^{-} assuming #pi mass");
1348   CreateEtaPtHisto2D("EtSimulatedProtonAssumingPion","Simulated E_{T} from p assuming #pi mass");
1349   CreateEtaPtHisto2D("EtSimulatedAntiProtonAssumingPion","Simulated E_{T} from #bar{p} assuming #pi mass");
1350   CreateEtaPtHisto2D("EtSimulatedKPlusAssumingProton","Simulated E_{T} from K^{+} assuming #pi mass");
1351   CreateEtaPtHisto2D("EtSimulatedKMinusAssumingProton","Simulated E_{T} from K^{-} assuming #pi mass");
1352   CreateEtaPtHisto2D("EtSimulatedPiPlusAssumingProton","Simulated E_{T} from p assuming #pi mass");
1353   CreateEtaPtHisto2D("EtSimulatedPiMinusAssumingProton","Simulated E_{T} from #bar{p} assuming #pi mass");
1354   CreateEtaPtHisto2D("EtSimulatedPiPlusAssumingKaon","Simulated E_{T} from K^{+} assuming #pi mass");
1355   CreateEtaPtHisto2D("EtSimulatedPiMinusAssumingKaon","Simulated E_{T} from K^{-} assuming #pi mass");
1356   CreateEtaPtHisto2D("EtSimulatedProtonAssumingKaon","Simulated E_{T} from p assuming #pi mass");
1357   CreateEtaPtHisto2D("EtSimulatedAntiProtonAssumingKaon","Simulated E_{T} from #bar{p} assuming #pi mass");
1358   if(fBaryonEnhancement){
1359     CreateEtaPtHisto2D("EtSimulatedProtonAssumingPionEnhanced","Simulated E_{T} from p assuming #pi mass");
1360     CreateEtaPtHisto2D("EtSimulatedAntiProtonAssumingPionEnhanced","Simulated E_{T} from #bar{p} assuming #pi mass");
1361   }
1362
1363   CreateEtaPtHisto2D("EtSimulatedLambda","Simulated E_{T} from #Lambda");
1364   CreateEtaPtHisto2D("EtSimulatedAntiLambda","Simulated E_{T} from #bar{#Lambda}");
1365   CreateEtaPtHisto2D("EtSimulatedK0S","Simulated E_{T} from K^{0}_{S}");
1366   CreateEtaPtHisto2D("EtSimulatedK0L","Simulated E_{T} from K^{0}_{L}");
1367   CreateEtaPtHisto2D("EtSimulatedLambdaReweighted","Simulated E_{T} from #Lambda");//These will also be used for baryon enhancement
1368   CreateEtaPtHisto2D("EtSimulatedAntiLambdaReweighted","Simulated E_{T} from #bar{#Lambda}");
1369   CreateEtaPtHisto2D("EtSimulatedK0SReweighted","Simulated E_{T} from K^{0}_{S}");
1370   CreateEtaPtHisto2D("EtSimulatedK0LReweighted","Simulated E_{T} from K^{0}_{L}");
1371   CreateEtaPtHisto2D("EtSimulatedNeutron","Simulated E_{T} from neutrons");
1372   CreateEtaPtHisto2D("EtSimulatedAntiNeutron","Simulated E_{T} from #bar{n}");
1373   CreateEtaPtHisto2D("EtSimulatedEPlus","Simulated E_{T} from e^{+}");
1374   CreateEtaPtHisto2D("EtSimulatedEMinus","Simulated E_{T} from e^{-}");
1375   CreateEtaPtHisto2D("EtSimulatedOmega","Simulated E_{T} from #Omega^{-}");
1376   CreateEtaPtHisto2D("EtSimulatedAntiOmega","Simulated E_{T} from #Omega^{+}");
1377   CreateEtaPtHisto2D("EtSimulatedXi","Simulated E_{T} from #Xi^{-}");
1378   CreateEtaPtHisto2D("EtSimulatedAntiXi","Simulated E_{T} from #Xi^{+}");
1379   CreateEtaPtHisto2D("EtSimulatedSigma","Simulated E_{T} from #Xi^{-}");
1380   CreateEtaPtHisto2D("EtSimulatedAntiSigma","Simulated E_{T} from #Xi^{+}");
1381   CreateEtaPtHisto2D("EtSimulatedXi0","Simulated E_{T} from #Xi^{0}");
1382   CreateEtaPtHisto2D("EtSimulatedAntiXi0","Simulated E_{T} from #Xi^{0}");
1383   CreateEtaPtHisto2D("EtSimulatedAllHadron","Simulated E_{T} from all hadrons");
1384
1385
1386   CreateEtaPtHisto2D("EtSimulatedLambdaDaughters","Simulated E_{T} from #Lambda Daughters");
1387   CreateEtaPtHisto2D("EtSimulatedAntiLambdaDaughters","Simulated E_{T} from #bar{#Lambda} Daughters");
1388   CreateEtaPtHisto2D("EtSimulatedK0SDaughters","Simulated E_{T} from K^{0}_{S} Daughters");
1389   CreateEtaPtHisto2D("EtSimulatedLambdaDaughtersReweighted","Simulated E_{T} from #Lambda Daughters");
1390   CreateEtaPtHisto2D("EtSimulatedAntiLambdaDaughtersReweighted","Simulated E_{T} from #bar{#Lambda} Daughters");
1391   CreateEtaPtHisto2D("EtSimulatedK0SDaughtersReweighted","Simulated E_{T} from K^{0}_{S} Daughters");
1392   CreateEtaPtHisto2D("EtSimulatedOmegaDaughters","Simulated E_{T} from #Omega^{-} Daughters");
1393   CreateEtaPtHisto2D("EtSimulatedAntiOmegaDaughters","Simulated E_{T} from #Omega^{+} Daughters");
1394   CreateEtaPtHisto2D("EtSimulatedXiDaughters","Simulated E_{T} from #Xi^{-} Daughters");
1395   CreateEtaPtHisto2D("EtSimulatedAntiXiDaughters","Simulated E_{T} from #Xi^{+} Daughters");
1396
1397
1398   CreateEtaPtHisto2D("EtSimulatedGamma","Simulated E_{T} from #gamma");
1399   CreateEtaPtHisto2D("EtSimulatedEta","Simulated E_{T} from #eta");
1400   CreateEtaPtHisto2D("EtSimulatedPi0","Simulated E_{T} from #pi^{0}");
1401   CreateEtaPtHisto2D("EtSimulatedOmega0","Simulated E_{T} from #omega");
1402
1403   TString *strTPC = new TString("TPC");
1404   TString *strITS = new TString("ITS");
1405   TString *strTPCITS = new TString("TPCITS");
1406   Int_t lastcutset = 1;
1407   if(fRequireITSHits) lastcutset = 2;
1408   for(Int_t i=0;i<=lastcutset;i++){
1409     TString *cutName = NULL;
1410     Float_t maxPtdEdx = 10;
1411     Float_t mindEdx = 35;
1412     Float_t maxdEdx = 150.0;
1413     switch(i){
1414     case 0:
1415       cutName = strTPC;
1416       break;
1417     case 1:
1418       cutName = strITS;
1419       maxPtdEdx = 5;
1420       maxdEdx = 500.0;
1421       break;
1422     case 2:
1423       cutName = strTPCITS;
1424       break;
1425     default:
1426       cerr<<"Error:  cannot make histograms!"<<endl;
1427       return;
1428     }
1429
1430     CreateEtaPtHisto2D(Form("EtReconstructed%sIdentifiedPiPlus",cutName->Data()),"Reconstructed E_{T} from identified #pi^{+}");
1431     CreateEtaPtHisto2D(Form("EtReconstructed%sIdentifiedPiMinus",cutName->Data()),"Reconstructed E_{T} from identified #pi^{-}");
1432     CreateEtaPtHisto2D(Form("EtReconstructed%sIdentifiedKPlus",cutName->Data()),"Reconstructed E_{T} from identified K^{+}");
1433     CreateEtaPtHisto2D(Form("EtReconstructed%sIdentifiedEMinus",cutName->Data()),"Reconstructed E_{T} from identified e^{-}");
1434     CreateEtaPtHisto2D(Form("EtReconstructed%sIdentifiedEPlus",cutName->Data()),"Reconstructed E_{T} from identified e^{+}");
1435     CreateEtaPtHisto2D(Form("EtReconstructed%sIdentifiedKMinus",cutName->Data()),"Reconstructed E_{T} from identified K^{-}");
1436     CreateEtaPtHisto2D(Form("EtReconstructed%sIdentifiedProton",cutName->Data()),"Reconstructed E_{T} from identified p");
1437     CreateEtaPtHisto2D(Form("EtReconstructed%sIdentifiedAntiProton",cutName->Data()),"Reconstructed E_{T} from identified #bar{p}");
1438     if(fBaryonEnhancement){
1439       CreateEtaPtHisto2D(Form("EtReconstructed%sIdentifiedProtonEnhanced",cutName->Data()),"Reconstructed E_{T} from identified p");
1440       CreateEtaPtHisto2D(Form("EtReconstructed%sIdentifiedAntiProtonEnhanced",cutName->Data()),"Reconstructed E_{T} from identified #bar{p}");
1441     }
1442     CreateEtaPtHisto2D(Form("EtNReconstructed%sUnidentified",cutName->Data()),"Number of Reconstructed unidentified particles");
1443     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedAssumingPion",cutName->Data()),"Reconstructed E_{T} from unidentified particles assuming pion mass");
1444     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedAssumingProton",cutName->Data()),"Reconstructed E_{T} from unidentified particles assuming pion mass");
1445     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedAssumingKaon",cutName->Data()),"Reconstructed E_{T} from unidentified particles assuming pion mass");
1446
1447     CreateEtaPtHisto2D(Form("EtNReconstructed%sUnidentifiedKaon",cutName->Data()),"Number of Reconstructed unidentified kaons particles");
1448     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedKaonAssumingPion",cutName->Data()),"Reconstructed E_{T} from unidentified kaons particles assuming pion mass");
1449     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedKaonAssumingProton",cutName->Data()),"Reconstructed E_{T} from unidentified kaons particles assuming pion mass");
1450     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedKaonAssumingKaon",cutName->Data()),"Reconstructed E_{T} from unidentified kaons particles assuming pion mass");
1451     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedKaon",cutName->Data()),"Reconstructed E_{T} from unidentified kaons particles assuming kaon mass");
1452     CreateEtaPtHisto2D(Form("EtNReconstructed%sUnidentifiedProton",cutName->Data()),"Number of Reconstructed unidentified proton particles");
1453     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedProtonAssumingPion",cutName->Data()),"Reconstructed E_{T} from unidentified proton particles assuming pion mass");
1454     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedProtonAssumingKaon",cutName->Data()),"Reconstructed E_{T} from unidentified proton particles assuming pion mass");
1455     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedPionAssumingKaon",cutName->Data()),"Reconstructed E_{T} from unidentified kaons particles assuming pion mass");
1456     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedProtonAssumingKaon",cutName->Data()),"Reconstructed E_{T} from unidentified proton particles assuming pion mass");
1457     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedPionAssumingProton",cutName->Data()),"Reconstructed E_{T} from unidentified kaons particles assuming pion mass");
1458     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedProtonAssumingProton",cutName->Data()),"Reconstructed E_{T} from unidentified proton particles assuming pion mass");
1459     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedProton",cutName->Data()),"Reconstructed E_{T} from unidentified proton particles assuming proton mass");
1460     if(fBaryonEnhancement){
1461       CreateEtaPtHisto2D(Form("EtNReconstructed%sUnidentifiedProtonEnhanced",cutName->Data()),"Number of Reconstructed unidentified proton particles");
1462       CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedProtonAssumingPionEnhanced",cutName->Data()),"Reconstructed E_{T} from unidentified proton particles assuming pion mass");
1463       CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedProtonEnhanced",cutName->Data()),"Reconstructed E_{T} from unidentified proton particles assuming proton mass");
1464     }
1465     CreateEtaPtHisto2D(Form("EtNReconstructed%sUnidentifiedPion",cutName->Data()),"Number of Reconstructed unidentified pions particles");
1466     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedPionAssumingPion",cutName->Data()),"Reconstructed E_{T} from unidentified pions particles assuming pion mass");
1467     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedPionAssumingKaon",cutName->Data()),"Reconstructed E_{T} from unidentified pions particles assuming pion mass");
1468     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedPionAssumingProton",cutName->Data()),"Reconstructed E_{T} from unidentified pions particles assuming pion mass");
1469     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentifiedPion",cutName->Data()),"Reconstructed E_{T} from unidentified pions particles assuming pion mass");
1470
1471     CreateEtaPtHisto2D(Form("EtReconstructed%sUnidentified",cutName->Data()),"Reconstructed E_{T} from unidentified particles using real mass");
1472     CreateEtaPtHisto2D(Form("EtReconstructed%sMisidentifiedElectrons",cutName->Data()),"Reconstructed E_{T} from misidentified electrons");
1473
1474
1475     CreateEtaPtHisto2D(Form("EtReconstructed%sPiPlus",cutName->Data()),"Reconstructed E_{T} from #pi^{+}");
1476     CreateEtaPtHisto2D(Form("EtReconstructed%sPiMinus",cutName->Data()),"Reconstructed E_{T} from #pi^{-}");
1477     CreateEtaPtHisto2D(Form("EtReconstructed%sKPlus",cutName->Data()),"Reconstructed E_{T} from K^{+}");
1478     CreateEtaPtHisto2D(Form("EtReconstructed%sKMinus",cutName->Data()),"Reconstructed E_{T} from K^{-}");
1479     CreateEtaPtHisto2D(Form("EtReconstructed%sProton",cutName->Data()),"Reconstructed E_{T} from p");
1480     CreateEtaPtHisto2D(Form("EtReconstructed%sAntiProton",cutName->Data()),"Reconstructed E_{T} from #bar{p}");
1481     if(fBaryonEnhancement){
1482       CreateEtaPtHisto2D(Form("EtReconstructed%sProtonEnhanced",cutName->Data()),"Reconstructed E_{T} from p");
1483       CreateEtaPtHisto2D(Form("EtReconstructed%sAntiProtonEnhanced",cutName->Data()),"Reconstructed E_{T} from #bar{p}");
1484     }
1485     CreateEtaPtHisto2D(Form("EtReconstructed%sChargedHadron",cutName->Data()),"Reconstructed E_{T} from charged hadrons");
1486     CreateEtaPtHisto2D(Form("EtNReconstructed%sPiPlus",cutName->Data()),"Reconstructed E_{T} from #pi^{+}");
1487     CreateEtaPtHisto2D(Form("EtNReconstructed%sPiMinus",cutName->Data()),"Reconstructed E_{T} from #pi^{-}");
1488     CreateEtaPtHisto2D(Form("EtNReconstructed%sKPlus",cutName->Data()),"Reconstructed E_{T} from K^{+}");
1489     CreateEtaPtHisto2D(Form("EtNReconstructed%sKMinus",cutName->Data()),"Reconstructed E_{T} from K^{-}");
1490     CreateEtaPtHisto2D(Form("EtNReconstructed%sProton",cutName->Data()),"Reconstructed E_{T} from p");
1491     CreateEtaPtHisto2D(Form("EtNReconstructed%sAntiProton",cutName->Data()),"Reconstructed E_{T} from #bar{p}");
1492     if(fBaryonEnhancement){
1493       CreateEtaPtHisto2D(Form("EtNReconstructed%sProtonEnhanced",cutName->Data()),"Reconstructed E_{T} from p");
1494       CreateEtaPtHisto2D(Form("EtNReconstructed%sAntiProtonEnhanced",cutName->Data()),"Reconstructed E_{T} from #bar{p}");
1495     }
1496     CreateEtaPtHisto2D(Form("EtNReconstructed%sChargedHadron",cutName->Data()),"Reconstructed E_{T} from charged hadrons");
1497     if(fDataSet==20100){//If this is Pb+Pb
1498       Int_t width = 5;
1499       if(fNCentBins<21) width = 10;
1500       for(Int_t j=0;j<fNCentBins;j++){
1501         CreateEtaPtHisto2D(Form("EtNReconstructed%sPiPlusCB%i",cutName->Data(),j),Form("Reconstructed E_{T} from #pi^{+} for %i-%i central",j*width,(j+1)*width));
1502         CreateEtaPtHisto2D(Form("EtNReconstructed%sPiMinusCB%i",cutName->Data(),j),Form("Reconstructed E_{T} from #pi^{-} for %i-%i central",j*width,(j+1)*width));
1503         CreateEtaPtHisto2D(Form("EtNReconstructed%sKPlusCB%i",cutName->Data(),j),Form("Reconstructed E_{T} from K^{+} for %i-%i central",j*width,(j+1)*width));
1504         CreateEtaPtHisto2D(Form("EtNReconstructed%sKMinusCB%i",cutName->Data(),j),Form("Reconstructed E_{T} from K^{-} for %i-%i central",j*width,(j+1)*width));
1505         CreateEtaPtHisto2D(Form("EtNReconstructed%sProtonCB%i",cutName->Data(),j),Form("Reconstructed E_{T} from p for %i-%i central",j*width,(j+1)*width));
1506         CreateEtaPtHisto2D(Form("EtNReconstructed%sAntiProtonCB%i",cutName->Data(),j),Form("Reconstructed E_{T} from #bar{p} for %i-%i central",j*width,(j+1)*width));
1507         CreateEtaPtHisto2D(Form("EtNReconstructed%sChargedHadronCB%i",cutName->Data(),j),Form("Reconstructed E_{T} from charged hadrons for %i-%i central",j*width,(j+1)*width));
1508       }
1509     }
1510
1511     CreateEtaPtHisto2D(Form("EtReconstructed%sChargedHadronAssumingPion",cutName->Data()),"Reconstructed E_{T} from charged hadrons assuming they are all pions");
1512     CreateEtaPtHisto2D(Form("EtReconstructed%sChargedHadronAssumingProton",cutName->Data()),"Reconstructed E_{T} from charged hadrons assuming they are all pions");
1513     CreateEtaPtHisto2D(Form("EtReconstructed%sChargedHadronAssumingKaon",cutName->Data()),"Reconstructed E_{T} from charged hadrons assuming they are all pions");
1514     CreateEtaPtHisto2D(Form("EtReconstructed%sKPlusAssumingPion",cutName->Data()),"Reconstructed E_{T} from K^{+} assuming #pi mass");
1515     CreateEtaPtHisto2D(Form("EtReconstructed%sKMinusAssumingPion",cutName->Data()),"Reconstructed E_{T} from K^{-} assuming #pi mass");
1516     CreateEtaPtHisto2D(Form("EtReconstructed%sProtonAssumingPion",cutName->Data()),"Reconstructed E_{T} from p assuming #pi mass");
1517     CreateEtaPtHisto2D(Form("EtReconstructed%sAntiProtonAssumingPion",cutName->Data()),"Reconstructed E_{T} from #bar{p} assuming #pi mass");
1518     CreateEtaPtHisto2D(Form("EtReconstructed%sPiPlusAssumingKaon",cutName->Data()),"Reconstructed E_{T} from K^{+} assuming #pi mass");
1519     CreateEtaPtHisto2D(Form("EtReconstructed%sPiMinusAssumingKaon",cutName->Data()),"Reconstructed E_{T} from K^{-} assuming #pi mass");
1520     CreateEtaPtHisto2D(Form("EtReconstructed%sKPlusAssumingKaon",cutName->Data()),"Reconstructed E_{T} from K^{+} assuming #pi mass");
1521     CreateEtaPtHisto2D(Form("EtReconstructed%sKMinusAssumingKaon",cutName->Data()),"Reconstructed E_{T} from K^{-} assuming #pi mass");
1522     CreateEtaPtHisto2D(Form("EtReconstructed%sProtonAssumingKaon",cutName->Data()),"Reconstructed E_{T} from p assuming #pi mass");
1523     CreateEtaPtHisto2D(Form("EtReconstructed%sAntiProtonAssumingKaon",cutName->Data()),"Reconstructed E_{T} from #bar{p} assuming #pi mass");
1524
1525     CreateEtaPtHisto2D(Form("EtReconstructed%sKPlusAssumingProton",cutName->Data()),"Reconstructed E_{T} from K^{+} assuming #pi mass");
1526     CreateEtaPtHisto2D(Form("EtReconstructed%sKMinusAssumingProton",cutName->Data()),"Reconstructed E_{T} from K^{-} assuming #pi mass");
1527     CreateEtaPtHisto2D(Form("EtReconstructed%sPiMinusAssumingProton",cutName->Data()),"Reconstructed E_{T} from p assuming #pi mass");
1528     CreateEtaPtHisto2D(Form("EtReconstructed%sPiPlusAssumingProton",cutName->Data()),"Reconstructed E_{T} from #bar{p} assuming #pi mass");
1529     CreateEtaPtHisto2D(Form("EtReconstructed%sProtonAssumingProton",cutName->Data()),"Reconstructed E_{T} from p assuming #pi mass");
1530     CreateEtaPtHisto2D(Form("EtReconstructed%sAntiProtonAssumingProton",cutName->Data()),"Reconstructed E_{T} from #bar{p} assuming #pi mass");
1531
1532     if(fBaryonEnhancement){
1533       CreateEtaPtHisto2D(Form("EtReconstructed%sProtonAssumingPionEnhanced",cutName->Data()),"Reconstructed E_{T} from p assuming #pi mass");
1534       CreateEtaPtHisto2D(Form("EtReconstructed%sAntiProtonAssumingPionEnhanced",cutName->Data()),"Reconstructed E_{T} from #bar{p} assuming #pi mass");
1535     }
1536
1537     CreateEtaPtHisto2D(Form("EtReconstructed%sEPlus",cutName->Data()),"Reconstructed E_{T} from e^{+}");
1538     CreateEtaPtHisto2D(Form("EtReconstructed%sEMinus",cutName->Data()),"Reconstructed E_{T} from e^{-}");
1539
1540
1541
1542     CreateEtaPtHisto2D(Form("EtReconstructed%sLambdaDaughters",cutName->Data()),"Reconstructed E_{T} from #Lambda Daughters");
1543   CreateEtaPtHisto2D(Form("EtReconstructed%sAntiLambdaDaughters",cutName->Data()),"Reconstructed E_{T} from #bar{#Lambda} Daughters");
1544   CreateEtaPtHisto2D(Form("EtReconstructed%sK0SDaughters",cutName->Data()),"Reconstructed E_{T} from K^{0}_{S} Daughters");
1545     CreateEtaPtHisto2D(Form("EtReconstructed%sLambdaDaughtersReweighted",cutName->Data()),"Reconstructed E_{T} from #Lambda Daughters");
1546   CreateEtaPtHisto2D(Form("EtReconstructed%sAntiLambdaDaughtersReweighted",cutName->Data()),"Reconstructed E_{T} from #bar{#Lambda} Daughters");
1547   CreateEtaPtHisto2D(Form("EtReconstructed%sK0SDaughtersReweighted",cutName->Data()),"Reconstructed E_{T} from K^{0}_{S} Daughters");
1548   CreateEtaPtHisto2D(Form("EtReconstructed%sOmegaDaughters",cutName->Data()),"Reconstructed E_{T} from #Omega^{-} Daughters");
1549   CreateEtaPtHisto2D(Form("EtReconstructed%sAntiOmegaDaughters",cutName->Data()),"Reconstructed E_{T} from #Omega^{+} Daughters");
1550   CreateEtaPtHisto2D(Form("EtReconstructed%sXiDaughters",cutName->Data()),"Reconstructed E_{T} from #Xi^{-} Daughters");
1551   CreateEtaPtHisto2D(Form("EtReconstructed%sAntiXiDaughters",cutName->Data()),"Reconstructed E_{T} from #Xi^{+} Daughters");
1552
1553     CreateIntHisto1D(Form("UnidentifiedPIDs%s",cutName->Data()),"PIDs of unidentified particles", "PID", "Number of particles",9, -4,4);
1554     CreateHisto2D(Form("MisidentifiedPIDs%s",cutName->Data()),"PIDs of misidentified particles", "PID real","PID identified",5, -.5,4.5,5, -.5,4.5);
1555     CreateHisto2D(Form("dEdxAll%s",cutName->Data()),"dE/dx for all particles","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
1556     CreateHisto2D(Form("dEdxPion%s",cutName->Data()),"dE/dx for #pi^{#pm}","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
1557     CreateHisto2D(Form("dEdxKaon%s",cutName->Data()),"dE/dx for K^{#pm}","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
1558     CreateHisto2D(Form("dEdxProton%s",cutName->Data()),"dE/dx for p(#bar{p})","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
1559     CreateHisto2D(Form("dEdxElectron%s",cutName->Data()),"dE/dx for e^{#pm}","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
1560     CreateHisto2D(Form("dEdxUnidentified%s",cutName->Data()),"dE/dx for unidentified particles","momentum (GeV/c)","dE/dx",400,0.0,maxPtdEdx,200,mindEdx,maxdEdx);
1561   }
1562   delete strTPC;
1563   delete strITS;
1564   delete strTPCITS;
1565
1566   Float_t minEt = 0.0;
1567   Float_t maxEt = 100.0;
1568   if(fDataSet==20100) maxEt=3000.0;
1569   Int_t nbinsEt = 100;
1570   char histoname[200];
1571   char histotitle[200];
1572   char xtitle[50];
1573   char ytitle[50];
1574   TString *sTPC = new TString("TPC");
1575   TString *sITS = new TString("ITS");
1576   TString *sTPCpt = new TString("0.15");
1577   TString *sITSpt = new TString("0.10");
1578   TString *sPID = new TString("");
1579   TString *sNoPID = new TString("NoPID");
1580   TString *sNoPIDString = new TString(", No PID");
1581   TString *sHadEt = new TString("HadEt");
1582   TString *sTotEt = new TString("TotEt");
1583   TString *sTotEtString = new TString("total E_{T}");
1584   TString *sHadEtString = new TString("hadronic E_{T}");
1585   TString *sFull = new TString("Full");
1586   TString *sEMCAL = new TString("EMCAL");
1587   TString *sPHOS = new TString("PHOS");
1588   float etDiff = 1.5;
1589
1590   for(int tpc = 0;tpc<lastcutset;tpc++){
1591     TString *detector = NULL;
1592     TString *ptstring = NULL;
1593     if(tpc==1) {detector = sTPC; ptstring = sTPCpt;}
1594     else{detector = sITS; ptstring = sITSpt;}
1595     for(int hadet = 0;hadet<2;hadet++){
1596       TString *et = NULL;
1597       TString *etstring = NULL;
1598       if(hadet==1) {et = sHadEt; etstring = sHadEtString;}
1599       else{et = sTotEt; etstring = sTotEtString;}
1600       for(int type = 0;type<3;type++){
1601         if(type==0 && !fInvestigateFull) continue;
1602         if(type==1 && !fInvestigateEMCal) continue;
1603         if(type==2 && !fInvestigatePHOS) continue;
1604         TString *acceptance = NULL;
1605         switch(type){
1606         case 0:
1607           acceptance = sFull;
1608           etDiff = 1.5;
1609           break;
1610         case 1:
1611           acceptance = sEMCAL;
1612           etDiff = 5;
1613           break;
1614         case 2:
1615           acceptance = sPHOS;
1616           etDiff = 5;
1617           break;
1618         default:
1619           acceptance = sFull;
1620         }
1621         for(int pid = 0;pid<2;pid++){
1622           TString *partid = NULL;
1623           TString *partidstring = NULL;
1624           if(pid==1){partid = sPID; partidstring = sPID;}
1625           else{partid = sNoPID; partidstring = sNoPIDString;}
1626           snprintf(histoname,200,"Sim%sMinusReco%s%sAcceptance%s%s",et->Data(),et->Data(),acceptance->Data(),detector->Data(),partid->Data());
1627           snprintf(histotitle,200,"(Simulated %s - reconstructed %s)/(Simulated %s) with %s acceptance for p_{T}>%s GeV/c%s",etstring->Data(),etstring->Data(),etstring->Data(),acceptance->Data(),ptstring->Data(),partidstring->Data());
1628           snprintf(ytitle,50,"(Simulated %s - reconstructed %s)/(Simulated %s)",etstring->Data(),etstring->Data(),etstring->Data());
1629           snprintf(xtitle,50,"Simulated %s",etstring->Data());
1630           CreateHisto2D(histoname,histotitle,xtitle,ytitle,nbinsEt,minEt,maxEt,nbinsEt,-etDiff,etDiff);
1631           if(hadet==0 && type==0 && fInvestigatePiKP){//we only want to do this once...  not the most elegant way of coding but hey...
1632             snprintf(histoname,200,"SimPiKPMinusRecoPiKP%sAcceptance%s%s",acceptance->Data(),detector->Data(),partid->Data());
1633             snprintf(histotitle,200,"(Sim PiKP - reco PiKP)/(Sim PiKP) with %s acceptance for p_{T}>%s GeV/c%s",acceptance->Data(),ptstring->Data(),partidstring->Data());
1634             snprintf(ytitle,50,"(Sim PiKP - reco PiKP)/(Sim PiKP)");
1635             snprintf(xtitle,50,"Simulated E_{T}^{#pi,K,p}");
1636             CreateHisto2D(histoname,histotitle,xtitle,ytitle,nbinsEt,minEt,maxEt,nbinsEt,-etDiff,etDiff);
1637           }
1638         }
1639       }
1640     }
1641   }
1642   CreateHisto1D("SimPiKPEt","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt,minEt,maxEt);
1643   CreateHisto1D("SimTotEt","Simulated Total E_{T}","Simulated Total E_{T}","Number of events",nbinsEt*4,minEt,maxEt);
1644   CreateHisto1D("SimHadEt","Simulated Hadronic E_{T}","Simulated Hadronic E_{T}","Number of events",nbinsEt*4,minEt,maxEt);
1645   if(fDataSet==20100){
1646     Int_t width = 5;
1647     if(fNCentBins<21) width = 10;
1648     for(Int_t j=0;j<fNCentBins;j++){
1649       CreateHisto1D(Form("SimTotEtCB%i",j),Form("Simulated Total E_{T} for %i-%i central",j*width,(j+1)*width),"Simulated Total E_{T}","Number of events",nbinsEt*4,minEt,maxEt);
1650       CreateHisto1D(Form("SimHadEtCB%i",j),Form("Simulated Hadronic E_{T} for %i-%i central",j*width,(j+1)*width),"Simulated Hadronic E_{T}","Number of events",nbinsEt*4,minEt,maxEt);
1651       CreateHisto1D(Form("SimPiKPEtCB%i",j),Form("Simulated #pi,K,p E_{T} for %i-%i central",j*width,(j+1)*width),"Simulated #pi,K,p E_{T}","Number of events",nbinsEt,minEt,maxEt);
1652     }
1653   }
1654
1655   etDiff = 0.15;
1656
1657   if(fInvestigateSmearing){
1658
1659     //======================================================================
1660
1661     snprintf(histoname,200,"SimPiKPEtMinusSimPtSmeared");
1662     snprintf(histotitle,200,"Simulated (true-smeared)/true for 0.5 percent momentum smearing");
1663     snprintf(ytitle,50,"(true-smeared)/true");
1664     snprintf(xtitle,50,"true p, K, p E_{T}");
1665     CreateHisto2D(histoname,histotitle,xtitle,ytitle,nbinsEt,minEt,maxEt,nbinsEt,-etDiff,etDiff);
1666     snprintf(histoname,200,"SimPiKPEtPtSmeared");
1667     snprintf(histotitle,200,"Simulated E_{T} for 0.5 percent momentum smearing");
1668     snprintf(ytitle,50,"Number of events");
1669     snprintf(xtitle,50,"p, K, p E_{T}");
1670     CreateHisto1D(histoname,histotitle,xtitle,ytitle,nbinsEt,minEt,maxEt);
1671
1672     //======================================================================
1673
1674     snprintf(histoname,200,"SimPiKPEtMinusSimEfficiencySmeared");
1675     snprintf(histotitle,200,"Simulated (true-smeared)/true for efficiency smearing");
1676     snprintf(ytitle,50,"(true-smeared)/true");
1677     snprintf(xtitle,50,"true p, K, p E_{T}");
1678     CreateHisto2D(histoname,histotitle,xtitle,ytitle,nbinsEt,minEt,maxEt,nbinsEt,-etDiff*5,etDiff*5);
1679     snprintf(histoname,200,"SimPiKPEtEfficiencySmeared");
1680     snprintf(histotitle,200,"Simulated E_{T} for efficiency smearing");
1681     snprintf(ytitle,50,"Number of events");
1682     snprintf(xtitle,50,"p, K, p E_{T}");
1683     CreateHisto1D(histoname,histotitle,xtitle,ytitle,nbinsEt,minEt,maxEt);
1684
1685     //======================================================================
1686
1687     snprintf(histoname,200,"SimPiKPEtMinusSimPtCutSmearedTPC");
1688     snprintf(histotitle,200,"Simulated (true-smeared)/true for p_{T}>0.15 GeV/c smearing");
1689     snprintf(ytitle,50,"(true-smeared)/true");
1690     snprintf(xtitle,50,"true p, K, p E_{T}");
1691     CreateHisto2D(histoname,histotitle,xtitle,ytitle,nbinsEt,minEt,maxEt,nbinsEt,-etDiff,etDiff);
1692     snprintf(histoname,200,"SimPiKPEtPtCutSmearedTPC");
1693     snprintf(histotitle,200,"Simulated E_{T} for p_{T}>0.15 GeV/c smearing");
1694     snprintf(ytitle,50,"Number of events");
1695     snprintf(xtitle,50,"p, K, p E_{T}");
1696     CreateHisto1D(histoname,histotitle,xtitle,ytitle,nbinsEt,minEt,maxEt);
1697
1698
1699     //======================================================================
1700
1701     snprintf(histoname,200,"SimPiKPEtMinusSimPtCutSmearedITS");
1702     snprintf(histotitle,200,"Simulated (true-smeared)/true for p_{T}>0.10 GeV/c smearing");
1703     snprintf(ytitle,50,"(true-smeared)/true");
1704     snprintf(xtitle,50,"true p, K, p E_{T}");
1705     CreateHisto2D(histoname,histotitle,xtitle,ytitle,nbinsEt,minEt,maxEt,nbinsEt,-etDiff,etDiff);
1706     snprintf(histoname,200,"SimPiKPEtPtCutSmearedITS");
1707     snprintf(histotitle,200,"Simulated E_{T} for p_{T}>0.10 GeV/c smearing");
1708     snprintf(ytitle,50,"Number of events");
1709     snprintf(xtitle,50,"p, K, p E_{T}");
1710     CreateHisto1D(histoname,histotitle,xtitle,ytitle,nbinsEt,minEt,maxEt);
1711
1712     //======================================================================
1713
1714     snprintf(histoname,200,"SimPiKPEtMinusSimPIDSmeared");
1715     snprintf(histotitle,200,"Simulated (true-smeared)/true for PID smearing");
1716     snprintf(ytitle,50,"(true-smeared)/true");
1717     snprintf(xtitle,50,"true p, K, p E_{T}");
1718     CreateHisto2D(histoname,histotitle,xtitle,ytitle,nbinsEt,minEt,maxEt,nbinsEt,-etDiff,etDiff);
1719     snprintf(histoname,200,"SimPiKPEtPIDSmeared");
1720     snprintf(histotitle,200,"Simulated E_{T} for PID smearing");
1721     snprintf(ytitle,50,"Number of events");
1722     snprintf(xtitle,50,"p, K, p E_{T}");
1723     CreateHisto1D(histoname,histotitle,xtitle,ytitle,nbinsEt,minEt,maxEt);
1724
1725     //======================================================================
1726
1727     snprintf(histoname,200,"SimPiKPEtMinusSimPIDSmearedNoID");
1728     snprintf(histotitle,200,"Simulated (true-smeared)/true for PID smearing No ID");
1729     snprintf(ytitle,50,"(true-smeared)/true");
1730     snprintf(xtitle,50,"true p, K, p E_{T}");
1731     CreateHisto2D(histoname,histotitle,xtitle,ytitle,nbinsEt,minEt,maxEt,nbinsEt,-etDiff,etDiff);
1732     snprintf(histoname,200,"SimPiKPEtPIDSmearedNoID");
1733     snprintf(histotitle,200,"Simulated E_{T} for PID smearing No ID");
1734     snprintf(ytitle,50,"Number of events");
1735     snprintf(xtitle,50,"p, K, p E_{T}");
1736     CreateHisto1D(histoname,histotitle,xtitle,ytitle,nbinsEt,minEt,maxEt);
1737   }
1738   delete sTPC;
1739   delete sITS;
1740   delete sTPCpt;
1741   delete sITSpt;
1742   delete sPID;
1743   delete sNoPID;
1744   delete sNoPIDString;
1745   delete sHadEt;
1746   delete sTotEt;
1747   delete sTotEtString;
1748   delete sHadEtString;
1749   delete sFull;
1750   delete sEMCAL;
1751   delete sPHOS;
1752   CreateIntHisto1D("NEvents","Number of events","number of events","Number of events",1,0,1);
1753   CreateResolutionPtHisto2D("presolutionTPC","p resolution","p^{rec}","(p^{sim}-p^{rec})/p^{rec}");
1754   CreateResolutionPtHisto2D("pTresolutionTPC","p_{T} resolution","p_{T}^{rec}","(p_{T}^{sim}-p_{T}^{rec})/p_{T}^{rec}");
1755   CreateResolutionPtHisto2D("ETresolutionTPC","E_{T} resolution","E_{T}^{rec}","(E_{T}^{sim}-E_{T}^{rec})/E_{T}^{rec}");
1756   CreateResolutionPtHisto2D("pTresolutionTPCITS","p_{T} resolution","p_{T}^{rec}","(p_{T}^{sim}-p_{T}^{rec})/p_{T}^{rec}");
1757   CreateResolutionPtHisto2D("ETresolutionTPCITS","E_{T} resolution","E_{T}^{rec}","(E_{T}^{sim}-E_{T}^{rec})/E_{T}^{rec}");
1758   CreateResolutionPtHisto2D("presolutionTPCITS","p resolution","p^{rec}","(p^{sim}-p^{rec})/p^{rec}");
1759   CreateResolutionPtHisto2D("pTresolutionITS","p_{T} resolution","p_{T}^{rec}","(p_{T}^{sim}-p_{T}^{rec})/p_{T}^{rec}");
1760   CreateResolutionPtHisto2D("ETresolutionITS","E_{T} resolution","E_{T}^{rec}","(E_{T}^{sim}-E_{T}^{rec})/E_{T}^{rec}");
1761   CreateResolutionPtHisto2D("presolutionITS","p resolution","p^{rec}","(p^{sim}-p^{rec})/p^{rec}");
1762   CreatePtHisto1D("pTsimITS","p_{T}^{sim}","p_{T}^{sim}","Number of particles");
1763   CreatePtHisto1D("pTsimTPC","p_{T}^{sim}","p_{T}^{sim}","Number of particles");
1764   CreatePtHisto1D("pTsimTPCITS","p_{T}^{sim}","p_{T}^{sim}","Number of particles");
1765   CreatePtHisto1D("pTrecITS","p_{T}^{rec}","p_{T}^{rec}","Number of particles");
1766   CreatePtHisto1D("pTrecTPC","p_{T}^{rec}","p_{T}^{rec}","Number of particles");
1767   CreatePtHisto1D("pTrecTPCITS","p_{T}^{rec}","p_{T}^{rec}","Number of particles");
1768   if(fDataSet==20100){
1769     Int_t width = 5;
1770     if(fNCentBins<21) width = 10;
1771     for(Int_t j=0;j<fNCentBins;j++){
1772       CreatePtHisto1D(Form("pTsimITSCB%i",j),Form("p_{T}^{sim} for %i-%i central",j*width,(j+1)*width),"p_{T}^{sim}","Number of particles");
1773       CreatePtHisto1D(Form("pTsimTPCITSCB%i",j),Form("p_{T}^{sim} for %i-%i central",j*width,(j+1)*width),"p_{T}^{sim}","Number of particles");
1774       CreatePtHisto1D(Form("pTsimTPCCB%i",j),Form("p_{T}^{sim} for %i-%i central",j*width,(j+1)*width),"p_{T}^{sim}","Number of particles");
1775       CreatePtHisto1D(Form("pTrecITSCB%i",j),Form("p_{T}^{rec} for %i-%i central",j*width,(j+1)*width),"p_{T}^{rec}","Number of particles");
1776       CreatePtHisto1D(Form("pTrecTPCITSCB%i",j),Form("p_{T}^{rec} for %i-%i central",j*width,(j+1)*width),"p_{T}^{rec}","Number of particles");
1777       CreatePtHisto1D(Form("pTrecTPCCB%i",j),Form("p_{T}^{rec} for %i-%i central",j*width,(j+1)*width),"p_{T}^{rec}","Number of particles");
1778     }
1779   }
1780
1781 }
1782