: AliAnalysisTask("AliPWG4HighPtQAMC", ""),
fESD(0),
fMC(0),
+ fStack(0),
+ fVtx(0x0),
fTrackCuts(0),
fTrackCutsITS(0),
fTrackType(0),
fh1PtHardTrials(0),
fPtAll(0),
fPtSel(0),
+ fPtSelFakes(0),
+ fNPointTPCFakes(0),
+ fPtSelLargeLabel(0),
+ fMultRec(0),
+ fNPointTPCMultRec(0),
+ fDeltaPtMultRec(0),
fPtAllminPtMCvsPtAll(0),
fPtAllminPtMCvsPtAllNPointTPC(0),
fPtAllminPtMCvsPtAllDCAR(0),
AliAnalysisTask(name,""),
fESD(0),
fMC(0),
+ fStack(0),
+ fVtx(0x0),
fTrackCuts(),
fTrackCutsITS(),
fTrackType(0),
fh1PtHardTrials(0),
fPtAll(0),
fPtSel(0),
+ fPtSelFakes(0),
+ fNPointTPCFakes(0),
+ fPtSelLargeLabel(0),
+ fMultRec(0),
+ fNPointTPCMultRec(0),
+ fDeltaPtMultRec(0),
fPtAllminPtMCvsPtAll(0),
fPtAllminPtMCvsPtAllNPointTPC(0),
fPtAllminPtMCvsPtAllDCAR(0),
} 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"));
- }
+ 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();
fNEventSel = new TH1F("fNEventSel","NEvent Selected for analysis",1,-0.5,0.5);
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("noMCEvent",0);
+ fNEventReject->Fill("NTracks<2",0);
+ fNEventReject->Fill("noVTX",0);
+ fNEventReject->Fill("VtxStatus",0);
+ fNEventReject->Fill("NCont<2",0);
+ fNEventReject->Fill("ZVTX>10",0);
fHistList->Add(fNEventReject);
fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
fHistList->Add(fPtAll);
fPtSel = new TH1F("fPtSel","PtSel",fgkNPtBins, fgkPtMin, fgkPtMax);
fHistList->Add(fPtSel);
-
+
+ fPtSelFakes = new TH1F("fPtSelFakes","PtSelFakes",fgkNPtBins, fgkPtMin, fgkPtMax);
+ fHistList->Add(fPtSelFakes);
+
+ fNPointTPCFakes = new TH1F("fNPointTPCFakes","fNPointTPCFakes",160,0.5,160.5);
+ fHistList->Add(fNPointTPCFakes);
+
+ fPtSelLargeLabel = new TH1F("fPtSelLargeLabel","PtSelLargeLabel",fgkNPtBins, fgkPtMin, fgkPtMax);
+ fHistList->Add(fPtSelLargeLabel);
+
+ fMultRec = new TH1F("fMultRec","Multiple reconstruction of tracks",20,0,20);
+ fHistList->Add(fMultRec);
+
+ fNPointTPCMultRec = new TH1F("fNPointTPCMultRec","Multiple reconstruction of tracks",160,0.5,160.5);
+ fHistList->Add(fNPointTPCMultRec);
+
+ fDeltaPtMultRec = new TH2F("fDeltaPtMultRec","Delta pT vs pT for multiple reconstructed tracks",100,0.,50.,100,-20.,20.);
+ fHistList->Add(fDeltaPtMultRec);
+
fPtAllminPtMCvsPtAll = new TH2F("fPtAllminPtMCvsPtAll","PtAllminPtMCvsPtAll",fgkNPtBins, fgkPtMin,fgkPtMax,fgkResPtBins,-1.,1.);
fPtAllminPtMCvsPtAll->SetXTitle("p_{t}^{MC}");
fPtAllminPtMCvsPtAll->SetYTitle("(1/p_{t}^{All}-1/p_{t}^{MC})/(1/p_{t}^{MC})");
TH1::AddDirectory(oldStatus);
}
+
//________________________________________________________________________
-void AliPWG4HighPtQAMC::Exec(Option_t *) {
- // Main loop
- // Called for each event
- AliDebug(2,Form(">> AliPWG4HighPtQATPConly::Exec \n"));
- // All events without selection
- fNEventAll->Fill(0.);
+Bool_t AliPWG4HighPtQAMC::SelectEvent() {
+ //
+ // Decide if event should be selected for analysis
+ //
+ // Checks following requirements:
+ // - fESD available
+ // - trigger info from AliPhysicsSelection
+ // - MCevent available
+ // - 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: fInputEvent not available\n"));
fNEventReject->Fill("noESD",1);
- PostData(0, fHistList);
- PostData(1, fHistListITS);
- 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);
- // Post output data
- PostData(0, fHistList);
- PostData(1, fHistListITS);
- return;
+ selectEvent = kFALSE;
+ return selectEvent;
}
-
- AliStack* stack = 0x0;
-
+
+ //MCEvent available?
+ //if yes: get stack
if(fMC) {
AliDebug(2,Form("MC particles: %d", fMC->GetNumberOfTracks()));
- stack = fMC->Stack(); //Particles Stack
- AliDebug(2,Form("MC particles stack: %d", stack->GetNtrack()));
+ fStack = fMC->Stack(); //Particles Stack
+ AliDebug(2,Form("MC particles stack: %d", fStack->GetNtrack()));
} else {
AliDebug(2,Form("ERROR: Could not retrieve MC eventHandler"));
fNEventReject->Fill("noMCEvent",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
+ if(fTrackType==1) fVtx = fESD->GetPrimaryVertexTPC();
+ else 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;
+ }
+
+ AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",fVtx->GetTitle(), fVtx->GetStatus(), fVtx->GetNContributors()));
+
+ return selectEvent;
+
+}
+
+//________________________________________________________________________
+void AliPWG4HighPtQAMC::Exec(Option_t *) {
+ // Main loop
+ // Called for each event
+ AliDebug(2,Form(">> AliPWG4HighPtQATPConly::Exec \n"));
+ // All events without selection
+ fNEventAll->Fill(0.);
+
+ if(!SelectEvent()) {
+ // Post output data
PostData(0, fHistList);
PostData(1, fHistListITS);
return;
}
- //___ 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
}
}
- const AliESDVertex *vtx = fESD->GetPrimaryVertex();
- // 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, fHistListITS);
- 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);
- // Post output data
- PostData(0, fHistList);
- PostData(1, fHistListITS);
- return;
- }
-
- AliDebug(2,Form("Vertex title %s, status %d, nCont %d\n",vtx->GetTitle(), vtx->GetStatus(), vtx->GetNContributors()));
-
- if(!fESD->GetNumberOfTracks() || fESD->GetNumberOfTracks()<2) {
- fNEventReject->Fill("NTracks<2",1);
- PostData(0, fHistList);
- PostData(1, fHistListITS);
- return;
- }
-
//Need to keep track of selected events
fNEventSel->Fill(0.);
Int_t nTracks = fESD->GetNumberOfTracks();
AliDebug(2,Form("nTracks ESD%d", nTracks));
- int nMCtracks = stack->GetNtrack();
+ int nMCtracks = fStack->GetNtrack();
Float_t pt = 0.;
Float_t ptMC = 0.;
AliESDtrack *track;
AliESDtrack *esdtrack = fESD->GetTrack(iTrack);
if(!esdtrack) continue;
- AliExternalTrackParam *trackTPC = (AliExternalTrackParam *)esdtrack->GetTPCInnerParam();
if(fTrackType==1)
track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
+ else if(fTrackType==2) {
+ track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack->GetID());
+ if(!track) {
+ delete 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) continue;
+ if(!track) {
+ if(fTrackType==1 || fTrackType==2) delete track;
+ continue;
+ }
Int_t label = TMath::Abs(track->GetLabel());
- if(label>=nMCtracks)continue;
- TParticle *particle = stack->Particle(label) ;
- if(!particle) continue;
+ if(label>=nMCtracks) {
+ if (fTrackCuts->AcceptTrack(track)) { fPtSelLargeLabel->Fill(pt); }
+ if(fTrackType==1 || fTrackType==2) delete track;
+ continue;
+ }
+ TParticle *particle = fStack->Particle(label) ;
+ if(!particle) {
+ if(fTrackType==1 || fTrackType==2) delete track;
+ continue;
+ }
ptMC = particle->Pt();
phi = track->Phi();
track->GetImpactParameters(dca2D,dcaZ);
}
- else if(fTrackType==1) { //TPConly
- pt = trackTPC->Pt();
- phi = trackTPC->Phi();
+ else if(fTrackType==1 || fTrackType==2) { //TPConly
+ pt = track->Pt();
+ phi = track->Phi();
track->GetImpactParametersTPC(dca2D,dcaZ);
}
else {continue;}
if (fTrackCuts->AcceptTrack(track)) {
fPtSel->Fill(pt);
-
+ if(track->GetLabel()<0) {
+ fPtSelFakes->Fill(pt);
+ fNPointTPCFakes->Fill(track->GetTPCNcls());
+ }
fPtSelMC->Fill(ptMC);
fPtAllminPtMCvsPtAll->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC) );
fPtAllminPtMCvsPtAllNPointTPC->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),track->GetTPCNcls());
fPtAllminPtMCvsPtAllNSigmaToVertex->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),nSigmaToVertex);
fPtAllminPtMCvsPtAllChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
fPtAllminPtMCvsPtAllRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
- }//fTrackCuts selection
-
-
+
+ //Check if track is reconstructed multiple times
+ int multCounter = 1;
+ for (Int_t iTrack2 = iTrack+1; iTrack2 < nTracks; iTrack2++) {
+ // AliESDtrack *track2 = GetTrackForAnalysis(iTrack2);
+ AliESDtrack *track2;
+ AliESDtrack *esdtrack2 = fESD->GetTrack(iTrack2);
+ if(!esdtrack2) continue;
+
+ if(fTrackType==1)
+ track2 = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack2->GetID());
+ else if(fTrackType==2) {
+ track2 = AliESDtrackCuts::GetTPCOnlyTrack(fESD,esdtrack2->GetID());
+ if(!track2) {
+ delete track2;
+ continue;
+ }
+ AliExternalTrackParam exParam2;
+ Bool_t relate = track2->RelateToVertexTPC(fVtx,fESD->GetMagneticField(),kVeryBig,&exParam2);
+ if( !relate ) {
+ delete track2;
+ continue;
+ }
+ track2->Set(exParam2.GetX(),exParam2.GetAlpha(),exParam2.GetParameter(),exParam2.GetCovariance());
+ }
+ else
+ track2 = esdtrack2;
+ if(!track2) {
+ if(fTrackType==1 || fTrackType==2) delete track2;
+ continue;
+ }
+
+ if (fTrackCuts->AcceptTrack(track2)) {
+ Int_t label2 = TMath::Abs(track2->GetLabel());
+ if(label==label2) {
+ fNPointTPCMultRec->Fill(track->GetTPCNcls());
+ fNPointTPCMultRec->Fill(track2->GetTPCNcls());
+ fDeltaPtMultRec->Fill(track->Pt(),track->Pt()-track2->Pt());
+ multCounter++;
+ }
+ }
+ if(fTrackType==1 || fTrackType==2) delete track2;
+ }//track2 loop
+ if(multCounter>1) fMultRec->Fill(multCounter);
+
+ }//fTrackCuts selection
+
//ITSrefit selection
if (fTrackCutsITS->AcceptTrack(track) && fTrackType==0) {
fPtITSminPtMCvsPtITSChi2C->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),chi2C);
fPtITSminPtMCvsPtITSRel1PtUncertainty->Fill(ptMC,(1./pt-1./ptMC)/(1./ptMC),relUncertainty1Pt);
}//fTrackCutsITS loop
-
+
+ if(fTrackType==1 || fTrackType==2) delete track;
+
}//ESD track loop
// Post output data