2 #include "AliMCEvent.h"
4 #include "AliVParticle.h"
5 #include "AliESDEvent.h"
6 #include "AliAODEvent.h"
7 #include "AliVVertex.h"
9 #include "AliGenEventHeader.h"
10 #include "AliInputEventHandler.h"
11 #include "AliAnalysisManager.h"
12 #include "AliAODMCHeader.h"
13 #include "AliAODMCParticle.h"
14 #include "AliMCParticle.h"
15 #include "AliVParticle.h"
17 #include "AliMCEventHandler.h"
18 #include "AliPIDResponse.h"
20 #include "AliSingleTrackEffCuts.h"
25 ClassImp(AliSingleTrackEffCuts)
27 //_____________________________________________________
28 AliSingleTrackEffCuts::AliSingleTrackEffCuts():
40 fTriggerMask(AliVEvent::kAny),
50 fParticlePid(AliPID::kPion),
65 // Default constructor
69 //________________________________________________________________________________
70 AliSingleTrackEffCuts::AliSingleTrackEffCuts(const char* name, const char* title):
71 AliAnalysisCuts(name,title),
82 fTriggerMask(AliVEvent::kAny),
92 fParticlePid(AliPID::kPion),
107 // Default constructor
111 //_________________________________________________________________________________
112 AliSingleTrackEffCuts::AliSingleTrackEffCuts(const AliSingleTrackEffCuts &source):
113 AliAnalysisCuts(source),
114 fisAOD(source.fisAOD),
115 fIsPdgCode(source.fIsPdgCode),
116 fPdgCode(source.fPdgCode),
117 fEtaMin(source.fEtaMin),
118 fEtaMax(source.fEtaMax),
121 fPtMin(source.fPtMin),
122 fPtMax(source.fPtMax),
123 fIsCharged(source.fIsCharged),
124 fTriggerMask(source.fTriggerMask),
125 fMinVtxType(source.fMinVtxType),
126 fMinVtxContr(source.fMinVtxContr),
127 fMaxVtxZ(source.fMaxVtxZ),
128 fCutOnZVertexSPD(source.fCutOnZVertexSPD),
129 fnClusITS(source.fnClusITS),
130 fnClusTPC(source.fnClusTPC),
131 fnClusTOF(source.fnClusTOF),
132 fnClusMUON(source.fnClusMUON),
133 fusePid(source.fusePid),
134 fParticlePid(source.fParticlePid),
135 fuseTPCPid(source.fuseTPCPid),
136 fnPTPCBins(source.fnPTPCBins),
137 fnPTPCBinLimits(source.fnPTPCBinLimits),
138 fPTPCBinLimits(source.fPTPCBinLimits),
139 fnSigmaTPC(source.fnSigmaTPC),
140 fPmaxTPC(source.fPmaxTPC),
141 fuseTOFPid(source.fuseTOFPid),
142 fnPTOFBins(source.fnPTOFBins),
143 fnPTOFBinLimits(source.fnPTOFBinLimits),
144 fPTOFBinLimits(source.fPTOFBinLimits),
145 fnSigmaTOF(source.fnSigmaTOF),
146 fPmaxTOF(source.fPmaxTOF)
154 //_________________________________________________________________________________________
155 AliSingleTrackEffCuts &AliSingleTrackEffCuts::operator=(const AliSingleTrackEffCuts &source)
158 // assignment operator
160 if(&source == this) return *this;
162 fisAOD = source.fisAOD;
164 fIsPdgCode = source.fIsPdgCode;
165 fPdgCode = source.fPdgCode;
166 fEtaMin = source.fEtaMin;
167 fEtaMax = source.fEtaMax;
168 fYMin = source.fYMin;
169 fYMax = source.fYMax;
170 fPtMin = source.fPtMin;
171 fPtMax = source.fPtMax;
172 fIsCharged = source.fIsCharged;
174 fTriggerMask = source.fTriggerMask;
175 fMinVtxType = source.fMinVtxType;
176 fMinVtxContr = source.fMinVtxContr;
177 fMaxVtxZ = source.fMaxVtxZ;
178 fCutOnZVertexSPD = source.fCutOnZVertexSPD;
180 fnClusITS = source.fnClusITS;
181 fnClusTPC = source.fnClusTPC;
182 fnClusTOF = source.fnClusTOF;
183 fnClusMUON = source.fnClusMUON;
185 fusePid = source.fusePid;
186 fParticlePid = source.fParticlePid;
187 fuseTPCPid = source.fuseTPCPid;
188 fnPTPCBins = source.fnPTPCBins;
189 fnPTPCBinLimits = source.fnPTPCBinLimits;
190 if(source.fPTPCBinLimits && source.fnSigmaTPC)
191 SetTPCSigmaPtBins(source.fnPTPCBins,source.fPTPCBinLimits,source.fnSigmaTPC);
192 fPmaxTPC = source.fPmaxTPC;
193 fuseTOFPid = source.fuseTOFPid;
194 fnPTOFBins = source.fnPTOFBins;
195 fnPTOFBinLimits = source.fnPTOFBinLimits;
196 if(source.fPTOFBinLimits && source.fnSigmaTOF)
197 SetTOFSigmaPtBins(source.fnPTOFBins,source.fPTOFBinLimits,source.fnSigmaTOF);
198 fPmaxTOF = source.fPmaxTOF;
203 //______________________________________________
204 AliSingleTrackEffCuts::~AliSingleTrackEffCuts()
211 delete [] fPTPCBinLimits;
215 delete [] fnSigmaTPC;
219 delete [] fPTOFBinLimits;
223 delete [] fnSigmaTOF;
229 //______________________________________________
230 Bool_t AliSingleTrackEffCuts::IsMCEventSelected(TObject* obj)
233 // Event selection at MC level
236 if(!obj) return kFALSE;
237 Bool_t isSelected = kTRUE;
240 AliGenEventHeader *genHeader=0; //ESDs
241 AliAODMCHeader *mcHeader=0; // AODs
242 Bool_t isAOD = obj->IsA()->InheritsFrom("AliAODEvent");
245 event = dynamic_cast<AliMCEvent*>(obj);
246 if (!event) return kFALSE;
247 genHeader = event->GenEventHeader();
248 if (!genHeader) return kFALSE;
250 AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*> (obj);
251 if (!aodEvent) return kFALSE;
252 mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
254 AliError("Could not find MC Header in AOD");
258 // Check for not null trigger mask to evict pPb MC buggy-events
259 Int_t runnumber = aodEvent->GetRunNumber();
260 if(aodEvent->GetTriggerMask()==0 &&
261 (runnumber>=195344 && runnumber<=195677)){
262 AliDebug(3,"Event rejected because of null trigger mask");
268 // Cut on Z Vertex ( |z-vtx| <= fMaxVtxZ )
270 Double_t zMCVertex=-1000;
272 genHeader->PrimaryVertex(vtxPos);
273 zMCVertex = vtxPos[2];
275 zMCVertex = mcHeader->GetVtxZ();
277 if( TMath::Abs(zMCVertex)>fMaxVtxZ ) isSelected = kFALSE;
283 //__________________________________________________________
284 Bool_t AliSingleTrackEffCuts::IsMCParticleGenerated(TObject* obj)
287 // Check generated particles (primary, charged, pdgcode)
290 if(!obj) return kFALSE;
291 if(!obj->InheritsFrom("AliVParticle")) {
292 AliError("object must derived from AliVParticle !");
295 AliVParticle* particle = dynamic_cast<AliVParticle *>(obj);
296 if(!particle) return kFALSE;
298 Bool_t isSelected = kTRUE;
300 // Check particle Pdg Code
301 if(fIsPdgCode && TMath::Abs( particle->PdgCode() )!= fPdgCode) isSelected = kFALSE;
304 if(fIsCharged && (particle->Charge()==0)) isSelected = kFALSE;
306 // Selection of Physical Primary particles
307 if(!fisAOD) { // check on ESDs
308 AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
309 AliMCEvent* mcevent = mcinfo->MCEvent();
310 AliStack* stack = ((AliMCEvent*)mcevent)->Stack();
311 AliMCParticle* mcPart = dynamic_cast<AliMCParticle *>(obj);
312 if (!mcPart) return kFALSE;
313 if(!stack->IsPhysicalPrimary(mcPart->GetLabel())) {
316 } else { // Check on AODs
317 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle *>(obj);
318 if (!mcPart) return kFALSE;
319 if(!mcPart->IsPhysicalPrimary()) isSelected = kFALSE;
326 //_____________________________________________________________________
327 Bool_t AliSingleTrackEffCuts::IsMCParticleInKineAcceptance(TObject *obj)
330 // Check generated particles (eta, y, pt)
333 if(!obj) return kFALSE;
334 if(!obj->InheritsFrom("AliVParticle")) AliError("object must derived from AliVParticle !");
335 AliVParticle* particle = dynamic_cast<AliVParticle *>(obj);
336 if(!particle) return kFALSE;
338 Bool_t isSelected = kTRUE;
341 if(particle->Eta()<fEtaMin || particle->Eta()>fEtaMax) isSelected = kFALSE;
344 Double_t energy = particle->E();
345 Double_t pz = particle->Pz();
346 Double_t particleY = (TMath::Abs(energy-pz)>1e-10) ? 0.5*TMath::Log( (energy+pz)/(energy-pz) ) : 1e6;
347 if(particleY<fYMin || particleY>fYMax) isSelected = kFALSE;
350 if(particle->Pt()<fPtMin || particle->Pt()>fPtMax) isSelected = kFALSE;
356 //_______________________________________________________________________
357 Bool_t AliSingleTrackEffCuts::IsMCParticleInReconstructable(TObject *obj)
360 // Check if particle has left enough hits in the detectors (only at ESD level)
363 if(!obj) return kFALSE;
364 TString className(obj->ClassName());
365 if (className.CompareTo("AliMCParticle") != 0) {
366 AliError("obj must point to an AliMCParticle !");
367 return kTRUE; // <===================================== FIX ME !!
370 AliMCParticle * part = dynamic_cast<AliMCParticle*>(obj);
371 if(!part) return kFALSE;
373 Bool_t isSelected = kTRUE;
375 Int_t nHitsITS=0, nHitsTPC=0, nHitsTRD=0, nHitsTOF=0, nHitsMUON=0;
376 for(Int_t iTrackRef=0; iTrackRef<part->GetNumberOfTrackReferences(); iTrackRef++) {
377 AliTrackReference * trackRef = part->GetTrackReference(iTrackRef);
379 Int_t detectorId = trackRef->DetectorId();
381 case AliTrackReference::kITS :
384 case AliTrackReference::kTPC :
387 case AliTrackReference::kTRD :
390 case AliTrackReference::kTOF :
393 case AliTrackReference::kMUON :
401 if(nHitsITS<fnClusITS) isSelected = kFALSE;
402 if(nHitsTPC<fnClusTPC) isSelected = kFALSE;
403 if(nHitsTOF<fnClusTOF) isSelected = kFALSE;
404 if(nHitsMUON<fnClusMUON) isSelected = kFALSE;
410 //_____________________________________________________________
411 Bool_t AliSingleTrackEffCuts::IsRecoEventSelected(TObject* obj)
414 // Event selection at reconstructed level (trigger, zvtx)
417 AliVEvent *event = dynamic_cast<AliVEvent*>(obj);
418 if(!event) return kFALSE;
420 Bool_t isSelected = kTRUE;
422 // Check if event is accepted by the Physics selection
423 UInt_t trigFired = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
424 Bool_t isEvtSelected = (trigFired & fTriggerMask);
425 if(!isEvtSelected) isSelected = kFALSE;
428 Bool_t isVtxSelected = IsVertexSelected(event);
429 if(!isVtxSelected) isSelected = kFALSE;
435 //______________________________________________________________________
436 Bool_t AliSingleTrackEffCuts::IsRecoParticleKineAcceptance(TObject *obj)
439 // Check if reconstructed particle is in the acceptance (eta, y, pt)
442 Bool_t isSelected = kTRUE;
444 AliVParticle *track = dynamic_cast<AliVParticle*>(obj);
445 if (!track) return kFALSE;
448 if(track->Eta()<fEtaMin || track->Eta()>fEtaMax) isSelected = kFALSE;
451 if(track->Y()<fYMin || track->Y()>fYMax) isSelected = kFALSE;
454 if(track->Pt()<fPtMin || track->Pt() >fPtMax) isSelected = kFALSE;
459 //_______________________________________________________________
460 Bool_t AliSingleTrackEffCuts::IsVertexSelected(AliVEvent *event)
463 // Check if the reconstructed vertex is selected
466 Bool_t accept = kTRUE;
467 Bool_t isAOD = event->IsA()->InheritsFrom("AliAODEvent");
469 const AliVVertex *vertex = event->GetPrimaryVertex();
476 // Cut on vertex type
477 TString title=vertex->GetTitle();
478 if(title.Contains("Z") && fMinVtxType>1){
480 } else if(title.Contains("3D") && fMinVtxType>2){
484 // cut on minimum number of contributors
485 if(vertex->GetNContributors()<fMinVtxContr){
486 AliInfo(Form("too few contributors %d",vertex->GetNContributors()));
490 // cut on absolute |z| of the vertex
491 if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) {
492 AliInfo("outside the Vtx range");
496 // cut on distance of SPD and TRK vertexes
497 if(fCutOnZVertexSPD==1) {
498 const AliVVertex *vSPD = NULL;
500 vSPD = ((AliAODEvent*)event)->GetPrimaryVertexSPD();
502 vSPD = ((AliESDEvent*)event)->GetPrimaryVertexSPD();
505 if(vSPD && vSPD->GetNContributors()>=fMinVtxContr) {
506 if(TMath::Abs(vSPD->GetZ()-vertex->GetZ())>0.5) accept = kFALSE;
513 //_________________________________________________________________________________________________
514 void AliSingleTrackEffCuts::SetTPCSigmaPtBins(Int_t nPtBins, Float_t *pBinLimits, Float_t *sigmaBin)
517 // Set TPC Pid number-of-Sigma in P bins
522 delete [] fPTPCBinLimits;
523 fPTPCBinLimits = NULL;
524 printf("Changing the TPC cut P bins\n");
527 delete [] fnSigmaTPC;
529 printf("Changing the TPC Sigma cut per p bin\n");
532 fnPTPCBins = nPtBins;
533 fnPTPCBinLimits = nPtBins+1;
534 fPTPCBinLimits = new Float_t[fnPTPCBinLimits];
535 for(Int_t ib=0; ib<nPtBins+1; ib++) fPTPCBinLimits[ib]=pBinLimits[ib];
536 fnSigmaTPC = new Float_t[fnPTPCBins];
537 for(Int_t ib=0; ib<nPtBins; ib++) fnSigmaTPC[ib]=sigmaBin[ib];
541 //_____________________________________________________________________________________________________
542 void AliSingleTrackEffCuts::SetTOFSigmaPtBins(Int_t nPtBins, Float_t *pBinLimits, Float_t *sigmaBin)
545 // Set TOF Pid number-of-Sigma in P bins
550 delete [] fPTOFBinLimits;
551 fPTOFBinLimits = NULL;
552 printf("Changing the TOF cut P bins\n");
555 delete [] fnSigmaTOF;
557 printf("Changing the TOF Sigma cut per p bin\n");
560 fnPTOFBins = nPtBins;
561 fnPTOFBinLimits = nPtBins+1;
562 fPTOFBinLimits = new Float_t[fnPTOFBinLimits];
563 for(Int_t ib=0; ib<nPtBins+1; ib++) fPTOFBinLimits[ib]=pBinLimits[ib];
564 fnSigmaTOF = new Float_t[fnPTOFBins];
565 for(Int_t ib=0; ib<nPtBins; ib++) fnSigmaTOF[ib]=sigmaBin[ib];
570 //___________________________________________________________________________
571 Bool_t AliSingleTrackEffCuts::CheckTPCPIDStatus(AliAODTrack *track) const{
573 // Check TPC PID status
576 if ((track->GetStatus() & AliESDtrack::kTPCin)==0) return kFALSE;
577 UShort_t nTPCClus=track->GetTPCClusterMap().CountBits();
578 if (nTPCClus<70) return kFALSE;
582 //___________________________________________________________________________
583 Bool_t AliSingleTrackEffCuts::CheckTOFPIDStatus(AliAODTrack *track) const{
585 // Check TOC PID status
588 if ((track->GetStatus()&AliESDtrack::kTOFout)==0) return kFALSE;
589 if ((track->GetStatus()&AliESDtrack::kTIME)==0) return kFALSE;
590 if ((track->GetStatus()&AliESDtrack::kTOFpid)==0) return kFALSE;
591 if ((track->GetStatus()&AliESDtrack::kTOFmismatch)!=0) return kFALSE;
595 //___________________________________________________________________________
596 Bool_t AliSingleTrackEffCuts::IsRecoParticlePID(TObject *obj)
599 // Check Particle PID (AOD only for now!)
602 Bool_t isSelected = kFALSE;
603 Bool_t isAOD = obj->IsA()->InheritsFrom("AliAODTrack");
605 if(!isAOD || !GetUsePid()) return isSelected;
606 if(!fuseTPCPid && !fuseTOFPid) return isSelected;
608 AliAODTrack *track = dynamic_cast<AliAODTrack*>(obj);
609 if(!track) { cout<<"No track found"<<endl; return isSelected; }
611 // AliAODPid *pid = track->GetDetPid();
612 // if(!pid) { cout<<"No AliAODPid found"<<endl; return isSelected; }
614 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
615 AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
616 AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
617 if(!pidResp) { cout<<"No PidResponse found"<<endl; return isSelected;}
619 Double_t pPart = track->P();
621 // Check detector status
622 Bool_t okTPC = CheckTPCPIDStatus(track);
623 Bool_t okTOF = CheckTOFPIDStatus(track);
625 // Check Number of Sigmas
626 Double_t nsigmaTPC=pidResp->NumberOfSigmasTPC((AliVParticle*)track,(AliPID::EParticleType)fParticlePid);
627 Double_t nsigmaTOF=pidResp->NumberOfSigmasTOF((AliVParticle*)track,(AliPID::EParticleType)fParticlePid);
628 Bool_t isTPCPid=false, isTOFPid=false;
630 // If use TPC and TPC infos are ok, check whether the sigma is ok in the given p range
631 if(fuseTPCPid && okTPC) {
632 for(Int_t j=0; j<fnPTPCBins; j++) {
633 // cout<<" checking bin: ("<<fPTPCBinLimits[j]<<","<<fPTPCBinLimits[j+1]<<") should be nsigma < "<<fnSigmaTPC[j]<<endl;
634 if ((pPart>fPTPCBinLimits[j]) && (pPart<fPTPCBinLimits[j+1]) && nsigmaTPC<fnSigmaTPC[j]) isTPCPid=true;
636 if(pPart>fPmaxTPC) isTPCPid=true;
639 // If use TPC and TPC infos are ok, check whether the sigma is ok in the given p range
640 if(fuseTOFPid && okTOF) {
641 for(Int_t j=0; j<fnPTOFBins; j++) {
642 if ((pPart>fPTOFBinLimits[j]) && (pPart<fPTOFBinLimits[j+1]) && nsigmaTOF<fnSigmaTOF[j]) isTPCPid=true;
644 if(pPart>fPmaxTOF) isTOFPid=true;
647 isSelected = (isTPCPid || isTOFPid) ? true : false;