X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWGLF%2FtotEt%2FAliAnalysisHadEtMonteCarlo.cxx;h=eef5cc3228b07d131f475a4a81ba18dc41693afe;hb=4b771b12d37c5e4f1247e00076b5efc780e78f9d;hp=7ec0bb7910d345c02f9d56f5d7be939218aaffc5;hpb=ba222433134f8cd9f328a3318abc971cd8c6a15c;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWGLF/totEt/AliAnalysisHadEtMonteCarlo.cxx b/PWGLF/totEt/AliAnalysisHadEtMonteCarlo.cxx index 7ec0bb7910d..eef5cc3228b 100644 --- a/PWGLF/totEt/AliAnalysisHadEtMonteCarlo.cxx +++ b/PWGLF/totEt/AliAnalysisHadEtMonteCarlo.cxx @@ -29,6 +29,13 @@ #include "AliCentrality.h" #include "AliLog.h" #include "AliPWG0Helper.h" +#include "AliPIDResponse.h" +#include "AliTPCPIDResponse.h" +#include "AliInputEventHandler.h" +#include "AliAnalysisManager.h" +#include "AliGenEventHeader.h" +#include "AliGenCocktailEventHeader.h" +#include "AliGenHijingEventHeader.h" //class AliPWG0Helper; //#include "$ALICE_ROOT/PWG0/AliPWG0Helper.h" @@ -41,7 +48,12 @@ Int_t AliAnalysisHadEtMonteCarlo::fgNumSmearWidths = 4; Float_t AliAnalysisHadEtMonteCarlo::fgSmearWidths[4] = {0.005,0.006,0.007,0.008}; AliAnalysisHadEtMonteCarlo::AliAnalysisHadEtMonteCarlo():AliAnalysisHadEt() + ,checkLabelForHIJING(kFALSE) + ,fNMCProducedMin(0) + ,fNMCProducedMax(0) ,fSimPiKPEt(0) + ,fSimRawEtTPC(0) + ,fSimRawEtITS(0) ,fSimHadEt(0) ,fSimTotEt(0) ,fSimPiKPEtShouldBeReco(0) @@ -74,6 +86,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 @@ -89,12 +103,17 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2) return 0; } AliStack *stack = mcEvent->Stack(); + + if(checkLabelForHIJING) SetGeneratorMinMaxParticles(mcEvent); fCentBin= -1; fGoodEvent = kTRUE;//for p+p collisions if we made it this far we have a good event - if(fDataSet==20100){//If this is Pb+Pb - AliCentrality *centrality = realEvent->GetCentrality(); - if(fNCentBins<21) fCentBin= centrality->GetCentralityClass10(fCentralityMethod); - else{ fCentBin= centrality->GetCentralityClass5(fCentralityMethod);} + if(fDataSet==20100 ||fDataSet==2011 ){//If this is Pb+Pb +// AliCentrality *centrality2 = realEvent->GetCentrality(); +// if(fNCentBins<21) fCentBin= centrality2->GetCentralityClass10(fCentralityMethod); +// else{ fCentBin= centrality2->GetCentralityClass5(fCentralityMethod);} +// cout<<"centrality "<GetCentrality(); + fCentBin = GetCentralityBin(fNCentBins, centrality); if(fCentBin ==-1) fGoodEvent = kFALSE;//but for Pb+Pb events we don't want to count events where we did not find a centrality } AnalyseEvent(ev); @@ -159,6 +178,9 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2) for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) { AliESDtrack *track = dynamic_cast (list->At(iTrack)); + UInt_t label = (UInt_t)TMath::Abs(track->GetLabel()); + //if(checkLabelForHIJING && !IsHIJINGLabel(label,mcEvent,stack) ){cout<<"I am rejecting this particle because it is not HIJING"<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); @@ -202,8 +232,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 +261,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 +324,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 +349,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 +508,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 +541,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 +575,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 +609,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 +643,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 +687,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 +706,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 +717,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 +740,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 +751,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 +775,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 +798,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 +817,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 +834,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 +851,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 +868,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 +881,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 +958,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); @@ -1016,6 +1072,7 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev) // Let's play with the stack! AliStack *stack = mcEvent->Stack(); + Int_t nPrim = stack->GetNtrack(); Float_t fSimPiKPEtPtSmeared = 0; @@ -1035,6 +1092,8 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev) 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 + if(checkLabelForHIJING && !IsHIJINGLabel(iPart,mcEvent,stack) ) continue; + if (!part) { Printf("ERROR: Could not get particle %d", iPart); @@ -1044,7 +1103,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 +1114,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 +1652,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 +1725,77 @@ Int_t AliAnalysisHadEtMonteCarlo::AnalyseEvent(AliVEvent* ev) FillHisto1D("SimTotEt",fSimTotEt,1.0); FillHisto1D("SimHadEt",fSimHadEt,1.0); FillHisto1D("SimPiKPEt",fSimPiKPEt,1.0); - if(fDataSet!=20100 && AliPWG0Helper::GetEventProcessType(mcEvent->Header()) == AliPWG0Helper::kND){ + FillHisto1D("SimRawEtTPC",fSimRawEtTPC,1.0); + FillHisto1D("SimRawEtITS",fSimRawEtITS,1.0); + if((fDataSet!=20100 || fDataSet==2011) && 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){ + if((fDataSet!=20100||fDataSet==2011) && 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){ + if((fDataSet!=20100 || fDataSet==2011) && 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 +1830,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(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^{-}"); @@ -1770,7 +1869,7 @@ void AliAnalysisHadEtMonteCarlo::CreateHistograms(){ CreateEtaPtHisto2D("EtNSimulatedAntiProtonEnhanced","Number of simulated #bar{p}"); } CreateEtaPtHisto2D("EtNSimulatedChargedHadron","Number of simulated charged hadrons"); - if(fDataSet==20100){//If this is Pb+Pb + if(fDataSet==20100 || fDataSet==2011){//If this is Pb+Pb Int_t width = 5; if(fNCentBins<21) width = 10; for(Int_t i=0;iData()),"Reconstructed E_{T} from #bar{p}"); } CreateEtaPtHisto2D(Form("EtNReconstructed%sChargedHadron",cutName->Data()),"Reconstructed E_{T} from charged hadrons"); - if(fDataSet==20100){//If this is Pb+Pb + if(fDataSet==20100||fDataSet==2011){//If this is Pb+Pb Int_t width = 5; if(fNCentBins<21) width = 10; for(Int_t j=0;jData()); 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 +2241,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); @@ -2146,7 +2281,7 @@ void AliAnalysisHadEtMonteCarlo::CreateHistograms(){ CreateHisto1D("SimHadEtDDV0AND","Simulated Hadronic E_{T}","Simulated Hadronic E_{T} for doubly diffractive events","Number of events",nbinsEt*4,minEt,maxEt); CreateHisto1D("SimTotEtDDMB","Simulated Total E_{T}","Simulated Total E_{T} for doubly diffractive events","Number of events",nbinsEt*4,minEt,maxEt); CreateHisto1D("SimHadEtDDMB","Simulated Hadronic E_{T}","Simulated Hadronic E_{T} for doubly diffractive events","Number of events",nbinsEt*4,minEt,maxEt); - if(fDataSet==20100){ + if(fDataSet==20100||fDataSet==2011){ Int_t width = 5; if(fNCentBins<21) width = 10; for(Int_t j=0;jGenEventHeader(); + + AliGenCocktailEventHeader *cocktail = dynamic_cast(eventHeader); + + if(!cocktail) return ; + + TList *genHeaders = cocktail->GetHeaders(); + + Int_t nGenerators = genHeaders->GetEntries(); + //printf("N generators %d \n", nGenerators); + + for(Int_t igen = 0; igen < nGenerators; igen++) + { + AliGenEventHeader * eventHeader2 = (AliGenEventHeader*)genHeaders->At(igen) ; + TString name = eventHeader2->GetName(); + + //printf("Generator %d: Class Name %s, Name %s, title %s \n",igen, eventHeader2->ClassName(), name.Data(), eventHeader2->GetTitle()); + + fNMCProducedMin = fNMCProducedMax; + fNMCProducedMax+= eventHeader2->NProduced(); + + if(name.Contains("Hijing",TString::kIgnoreCase)){ + //cout<<"Found HIJING event and set range "<GenEventHeader(); + AliGenCocktailEventHeader *cocktail = dynamic_cast(eventHeader); + + if(!cocktail) return 0x0 ; + + TList *genHeaders = cocktail->GetHeaders(); + + Int_t nGenerators = genHeaders->GetEntries(); + //printf("N generators %d \n", nGenerators); + + for(Int_t igen = 0; igen < nGenerators; igen++) + { + AliGenEventHeader * eventHeader2 = (AliGenEventHeader*)genHeaders->At(igen) ; + TString name = eventHeader2->GetName(); + + //printf("Generator %d: Class Name %s, Name %s, title %s \n",igen, eventHeader2->ClassName(), name.Data(), eventHeader2->GetTitle()); + + if(name.Contains("Hijing",TString::kIgnoreCase)) return eventHeader2 ; + } + + return 0x0; + +} +Bool_t AliAnalysisHadEtMonteCarlo::IsHIJINGLabel(Int_t label,AliMCEvent *eventMC,AliStack *stack) +{ + + // Find if cluster/track was generated by HIJING + + AliGenHijingEventHeader* hijingHeader = dynamic_cast (GetGenEventHeader(eventMC)); + + //printf("header %p, label %d\n",hijingHeader,label); + + if(!hijingHeader || label < 0 ) return kFALSE; + + + //printf("pass a), N produced %d\n",nproduced); + + if(label >= fNMCProducedMin && label < fNMCProducedMax) + { + //printf(" accept!, label is smaller than produced, N %d\n",nproduced); + + return kTRUE; + } + + if(!stack) return kFALSE; + + Int_t nprimaries = stack->GetNtrack(); + + if(label > nprimaries) return kFALSE; + + TParticle * mom = stack->Particle(label); + + Int_t iMom = label; + Int_t iParent = mom->GetFirstMother(); + while(iParent!=-1){ + if(iParent >= fNMCProducedMin && iParent < fNMCProducedMax){ + return kTRUE; + } + + iMom = iParent; + mom = stack->Particle(iMom); + iParent = mom->GetFirstMother(); + } + + return kFALSE ; + +} + + + +