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