**************************************************************************/
//-----------------------------------------------------------------------
-// Example of task (running locally, on AliEn and CAF),
-// which provides standard way of calculating acceptance and efficiency
-// between different steps of the procedure.
+// Efficiency between different steps of the procedure.
// The ouptut of the task is a AliCFContainer from which the efficiencies
// can be calculated
//-----------------------------------------------------------------------
#include "AliESDtrack.h"
#include "AliESDtrackCuts.h"
#include "AliExternalTrackParam.h"
+#include "AliCentrality.h"
#include "AliLog.h"
fReadAODData(0),
fCFManagerPos(0x0),
fCFManagerNeg(0x0),
- fESD(0),
- fTrackCuts(0),
- fTrackCutsTPConly(0),
+ fESD(0x0),
+ fMC(0x0),
+ fStack(0x0),
+ fVtx(0x0),
+ fIsPbPb(0),
+ fCentClass(10),
+ fTrackType(0),
+ fTrackCuts(0x0),
+ fSigmaConstrainedMax(5.),
fAvgTrials(1),
fHistList(0),
fNEventAll(0),
fNEventSel(0),
fNEventReject(0),
+ fh1Centrality(0x0),
fh1Xsec(0),
fh1Trials(0),
fh1PtHard(0),
fReadAODData(0),
fCFManagerPos(0x0),
fCFManagerNeg(0x0),
- fESD(0),
- fTrackCuts(),
- fTrackCutsTPConly(0),
+ fESD(0x0),
+ fMC(0x0),
+ fStack(0x0),
+ fVtx(0x0),
+ fIsPbPb(0),
+ fCentClass(10),
+ fTrackType(0),
+ fTrackCuts(0x0),
+ fSigmaConstrainedMax(5.),
fAvgTrials(1),
fHistList(0),
fNEventAll(0),
fNEventSel(0),
fNEventReject(0),
+ fh1Centrality(0x0),
fh1Xsec(0),
fh1Trials(0),
fh1PtHard(0),
// Only called once at beginning
//
PostData(3,fTrackCuts);
- PostData(4,fTrackCutsTPConly);
}
//________________________________________________________________________
TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
if (!tree) {
- AliDebug(2,Form("ERROR: Could not read chain from input slot 0"));
- } else {
-
- AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
-
- if (!esdH) {
- AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
- } else {
- fESD = esdH->GetEvent();
- }
+ AliDebug(2,Form( "ERROR: Could not read chain from input slot 0 \n"));
+ return;
}
+
+ AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+
+ if (!esdH) {
+ AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
+ return;
+ } else
+ fESD = esdH->GetEvent();
+ AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+ if (!eventHandler) {
+ AliDebug(2,Form( "ERROR: Could not retrieve MC event handler \n"));
+ }
+ else
+ fMC = eventHandler->MCEvent();
+
}
-//_________________________________________________
-void AliPWG4HighPtSpectra::Exec(Option_t *)
-{
+
+//________________________________________________________________________
+Bool_t AliPWG4HighPtSpectra::SelectEvent() {
//
- // Main loop function
+ // Decide if event should be selected for analysis
//
- AliDebug(2,Form(">> AliPWG4HighPtSpectra::Exec \n"));
- // All events without selection
- fNEventAll->Fill(0.);
+ // Checks following requirements:
+ // - fESD available
+ // - trigger info from AliPhysicsSelection
+ // - number of reconstructed tracks > 1
+ // - primary vertex reconstructed
+ // - z-vertex < 10 cm
+
+ Bool_t selectEvent = kTRUE;
+ //fESD object available?
if (!fESD) {
- AliDebug(2,Form("ERROR: fESD not available"));
+ AliDebug(2,Form("ERROR: fInputEvent not available\n"));
fNEventReject->Fill("noESD",1);
- PostData(0,fHistList);
- PostData(1,fCFManagerPos->GetParticleContainer());
- PostData(2,fCFManagerNeg->GetParticleContainer());
- return;
+ selectEvent = kFALSE;
+ return selectEvent;
}
+ //Trigger
UInt_t isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
if(!(isSelected&AliVEvent::kMB)) { //Select collison candidates
AliDebug(2,Form(" Trigger Selection: event REJECTED ... "));
fNEventReject->Fill("Trigger",1);
+ selectEvent = kFALSE;
+ return selectEvent;
+ }
+
+ //Check if number of reconstructed tracks is larger than 1
+ if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2) {
+ fNEventReject->Fill("NTracks<2",1);
+ selectEvent = kFALSE;
+ return selectEvent;
+ }
+
+ //Check if vertex is reconstructed
+ fVtx = fESD->GetPrimaryVertexSPD();
+
+ if(!fVtx) {
+ fNEventReject->Fill("noVTX",1);
+ selectEvent = kFALSE;
+ return selectEvent;
+ }
+
+ if(!fVtx->GetStatus()) {
+ fNEventReject->Fill("VtxStatus",1);
+ selectEvent = kFALSE;
+ return selectEvent;
+ }
+
+ // Need vertex cut
+ // TString vtxName(fVtx->GetName());
+ if(fVtx->GetNContributors()<2) {
+ fNEventReject->Fill("NCont<2",1);
+ selectEvent = kFALSE;
+ return selectEvent;
+ }
+
+ //Check if z-vertex < 10 cm
+ double primVtx[3];
+ fVtx->GetXYZ(primVtx);
+ if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
+ fNEventReject->Fill("ZVTX>10",1);
+ selectEvent = kFALSE;
+ return selectEvent;
+ }
+
+ //Centrality selection should only be done in case of PbPb
+ if(IsPbPb()) {
+ Float_t cent = 0.;
+ if(fCentClass!=CalculateCentrality(fESD) && fCentClass!=10) {
+ fNEventReject->Fill("cent",1);
+ selectEvent = kFALSE;
+ return selectEvent;
+ }
+ else {
+ if(dynamic_cast<AliESDEvent*>(fESD)->GetCentrality()) {
+ cent = dynamic_cast<AliESDEvent*>(fESD)->GetCentrality()->GetCentralityPercentile("V0M");
+ }
+ if(cent>90.) {
+ fNEventReject->Fill("cent>90",1);
+ selectEvent = kFALSE;
+ return selectEvent;
+ }
+ fh1Centrality->Fill(cent);
+ }
+ }
+
+ return selectEvent;
+
+}
+
+//________________________________________________________________________
+Int_t AliPWG4HighPtSpectra::CalculateCentrality(AliESDEvent *esd){
+
+
+ Float_t cent = 999;
+
+ if(esd){
+ if(esd->GetCentrality()){
+ cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
+ }
+ }
+
+ if(cent<0) return 5;
+ if(cent>80)return 4;
+ if(cent>50)return 3;
+ if(cent>30)return 2;
+ if(cent>10)return 1;
+ return 0;
+
+}
+
+//_________________________________________________
+void AliPWG4HighPtSpectra::Exec(Option_t *)
+{
+ //
+ // Main loop function
+ //
+ AliDebug(2,Form(">> AliPWG4HighPtSpectra::Exec \n"));
+
+ // All events without selection
+ fNEventAll->Fill(0.);
+
+ if(!SelectEvent()) {
+ fNEventReject->Fill("NTracks<2",1);
+ // Post output data
PostData(0,fHistList);
PostData(1,fCFManagerPos->GetParticleContainer());
PostData(2,fCFManagerNeg->GetParticleContainer());
return;
}
- // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
- // This handler can return the current MC event
-
- AliMCEventHandler *eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
-
- AliStack* stack = 0x0;
- AliMCEvent* mcEvent = 0x0;
-
- if(eventHandler) {
- mcEvent = eventHandler->MCEvent();
- if (!mcEvent) {
- AliDebug(2,Form("ERROR: Could not retrieve MC event"));
- fNEventReject->Fill("noMCEvent",1);
- PostData(0,fHistList);
- PostData(1,fCFManagerPos->GetParticleContainer());
- PostData(2,fCFManagerNeg->GetParticleContainer());
- return;
- }
-
- AliDebug(2,Form("MC particles: %d", mcEvent->GetNumberOfTracks()));
-
- stack = mcEvent->Stack(); //Particles Stack
-
- AliDebug(2,Form("MC particles stack: %d", stack->GetNtrack()));
+ //MCEvent available?
+ //if yes: get stack
+ if(fMC) {
+ AliDebug(2,Form("MC particles: %d", fMC->GetNumberOfTracks()));
+ fStack = fMC->Stack(); //Particles Stack
+ AliDebug(2,Form("MC particles stack: %d", fStack->GetNtrack()));
}
- //___ get MC information __________________________________________________________________
-
+ // ---- Get MC Header information (for MC productions in pThard bins) ----
Double_t ptHard = 0.;
Double_t nTrials = 1; // trials for MC trigger weight for real data
- if(mcEvent){
- AliGenPythiaEventHeader* pythiaGenHeader = GetPythiaEventHeader(mcEvent);
+ if(fMC){
+ AliGenPythiaEventHeader* pythiaGenHeader = GetPythiaEventHeader(fMC);
if(pythiaGenHeader){
nTrials = pythiaGenHeader->Trials();
ptHard = pythiaGenHeader->GetPtHard();
}
}
- const AliESDVertex *vtx = fESD->GetPrimaryVertex();
- AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
- // Need vertex cut
- TString vtxName(vtx->GetName());
- if(vtx->GetNContributors() < 2 || (vtxName.Contains("TPCVertex")) ) {
- // SPD vertex
- vtx = fESD->GetPrimaryVertexSPD();
- if(vtx->GetNContributors()<2) {
- vtx = 0x0;
- fNEventReject->Fill("noVTX",1);
- // Post output data
- PostData(0,fHistList);
- PostData(1,fCFManagerPos->GetParticleContainer());
- PostData(2,fCFManagerNeg->GetParticleContainer());
- return;
- }
- }
-
- double primVtx[3];
- vtx->GetXYZ(primVtx);
- if(TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2]>10.)){
- fNEventReject->Fill("ZVTX>10",1);
- PostData(0,fHistList);
- PostData(1,fCFManagerPos->GetParticleContainer());
- PostData(2,fCFManagerNeg->GetParticleContainer());
- return;
- }
-
- if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2){
- fNEventReject->Fill("NTracks<2",1);
- // Post output data
- PostData(0,fHistList);
- PostData(1,fCFManagerPos->GetParticleContainer());
- PostData(2,fCFManagerNeg->GetParticleContainer());
- return;
- }
Int_t nTracks = fESD->GetNumberOfTracks();
AliDebug(2,Form("nTracks %d", nTracks));
Double_t containerInputRec[3] = {0.,0.,0.};
- Double_t containerInputTPConly[3] = {0.,0.,0.};
Double_t containerInputMC[3] = {0.,0.,0.};
- Double_t containerInputRecMC[3] = {0.,0.,0.};
- Double_t containerInputTPConlyMC[3] = {0.,0.,0.};
+ Double_t containerInputRecMC[3] = {0.,0.,0.}; //reconstructed yield as function of MC variable
//Now go to rec level
for (Int_t iTrack = 0; iTrack<nTracks; iTrack++)
{
- if(!fESD->GetTrack(iTrack) ) continue;
- AliESDtrack* track = fESD->GetTrack(iTrack);
- if(!(AliExternalTrackParam *)track->GetTPCInnerParam()) continue;
- AliExternalTrackParam *trackTPC = (AliExternalTrackParam *)track->GetTPCInnerParam();
- if(!track || !trackTPC) continue;
+ //Get track for analysis
+ AliESDtrack *track = 0x0;
+ AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
+ if(!esdtrack) continue;
+
+ if(fTrackType==1)
+ track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
+ else if(fTrackType==2) {
+ track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
+ if(!track) continue;
+
+ AliExternalTrackParam exParam;
+ Bool_t relate = track->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam);
+ if( !relate ) {
+ delete track;
+ continue;
+ }
+ track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
+ }
+ else
+ track = esdtrack;
+
+ if(!track) {
+ if(fTrackType==1 || fTrackType==2) delete track;
+ continue;
+ }
+
+ if(fTrackType==2) {
+ //Cut on chi2 of constrained fit
+ if(track->GetConstrainedChi2TPC() > fSigmaConstrainedMax*fSigmaConstrainedMax) {
+ delete track;
+ continue;
+ }
+ }
+
//fill the container
containerInputRec[0] = track->Pt();
containerInputRec[1] = track->Phi();
containerInputRec[2] = track->Eta();
- //Store TPC Inner Params for TPConly tracks
- containerInputTPConly[0] = trackTPC->Pt();
- containerInputTPConly[1] = trackTPC->Phi();
- containerInputTPConly[2] = trackTPC->Eta();
-
- AliESDtrack* trackTPCESD = fTrackCutsTPConly->GetTPCOnlyTrack(fESD, iTrack);
- if(trackTPCESD) {
- if (fTrackCutsTPConly->AcceptTrack(trackTPCESD)) {
- if(trackTPC->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
- if(trackTPC->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputTPConly,kStepReconstructedTPCOnly);
-
- //Only fill the MC containers if MC information is available
- if(eventHandler) {
- Int_t label = TMath::Abs(track->GetLabel());
- TParticle *particle = stack->Particle(label) ;
- if(!particle) continue;
-
- containerInputTPConlyMC[0] = particle->Pt();
- containerInputTPConlyMC[1] = particle->Phi();
- containerInputTPConlyMC[2] = particle->Eta();
-
- //Container with primaries
- if(stack->IsPhysicalPrimary(label)) {
- if(particle->GetPDG()->Charge()>0.) {
- fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedTPCOnlyMC);
- }
- if(particle->GetPDG()->Charge()<0.) {
- fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedTPCOnlyMC);
- }
- }
- }
-
- }
- }
-
if (fTrackCuts->AcceptTrack(track)) {
if(track->GetSign()>0.) fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
if(track->GetSign()<0.) fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepReconstructed);
//Only fill the MC containers if MC information is available
- if(eventHandler) {
+ if(fMC) {
Int_t label = TMath::Abs(track->GetLabel());
- TParticle *particle = stack->Particle(label) ;
- if(!particle) continue;
-
+ TParticle *particle = fStack->Particle(label) ;
+ if(!particle) {
+ if(fTrackType==1 || fTrackType==2)
+ delete track;
+ continue;
+ }
containerInputRecMC[0] = particle->Pt();
containerInputRecMC[1] = particle->Phi();
containerInputRecMC[2] = particle->Eta();
//Container with primaries
- if(stack->IsPhysicalPrimary(label)) {
+ if(fStack->IsPhysicalPrimary(label)) {
if(particle->GetPDG()->Charge()>0.) {
fCFManagerPos->GetParticleContainer()->Fill(containerInputRecMC,kStepReconstructedMC);
}
}
//Container with secondaries
- if (!stack->IsPhysicalPrimary(label) ) {
+ if (!fStack->IsPhysicalPrimary(label) ) {
if(particle->GetPDG()->Charge()>0.) {
- fCFManagerPos->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
+ fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
}
if(particle->GetPDG()->Charge()<0.) {
- fCFManagerNeg->GetParticleContainer()->Fill(containerInputRec,kStepSecondaries);
+ fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepSecondaries);
}
}
}
- }//trackCuts
+ }//trackCuts global tracks
- delete trackTPCESD;
+ if(fTrackType==1 || fTrackType==2)
+ delete track;
}//track loop
//Fill MC containters if particles are findable
- if(eventHandler) {
- for(int iPart = 1; iPart<(mcEvent->GetNumberOfPrimaries()); iPart++)//stack->GetNprimary();
- {
- AliMCParticle *mcPart = (AliMCParticle*)mcEvent->GetTrack(iPart);
- if(!mcPart) continue;
- //fill the container
- containerInputMC[0] = mcPart->Pt();
- containerInputMC[1] = mcPart->Phi();
- containerInputMC[2] = mcPart->Eta();
-
- if(stack->IsPhysicalPrimary(iPart)) {
- if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
- if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
- }
+ if(fMC) {
+ for(int iPart = 1; iPart<(fMC->GetNumberOfPrimaries()); iPart++) {
+ AliMCParticle *mcPart = (AliMCParticle*)fMC->GetTrack(iPart);
+ if(!mcPart) continue;
+ //fill the container
+ containerInputMC[0] = mcPart->Pt();
+ containerInputMC[1] = mcPart->Phi();
+ containerInputMC[2] = mcPart->Eta();
+
+ if(fStack->IsPhysicalPrimary(iPart)) {
+ if(mcPart->Charge()>0. && fCFManagerPos->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerPos->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
+ if(mcPart->Charge()<0. && fCFManagerNeg->CheckParticleCuts(kStepMCAcceptance,mcPart)) fCFManagerNeg->GetParticleContainer()->Fill(containerInputMC,kStepMCAcceptance);
}
+ }
}
PostData(0,fHistList);
fHistList->Add(fNEventSel);
fNEventReject = new TH1F("fNEventReject","Reason events are rejectected for analysis",20,0,20);
+ //Set labels
+ fNEventReject->Fill("noESD",0);
+ fNEventReject->Fill("Trigger",0);
+ fNEventReject->Fill("NTracks<2",0);
+ fNEventReject->Fill("noVTX",0);
+ fNEventReject->Fill("VtxStatus",0);
+ fNEventReject->Fill("NCont<2",0);
+ fNEventReject->Fill("ZVTX>10",0);
+ fNEventReject->Fill("cent",0);
+ fNEventReject->Fill("cent>90",0);
fHistList->Add(fNEventReject);
+ fh1Centrality = new TH1F("fh1Centrality","fh1Centrality; Centrality %",100,0,100);
+ fHistList->Add(fh1Centrality);
+
fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
fHistList->Add(fh1Xsec);
TH1::AddDirectory(oldStatus);
+ PostData(0,fHistList);
+ PostData(1,fCFManagerPos->GetParticleContainer());
+ PostData(2,fCFManagerNeg->GetParticleContainer());
+
}
#endif
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff