]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGLF/totEt/AliAnalysisHadEtMonteCarlo.cxx
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGLF / totEt / AliAnalysisHadEtMonteCarlo.cxx
index 7ec0bb7910d345c02f9d56f5d7be939218aaffc5..77f31b54a8f750ccd5d15bff3ea67799d89631f9 100644 (file)
 #include "AliCentrality.h"
 #include "AliLog.h"
 #include "AliPWG0Helper.h"
+#include "AliPIDResponse.h"
+#include "AliTPCPIDResponse.h" 
+#include "AliInputEventHandler.h"
+#include "AliAnalysisManager.h"
 //class AliPWG0Helper;
 //#include "$ALICE_ROOT/PWG0/AliPWG0Helper.h"
 
@@ -42,6 +46,8 @@ Float_t AliAnalysisHadEtMonteCarlo::fgSmearWidths[4] = {0.005,0.006,0.007,0.008}
 
 AliAnalysisHadEtMonteCarlo::AliAnalysisHadEtMonteCarlo():AliAnalysisHadEt()
                                                        ,fSimPiKPEt(0)
+                                                       ,fSimRawEtTPC(0)
+                                                       ,fSimRawEtITS(0)
                                                        ,fSimHadEt(0)
                                                        ,fSimTotEt(0) 
                                                        ,fSimPiKPEtShouldBeReco(0)
@@ -74,6 +80,8 @@ void AliAnalysisHadEtMonteCarlo::ResetEventValues(){//resetting event variables
     fSimHadEt=0.0;
     fSimTotEt=0.0;
     fSimPiKPEt=0.0;
+    fSimRawEtTPC=0.0;
+    fSimRawEtITS=0.0;
 }
 Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
 { // analyse MC and real event info
@@ -166,19 +174,27 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
          }
        else{
          Float_t nSigmaPion,nSigmaProton,nSigmaKaon,nSigmaElectron;
-         pID->MakeTPCPID(track);
-         pID->MakeITSPID(track);
+//       pID->MakeTPCPID(track);
+//       pID->MakeITSPID(track);
          if(cutset!=1){
-           nSigmaPion = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kPion));
-           nSigmaProton = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kProton));
-           nSigmaKaon = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kKaon));
-           nSigmaElectron = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kElectron));
+           nSigmaPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion)); 
+           nSigmaProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton)); 
+           nSigmaKaon =TMath::Abs( fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)); 
+           nSigmaElectron =TMath::Abs( fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron)); 
+//         nSigmaPion = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kPion));
+//         nSigmaProton = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kProton));
+//         nSigmaKaon = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kKaon));
+//         nSigmaElectron = TMath::Abs(pID->NumberOfSigmasTPC(track,AliPID::kElectron));
          }
          else{
-           nSigmaPion = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kPion));
-           nSigmaProton = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kProton));
-           nSigmaKaon = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kKaon));
-           nSigmaElectron = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kElectron));
+           nSigmaPion = TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kPion)); 
+           nSigmaProton = TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kProton)); 
+           nSigmaKaon = TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kKaon)); 
+           nSigmaElectron = TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kElectron)); 
+//         nSigmaPion = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kPion));
+//         nSigmaProton = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kProton));
+//         nSigmaKaon = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kKaon));
+//         nSigmaElectron = TMath::Abs(pID->NumberOfSigmasITS(track,AliPID::kElectron));
          }
 //       bool isPion = (nSigmaPion<3.0 && nSigmaProton>2.0 && nSigmaKaon>2.0);
 //       bool isElectron = (nSigmaElectron<2.0 && nSigmaPion>4.0 && nSigmaProton>3.0 && nSigmaKaon>3.0);
@@ -203,7 +219,7 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
          FillHisto2D(Form("dEdxAll%s",cutName->Data()),track->P(),dEdx,1.0);
 
          UInt_t label = (UInt_t)TMath::Abs(track->GetLabel());
-         TParticle  *simPart  = stack->Particle(label);
+           TParticle  *simPart  = stack->Particle(label);
          if(!simPart) {
            Printf("no MC particle\n");                 
            continue;           
@@ -232,23 +248,24 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
              //for calculating et as it's done in the reconstructed data
              Float_t corrBkgd=0.0;
              Float_t corrNotID=0.0;
-             Float_t corrNoID=0.0;// = fHadEtReco->GetCorrections()->GetNotIDCorrectionNoPID(track->Pt());
+             //Float_t corrNoID=0.0;// = fHadEtReco->GetCorrections()->GetNotIDCorrectionNoPID(track->Pt());
              Float_t corrEff = 0.0;
-             Float_t corrEffNoID = 0.0;
+             //Float_t corrEffNoID = 0.0;
              Float_t et = 0.0;
              if(cutset==2){//TPC
                corrBkgd = fHadEtReco->GetCorrections()->GetBackgroundCorrectionTPC(track->Pt());
-               corrEffNoID = fHadEtReco->GetCorrections()->GetTPCEfficiencyCorrectionHadron(track->Pt(),fCentBin);
+               //corrEffNoID = fHadEtReco->GetCorrections()->GetTPCEfficiencyCorrectionHadron(track->Pt(),fCentBin);
                corrNotID = fHadEtReco->GetCorrections()->GetNotIDConstCorrectionTPC();
-               corrNoID = fHadEtReco->GetCorrections()->GetNotIDConstCorrectionTPCNoID();
+               //corrNoID = fHadEtReco->GetCorrections()->GetNotIDConstCorrectionTPCNoID();
              }
              if(cutset==1){//ITS
                corrBkgd = fHadEtReco->GetCorrections()->GetBackgroundCorrectionITS(track->Pt());
-               corrEffNoID = fHadEtReco->GetCorrections()->GetITSEfficiencyCorrectionHadron(track->Pt(),fCentBin);
+               //corrEffNoID = fHadEtReco->GetCorrections()->GetITSEfficiencyCorrectionHadron(track->Pt(),fCentBin);
                corrNotID = fHadEtReco->GetCorrections()->GetNotIDConstCorrectionITS();
-               corrNoID = fHadEtReco->GetCorrections()->GetNotIDConstCorrectionITSNoID();
+               //corrNoID = fHadEtReco->GetCorrections()->GetNotIDConstCorrectionITSNoID();
              }
              
+             
              bool isprimary = stack->IsPhysicalPrimary(label);
              if (TMath::Abs(track->Eta()) < fHadEtReco->GetCorrections()->GetEtaCut()){
                  if(isPion){
@@ -294,16 +311,21 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                  if(!isPion && !isProton && !isKaon && !unidentified){
                      eTBkgdAsReconstructed += et*corrBkgd*corrEff*corrNotID;
                  }
-                 Int_t pdgCode =  simPart->GetPDG(0)->PdgCode();
-                 if(pdgCode==fgPiPlusCode ||pdgCode==fgPiMinusCode){eTtotalAsReconstructedPi+=et*corrBkgd*corrEff*corrNotID;}
-                 if(pdgCode==fgKPlusCode ||pdgCode==fgKMinusCode){eTtotalAsReconstructedK+=et*corrBkgd*corrEff*corrNotID;}
-                 if(pdgCode==fgProtonCode ||pdgCode==fgAntiProtonCode){eTtotalAsReconstructedP+=et*corrBkgd*corrEff*corrNotID;}
+                 TParticlePDG *pdg = simPart->GetPDG(0);
+                 if(pdg){
+                   Int_t pdgCode =  simPart->GetPDG(0)->PdgCode();
+                   if(pdgCode==fgPiPlusCode ||pdgCode==fgPiMinusCode){eTtotalAsReconstructedPi+=et*corrBkgd*corrEff*corrNotID;}
+                   if(pdgCode==fgKPlusCode ||pdgCode==fgKMinusCode){eTtotalAsReconstructedK+=et*corrBkgd*corrEff*corrNotID;}
+                   if(pdgCode==fgProtonCode ||pdgCode==fgAntiProtonCode){eTtotalAsReconstructedP+=et*corrBkgd*corrEff*corrNotID;}
+                 }
                }
            }
 
            if(cutset==2) eTtotalSimAll += Et(simPart);
            if(stack->IsPhysicalPrimary(label)){
              if (TMath::Abs(simPart->Eta()) < fHadEtReco->GetCorrections()->GetEtaCut()){
+               TParticlePDG *pdg = simPart->GetPDG(0);
+               if(!pdg) continue;
                Int_t pdgCode =  simPart->GetPDG(0)->PdgCode();
                Int_t mypid = 0;
                if(pdgCode==AliAnalysisHadEt::fgPiPlusCode) mypid = 1;
@@ -314,7 +336,7 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                if(pdgCode==fgAntiProtonCode) mypid = 2;
                if(pdgCode==fgKMinusCode) mypid = 3;
                if(pdgCode==fgEMinusCode) mypid = 4;
-               bool filled = false;      
+               //bool filled = false;      
                //for smearing investigations
                if(fInvestigateSmearing && cutset==2){
                  pTtotalReco += simPart->Pt();
@@ -473,7 +495,7 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                    FillHisto2D(Form("EtReconstructed%sPiPlusAssumingKaon",cutName->Data()),pT,eta,myEtK);
                    FillHisto2D(Form("EtReconstructed%sPiPlusAssumingProton",cutName->Data()),pT,eta,myEtP);
                  }
-                 filled = true;
+                 //filled = true;
                }
                if(pdgCode == fgPiMinusCode){
                  float myEt = Et(simPart);
@@ -506,7 +528,7 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                    FillHisto2D(Form("EtReconstructed%sPiMinusAssumingKaon",cutName->Data()),pT,eta,myEtK);
                    FillHisto2D(Form("EtReconstructed%sPiMinusAssumingProton",cutName->Data()),pT,eta,myEtP);
                  }
-                 filled = true;
+                 //filled = true;
                }
                if(pdgCode == fgKPlusCode){
                  float myEt = Et(simPart);
@@ -540,7 +562,7 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                    FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingProton",cutName->Data()),pT,eta,myEtP);
                    FillHisto2D(Form("EtReconstructed%sKPlusAssumingProton",cutName->Data()),pT,eta,myEtP);
                  }
-                 filled = true;
+                 //filled = true;
                }
                if(pdgCode == fgKMinusCode){
                  float myEt = Et(simPart);
@@ -574,7 +596,7 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                    FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingProton",cutName->Data()),pT,eta,myEtP);
                    FillHisto2D(Form("EtReconstructed%sKMinusAssumingProton",cutName->Data()),pT,eta,myEtP);
                  }
-                 filled = true;
+                 //filled = true;
                }
                if(pdgCode == fgProtonCode){
                  float myEt = Et(simPart);
@@ -608,7 +630,7 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                    FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingProton",cutName->Data()),pT,eta,myEt);
                    FillHisto2D(Form("EtReconstructed%sProtonAssumingProton",cutName->Data()),pT,eta,myEt);
                  }
-                 filled = true;
+                 //filled = true;
 
                  if( !fRunLightweight){
                    if(fBaryonEnhancement){
@@ -652,7 +674,7 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                    FillHisto2D(Form("EtReconstructed%sChargedHadronAssumingProton",cutName->Data()),pT,eta,myEt);
                    FillHisto2D(Form("EtReconstructed%sAntiProtonAssumingProton",cutName->Data()),pT,eta,myEt);
                  }
-                 filled = true;
+                 //filled = true;
                  if( !fRunLightweight){
                    if(fBaryonEnhancement){
                      float enhancement = ProtonBaryonEnhancement(track->Pt());
@@ -671,7 +693,7 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                      FillHisto2D(Form("EtReconstructed%sMisidentifiedElectrons",cutName->Data()),track->Pt(),track->Eta(),myEtPi);
                    }
                  }
-                 filled = true;
+                 //filled = true;
                }
                if(pdgCode == fgEMinusCode){
                  if( !fRunLightweight){
@@ -682,7 +704,7 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                    float myEt = Et(simPart);
                    FillHisto2D(Form("EtReconstructed%sEMinus",cutName->Data()),simPart->Pt(),simPart->Eta(),myEt);
                  }
-                 filled = true;
+                 //filled = true;
                }
                if( !fRunLightweight){
                  if(myEtReco>0.0){FillHisto2D(Form("ETresolution%s",cutName->Data()),myEtReco,(myEtSim-myEtReco)/myEtReco,1.0);}
@@ -705,7 +727,8 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
              if(isProton) myrecoEt = Et(track->P(),track->Theta(),fgProtonCode,track->Charge());
              if(isKaon) myrecoEt = Et(track->P(),track->Theta(),fgKPlusCode,track->Charge());
              if (TMath::Abs(simPart->Eta()) < fHadEtReco->GetCorrections()->GetEtaCut()){
-               TParticle *mom = stack->Particle(simPart->GetFirstMother());
+               TParticle *mom = NULL;
+               if(simPart->GetFirstMother()>=0) mom= stack->Particle(simPart->GetFirstMother());
                if(mom){
                  TParticlePDG *pc = mom->GetPDG(0);
                  if(pc){
@@ -715,6 +738,8 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                      float myEt = Et(simPart);
                      float pT = simPart->Pt();
                      float eta = simPart->Eta();
+                     TParticlePDG *simpdg = simPart->GetPDG(0);
+                     if(!simpdg) continue;
                      eTtotalRecoBkgd+=myEt;
                      if(fUseRecoPt){//Then we switch the pT and the Et
                        myEt = Et(track->P(),track->Theta(),simPart->GetPDG(0)->PdgCode(),track->Charge());
@@ -737,6 +762,8 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                      float pT = simPart->Pt();
                      float eta = simPart->Eta();
                      eTtotalRecoBkgd+=myEt;
+                     TParticlePDG *simpdg = simPart->GetPDG(0);
+                     if(!simpdg) continue;
                      if(fUseRecoPt){//Then we switch the pT and the Et
                        myEt = Et(track->P(),track->Theta(),simPart->GetPDG(0)->PdgCode(),track->Charge());
                        pT = track->Pt();
@@ -758,6 +785,8 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                      float pT = simPart->Pt();
                      float eta = simPart->Eta();
                      eTtotalRecoBkgd+=myEt;
+                     TParticlePDG *simpdg = simPart->GetPDG(0);
+                     if(!simpdg) continue;
                      if(fUseRecoPt){//Then we switch the pT and the Et
                        myEt = Et(track->P(),track->Theta(),simPart->GetPDG(0)->PdgCode(),track->Charge());
                        pT = track->Pt();
@@ -775,6 +804,8 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                      float pT = simPart->Pt();
                      float eta = simPart->Eta();
                      eTtotalRecoBkgd+=myEt;
+                     TParticlePDG *simpdg = simPart->GetPDG(0);
+                     if(!simpdg) continue;
                      if(fUseRecoPt){//Then we switch the pT and the Et
                        myEt = Et(track->P(),track->Theta(),simPart->GetPDG(0)->PdgCode(),track->Charge());
                        pT = track->Pt();
@@ -790,6 +821,8 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                      float pT = simPart->Pt();
                      float eta = simPart->Eta();
                      eTtotalRecoBkgd+=myEt;
+                     TParticlePDG *simpdg = simPart->GetPDG(0);
+                     if(!simpdg) continue;
                      if(fUseRecoPt){//Then we switch the pT and the Et
                        myEt = Et(track->P(),track->Theta(),simPart->GetPDG(0)->PdgCode(),track->Charge());
                        pT = track->Pt();
@@ -805,6 +838,8 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                      float pT = simPart->Pt();
                      float eta = simPart->Eta();
                      eTtotalRecoBkgd+=myEt;
+                     TParticlePDG *simpdg = simPart->GetPDG(0);
+                     if(!simpdg) continue;
                      if(fUseRecoPt){//Then we switch the pT and the Et
                        myEt = Et(track->P(),track->Theta(),simPart->GetPDG(0)->PdgCode(),track->Charge());
                        pT = track->Pt();
@@ -820,6 +855,8 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                      float pT = simPart->Pt();
                      float eta = simPart->Eta();
                      eTtotalRecoBkgd+=myEt;
+                     TParticlePDG *simpdg = simPart->GetPDG(0);
+                     if(!simpdg) continue;
                      if(fUseRecoPt){//Then we switch the pT and the Et
                        myEt = Et(track->P(),track->Theta(),simPart->GetPDG(0)->PdgCode(),track->Charge());
                        pT = track->Pt();
@@ -831,10 +868,15 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                    }
 
                    if(mom->GetFirstMother()>0){
-                     TParticle *grandma = stack->Particle(mom->GetFirstMother());
+                     TParticle *grandma = NULL;
+                     if(mom->GetFirstMother()>=0) stack->Particle(mom->GetFirstMother());
                      if(grandma){
+                       TParticlePDG *mompdg = mom->GetPDG(0);
+                       if(!mompdg) continue;
                        Int_t pdgCodeMom =  mom->GetPDG(0)->PdgCode();
                        if(pdgCodeMom==fgPiPlusCode || pdgCodeMom==fgPiMinusCode || pdgCodeMom==fgProtonCode ||pdgCodeMom==fgAntiProtonCode || pdgCodeMom==fgKPlusCode || pdgCode==fgKMinusCode){
+                         TParticlePDG *grandmapdg = grandma->GetPDG(0);
+                         if(!grandmapdg) continue;
                          Int_t pdgCodeGrandma =  grandma->GetPDG(0)->PdgCode();
                      
                          if(pdgCodeGrandma == fgXiCode){
@@ -903,6 +945,7 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                    }
                    if(!written){
                      int mycode = simPart->GetPDG(0)->PdgCode();
+                     if(!mycode) continue;
                      if( (pdgCode == fgGammaCode || pdgCode == fgPi0Code) && (mycode==fgEPlusCode||mycode==fgEMinusCode)){
                        written = true;
                        float myEt = Et(simPart);
@@ -1044,7 +1087,8 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
        if (stack->IsPhysicalPrimary(iPart)){//primaries
 
          if (TMath::Abs(part->Eta()) < fHadEtReco->GetCorrections()->GetEtaCut())          {
-
+           TParticlePDG *pdg = part->GetPDG(0);
+           if(!pdg) continue;
            Int_t pdgCode =  part->GetPDG(0)->PdgCode();
            bool filled = false;
            //Investigating smearing...
@@ -1054,6 +1098,8 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
                //To investigate Smearing...
                Float_t myet = Et(part);
                fSimPiKPEt += myet;
+               if(part->Pt()>0.150) fSimRawEtTPC += myet;
+               if(part->Pt()>0.100) fSimRawEtITS += myet;
                Float_t theta = part->Theta();
                Short_t charge = 1;
                Float_t momentum = part->P();
@@ -1590,11 +1636,13 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
              Int_t pdgCodeMom = -99999999;
              float momEta = -30;
              float mompT = -5;
-             if(part->GetFirstMother()){
+             if(part->GetFirstMother()>=0){
                mom = stack->Particle(part->GetFirstMother());
-               pdgCodeMom =  mom->GetPDG(0)->PdgCode();
-               momEta = mom->Eta();
-               mompT = mom->Pt();
+               if(mom->GetPDG(0)){
+                 pdgCodeMom =  mom->GetPDG(0)->PdgCode();
+                 momEta = mom->Eta();
+                 mompT = mom->Pt();
+               }
              }
              //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
              if(pdgCodeMom == fgEtaCode){
@@ -1661,58 +1709,77 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
     FillHisto1D("SimTotEt",fSimTotEt,1.0);
     FillHisto1D("SimHadEt",fSimHadEt,1.0);
     FillHisto1D("SimPiKPEt",fSimPiKPEt,1.0);
+    FillHisto1D("SimRawEtTPC",fSimRawEtTPC,1.0);
+    FillHisto1D("SimRawEtITS",fSimRawEtITS,1.0);
     if(fDataSet!=20100 && AliPWG0Helper::GetEventProcessType(mcEvent->Header()) == AliPWG0Helper::kND){
       FillHisto1D("SimHadEtND",fSimHadEt,1.0);
-      FillHisto1D("SimTotEtND",fSimHadEt,1.0);
+      FillHisto1D("SimTotEtND",fSimTotEt,1.0);
       FillHisto1D("NEventsND",0.5,1);
       FillHisto1D("SimPiKPEtND",fSimPiKPEt,1.0);
+      FillHisto1D("SimRawEtNDTPC",fSimRawEtTPC,1.0);
+      FillHisto1D("SimRawEtNDITS",fSimRawEtITS,1.0);
       if(kIsOfflineV0AND){
        FillHisto1D("SimHadEtNDV0AND",fSimHadEt,1.0);
-       FillHisto1D("SimTotEtNDV0AND",fSimHadEt,1.0);
+       FillHisto1D("SimTotEtNDV0AND",fSimTotEt,1.0);
        FillHisto1D("NEventsNDV0AND",0.5,1);
        FillHisto1D("SimPiKPEtNDV0AND",fSimPiKPEt,1.0);
+       FillHisto1D("SimRawEtNDV0ANDTPC",fSimRawEtTPC,1.0);
+       FillHisto1D("SimRawEtNDV0ANDITS",fSimRawEtITS,1.0);
       }
       if(kIsOfflineMB){
        FillHisto1D("SimHadEtNDMB",fSimHadEt,1.0);
-       FillHisto1D("SimTotEtNDMB",fSimHadEt,1.0);
+       FillHisto1D("SimTotEtNDMB",fSimTotEt,1.0);
        FillHisto1D("NEventsNDMB",0.5,1);
        FillHisto1D("SimPiKPEtNDMB",fSimPiKPEt,1.0);
+       FillHisto1D("SimRawEtNDMBTPC",fSimRawEtTPC,1.0);
+       FillHisto1D("SimRawEtNDMBITS",fSimRawEtITS,1.0);
       }
     }
     if(fDataSet!=20100 && AliPWG0Helper::GetEventProcessType(mcEvent->Header()) == AliPWG0Helper::kSD){
       FillHisto1D("SimHadEtSD",fSimHadEt,1.0);
-      FillHisto1D("SimTotEtSD",fSimHadEt,1.0);
+      FillHisto1D("SimTotEtSD",fSimTotEt,1.0);
       FillHisto1D("NEventsSD",0.5,1);
       FillHisto1D("SimPiKPEtSD",fSimPiKPEt,1.0);
+      FillHisto1D("SimRawEtSDTPC",fSimRawEtTPC,1.0);
+      FillHisto1D("SimRawEtSDITS",fSimRawEtITS,1.0);
       if(kIsOfflineV0AND){
        FillHisto1D("SimHadEtSDV0AND",fSimHadEt,1.0);
-       FillHisto1D("SimTotEtSDV0AND",fSimHadEt,1.0);
+       FillHisto1D("SimTotEtSDV0AND",fSimTotEt,1.0);
        FillHisto1D("NEventsSDV0AND",0.5,1);
        FillHisto1D("SimPiKPEtSDV0AND",fSimPiKPEt,1.0);
+       FillHisto1D("SimRawEtSDV0ANDTPC",fSimRawEtTPC,1.0);
+       FillHisto1D("SimRawEtSDV0ANDITS",fSimRawEtITS,1.0);
       }
       if(kIsOfflineMB){
        FillHisto1D("SimHadEtSDMB",fSimHadEt,1.0);
-       FillHisto1D("SimTotEtSDMB",fSimHadEt,1.0);
+       FillHisto1D("SimTotEtSDMB",fSimTotEt,1.0);
        FillHisto1D("NEventsSDMB",0.5,1);
        FillHisto1D("SimPiKPEtSDMB",fSimPiKPEt,1.0);
+       FillHisto1D("SimRawEtSDMBTPC",fSimRawEtTPC,1.0);
+       FillHisto1D("SimRawEtSDMBITS",fSimRawEtITS,1.0);
       }
     }
     if(fDataSet!=20100 && AliPWG0Helper::GetEventProcessType(mcEvent->Header()) == AliPWG0Helper::kDD){
       FillHisto1D("SimHadEtDD",fSimHadEt,1.0);
-      FillHisto1D("SimTotEtDD",fSimHadEt,1.0);
+      FillHisto1D("SimTotEtDD",fSimTotEt,1.0);
       FillHisto1D("NEventsDD",0.5,1);
       FillHisto1D("SimPiKPEtDD",fSimPiKPEt,1.0);
+      FillHisto1D("SimRawEtDDTPC",fSimRawEtTPC,1.0);
       if(kIsOfflineV0AND){
        FillHisto1D("SimHadEtDDV0AND",fSimHadEt,1.0);
-       FillHisto1D("SimTotEtDDV0AND",fSimHadEt,1.0);
+       FillHisto1D("SimTotEtDDV0AND",fSimTotEt,1.0);
        FillHisto1D("NEventsDDV0AND",0.5,1);
        FillHisto1D("SimPiKPEtDDV0AND",fSimPiKPEt,1.0);
+       FillHisto1D("SimRawEtDDV0ANDTPC",fSimRawEtTPC,1.0);
+       FillHisto1D("SimRawEtDDV0ANDITS",fSimRawEtITS,1.0);
       }
       if(kIsOfflineMB){
        FillHisto1D("SimHadEtDDMB",fSimHadEt,1.0);
-       FillHisto1D("SimTotEtDDMB",fSimHadEt,1.0);
+       FillHisto1D("SimTotEtDDMB",fSimTotEt,1.0);
        FillHisto1D("NEventsDDMB",0.5,1);
        FillHisto1D("SimPiKPEtDDMB",fSimPiKPEt,1.0);
+       FillHisto1D("SimRawEtDDMBTPC",fSimRawEtTPC,1.0);
+       FillHisto1D("SimRawEtDDMBITS",fSimRawEtITS,1.0);
       }
     }
     if(fCentBin != -1){//if we have Pb+Pb and a centrality bin was found
@@ -1747,6 +1814,22 @@ void AliAnalysisHadEtMonteCarlo::Init()
 }
 void AliAnalysisHadEtMonteCarlo::CreateHistograms(){
   //for simulated Et only (no reconstruction)
+
+  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
+  if (!man) {
+    AliFatal("Analysis manager needed");
+    return;
+  }
+  AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
+  if (!inputHandler) {
+    AliFatal("Input handler needed");
+    return;
+  }
+
+  //pid response object
+  fPIDResponse=inputHandler->GetPIDResponse();
+  if (!fPIDResponse) AliError("PIDResponse object was not created");
+
   if( !fRunLightweight){
     CreateEtaPtHisto2D(TString("EtSimulatedPiPlus"),TString("Simulated E_{T} from #pi^{+}"));
     CreateEtaPtHisto2D("EtSimulatedPiMinus","Simulated E_{T} from #pi^{-}");
@@ -2018,8 +2101,15 @@ void AliAnalysisHadEtMonteCarlo::CreateHistograms(){
   Float_t minEt = 0.0;
   Float_t maxEt = 100.0;
   Float_t minEtPiKP = 0.0;
-  Float_t maxEtPiKP = 50.0;
-  if(fDataSet==20100) maxEt=4000.0;
+  Float_t maxEtPiKP = 100.0;
+  if(fDataSet==20100){
+    maxEt=4000.0;
+    maxEtPiKP = 2500;
+  }
+  if(fDataSet==2013){
+    maxEt=100.0;
+    maxEtPiKP = 100.0;
+  }
   Int_t nbinsEt = 100;
   char histoname[200];
   char histotitle[200];
@@ -2035,9 +2125,11 @@ void AliAnalysisHadEtMonteCarlo::CreateHistograms(){
   TString *sHadEt = new TString("HadEt");
   TString *sTotEt = new TString("TotEt");
   TString *sPiKPEt =  new TString("PiKPEt");
+  TString *sRawEt =  new TString("RawEt");
   TString *sTotEtString = new TString("total E_{T}");
   TString *sHadEtString = new TString("hadronic E_{T}");
   TString *sPiKPEtString = new TString("E_{T}^{#pi,K,p}");
+  TString *sRawEtString = new TString("E_{T}^{raw}");
   TString *sFull = new TString("Full");
   TString *sEMCAL = new TString("EMCAL");
   TString *sPHOS = new TString("PHOS");
@@ -2102,6 +2194,13 @@ void AliAnalysisHadEtMonteCarlo::CreateHistograms(){
              snprintf(ytitle,50,"Reconstructed %s",sPiKPEtString->Data());
              snprintf(xtitle,50,"Simulated %s",sPiKPEtString->Data());
              CreateHisto2D(histoname,histotitle,xtitle,ytitle,nbinsEt*4,minEtPiKP,maxEtPiKP,nbinsEt*4,minEtPiKP,maxEtPiKP);
+
+             //And for the raw ET
+             snprintf(histoname,200,"Sim%sVsReco%s%sAcceptance%s%s",sRawEt->Data(),sRawEt->Data(),acceptance->Data(),detector->Data(),partid->Data());
+             snprintf(histotitle,200,"Simulated %s vs reconstructed %s with %s acceptance for p_{T}>%s GeV/c%s",sRawEtString->Data(),sRawEtString->Data(),acceptance->Data(),ptstring->Data(),partidstring->Data());
+             snprintf(ytitle,50,"Reconstructed %s",sRawEtString->Data());
+             snprintf(xtitle,50,"Simulated %s",sRawEtString->Data());
+             CreateHisto2D(histoname,histotitle,xtitle,ytitle,nbinsEt*4,minEtPiKP,maxEtPiKP,nbinsEt*4,minEtPiKP,maxEtPiKP);
            }
          }
 
@@ -2126,6 +2225,26 @@ void AliAnalysisHadEtMonteCarlo::CreateHistograms(){
   CreateHisto1D("SimPiKPEtNDMB","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
   CreateHisto1D("SimPiKPEtDDMB","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
   CreateHisto1D("SimPiKPEtSDMB","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtTPC","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtNDTPC","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtDDTPC","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtSDTPC","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtNDV0ANDTPC","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtDDV0ANDTPC","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtSDV0ANDTPC","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtNDMBTPC","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtDDMBTPC","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtSDMBTPC","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtITS","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtNDITS","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtDDITS","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtSDITS","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtNDV0ANDITS","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtDDV0ANDITS","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtSDV0ANDITS","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtNDMBITS","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtDDMBITS","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
+  CreateHisto1D("SimRawEtSDMBITS","Simulated #pi,K,p E_{T}","Simulated #pi,K,p E_{T}","Number of events",nbinsEt*4,minEtPiKP,maxEtPiKP);
   CreateHisto1D("SimTotEt","Simulated Total E_{T}","Simulated Total E_{T}","Number of events",nbinsEt*4,minEt,maxEt);
   CreateHisto1D("SimHadEt","Simulated Hadronic E_{T}","Simulated Hadronic E_{T}","Number of events",nbinsEt*4,minEt,maxEt);
   CreateHisto1D("SimTotEtND","Simulated Total E_{T}","Simulated Total E_{T} for non-diffractive events","Number of events",nbinsEt*4,minEt,maxEt);