1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 // Extra cuts implemented by the ALICE Heavy Flavour Electron Group
19 // - TPC cluster ratio
23 // Markus Fasel <M.Fasel@gsi.de>
33 #include "AliAODEvent.h"
34 #include "AliAODTrack.h"
35 #include "AliAODPid.h"
36 #include "AliAODVertex.h"
37 #include "AliAODTrack.h"
38 #include "AliESDEvent.h"
39 #include "AliESDVertex.h"
40 #include "AliESDtrack.h"
42 #include "AliMCParticle.h"
43 #include "AliVEvent.h"
44 #include "AliVTrack.h"
45 #include "AliVParticle.h"
46 #include "AliVertexerTracks.h"
47 #include "AliVVertex.h"
48 #include "AliExternalTrackParam.h"
50 #include "AliHFEextraCuts.h"
52 ClassImp(AliHFEextraCuts)
54 const Int_t AliHFEextraCuts::fgkNQAhistos = 10;
56 //______________________________________________________
57 AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
58 AliCFCutBase(name, title),
63 fMinNClustersTPCPID(0),
68 fTRDtrackletsExact(0),
72 fTPCclusterRatioDef(0),
73 fFractionTPCShared(-1.0),
74 fOppSideIPcut(kFALSE),
79 fListKinkMothers(1000),
80 fNumberKinkMothers(0),
86 // Default Constructor
88 //printf("Set the number of min ITS clusters %d\n",(Int_t)fMinNbITScls);
89 memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
90 memset(fIPcutParam, 0, sizeof(Float_t) * 4);
93 //______________________________________________________
94 AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
97 fCutCorrelation(c.fCutCorrelation),
98 fRequirements(c.fRequirements),
99 fMinNClustersTPC(c.fMinNClustersTPC),
100 fMinNClustersTPCPID(c.fMinNClustersTPCPID),
101 fClusterRatioTPC(c.fClusterRatioTPC),
102 fMinTrackletsTRD(c.fMinTrackletsTRD),
103 fMaxChi2TRD(c.fMaxChi2TRD),
104 fMinNbITScls(c.fMinNbITScls),
105 fTRDtrackletsExact(c.fTRDtrackletsExact),
106 fPixelITS(c.fPixelITS),
107 fDriftITS(c.fDriftITS),
108 fTPCclusterDef(c.fTPCclusterDef),
109 fTPCclusterRatioDef(c.fTPCclusterRatioDef),
110 fFractionTPCShared(c.fFractionTPCShared),
111 fOppSideIPcut(c.fOppSideIPcut),
112 fTOFsignalDx(c.fTOFsignalDx),
113 fTOFsignalDz(c.fTOFsignalDz),
114 fMagField(c.fMagField),
115 fAODFilterBit(c.fAODFilterBit),
116 fListKinkMothers(c.fListKinkMothers),
117 fNumberKinkMothers(c.fNumberKinkMothers),
124 // Performs a deep copy
126 memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
127 memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
130 fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
131 if(fQAlist) fQAlist->SetOwner();
135 //______________________________________________________
136 AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
138 // Assignment operator
141 AliCFCutBase::operator=(c);
143 fCutCorrelation = c.fCutCorrelation;
144 fRequirements = c.fRequirements;
145 fClusterRatioTPC = c.fClusterRatioTPC;
146 fMinNClustersTPC = c.fMinNClustersTPC;
147 fMinNClustersTPCPID = c.fMinNClustersTPCPID;
148 fMinTrackletsTRD = c.fMinTrackletsTRD;
149 fMaxChi2TRD = c.fMaxChi2TRD;
150 fMinNbITScls = c.fMinNbITScls;
151 fTRDtrackletsExact = c.fTRDtrackletsExact;
152 fPixelITS = c.fPixelITS;
153 fDriftITS = c.fDriftITS;
154 fTPCclusterDef = c.fTPCclusterDef;
155 fTPCclusterRatioDef = c.fTPCclusterRatioDef;
156 fFractionTPCShared = c.fFractionTPCShared;
157 fOppSideIPcut = c.fOppSideIPcut;
158 fTOFsignalDx = c.fTOFsignalDx;
159 fTOFsignalDz = c.fTOFsignalDz;
160 fMagField = c.fMagField;
161 fAODFilterBit = c.fAODFilterBit;
162 fListKinkMothers = c.fListKinkMothers;
163 fNumberKinkMothers = c.fNumberKinkMothers;
165 fDebugLevel = c.fDebugLevel;
166 memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
167 memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
170 fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
171 if(fQAlist) fQAlist->SetOwner();
177 //______________________________________________________
178 AliHFEextraCuts::~AliHFEextraCuts(){
182 if(fQAlist) delete fQAlist;
185 //______________________________________________________
186 void AliHFEextraCuts::SetRecEventInfo(const TObject *event){
188 // Set Virtual event an make a copy
191 AliError("Pointer to AliVEvent !");
194 TString className(event->ClassName());
195 if (! (className.CompareTo("AliESDEvent")==0 || className.CompareTo("AliAODEvent")==0)) {
196 AliError("argument must point to an AliESDEvent or AliAODEvent !");
199 fEvent = (AliVEvent*) event;
201 AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
203 // Initialize lookup for kink mother rejection
204 if(aodevent->GetNumberOfVertices() > fListKinkMothers.GetSize()) fListKinkMothers.Set(aodevent->GetNumberOfVertices());
205 fNumberKinkMothers = 0;
206 for(Int_t ivtx = 0; ivtx < aodevent->GetNumberOfVertices(); ivtx++){
207 AliAODVertex *aodvtx = aodevent->GetVertex(ivtx);
208 if(aodvtx->GetType() != AliAODVertex::kKink) continue;
209 AliAODTrack *mother = (AliAODTrack *) aodvtx->GetParent();
210 if(!mother) continue;
211 Int_t idmother = mother->GetID();
212 fListKinkMothers[fNumberKinkMothers++] = idmother;
217 //______________________________________________________
218 Bool_t AliHFEextraCuts::IsSelected(TObject *o){
220 // Steering function for the track selection
222 TClass *type = o->IsA();
223 AliDebug(2, Form("Object type %s", type->GetName()));
224 if((type == AliESDtrack::Class()) || (type == AliAODTrack::Class())){
225 return CheckRecCuts(dynamic_cast<AliVTrack *>(o));
227 return CheckMCCuts(dynamic_cast<AliVParticle *>(o));
230 //______________________________________________________
231 Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
233 // Checks cuts on reconstructed tracks
234 // returns true if track is selected
235 // QA histograms are filled before track selection and for
236 // selected tracks after track selection
238 AliDebug(1, Form("%s: Called", GetName()));
239 ULong64_t survivedCut = 0; // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
240 if(IsQAOn()) FillQAhistosRec(track, kBeforeCuts);
242 Float_t impactR = -999.;
243 Float_t impactZ = -999.;
244 Double_t hfeimpactR, hfeimpactnsigmaR;
245 Double_t hfeimpactRcut, hfeimpactnsigmaRcut;
246 Double_t maximpactRcut;
247 GetImpactParameters(track, impactR, impactZ);
248 if(TESTBIT(fRequirements, kMinHFEImpactParamR) || TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)||TESTBIT(fRequirements, kMinHFEImpactParamRcharge)){
249 // Protection for PbPb
250 GetHFEImpactParameterCuts(track, hfeimpactRcut, hfeimpactnsigmaRcut);
251 GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
253 Int_t nclsITS = GetITSNbOfcls(track);
254 UInt_t nclsTPC = GetTPCncls(track);
255 UInt_t nclsTPCPID = track->GetTPCsignalN();
256 // printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
257 Float_t fractionSharedClustersTPC = GetTPCsharedClustersRatio(track);
258 Double_t ratioTPC = GetTPCclusterRatio(track);
259 UChar_t trdTracklets;
260 trdTracklets = track->GetTRDntrackletsPID();
261 Float_t trdchi2=-999.;
262 trdchi2=GetTRDchi(track);
263 UChar_t itsPixel = track->GetITSClusterMap();
264 Int_t status1 = GetITSstatus(track, 0);
265 Int_t status2 = GetITSstatus(track, 1);
266 Bool_t statusL0 = CheckITSstatus(status1);
267 Bool_t statusL1 = CheckITSstatus(status2);
268 Double_t tofsignalDx = 0.0;
269 Double_t tofsignalDz = 0.0;
270 GetTOFsignalDxDz(track,tofsignalDx,tofsignalDz);
272 if(TESTBIT(fRequirements, kTPCfractionShared)) {
273 // cut on max fraction of shared TPC clusters
274 if(TMath::Abs(fractionSharedClustersTPC) < fFractionTPCShared) SETBIT(survivedCut, kTPCfractionShared);
276 if(TESTBIT(fRequirements, kMinImpactParamR)){
277 // cut on min. Impact Parameter in Radial direction
278 if(TMath::Abs(impactR) >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
280 if(TESTBIT(fRequirements, kMinImpactParamZ)){
281 // cut on min. Impact Parameter in Z direction
282 if(TMath::Abs(impactZ) >= fImpactParamCut[1]) SETBIT(survivedCut, kMinImpactParamZ);
284 if(TESTBIT(fRequirements, kMaxImpactParamR)){
285 // cut on max. Impact Parameter in Radial direction
286 if(TMath::Abs(impactR) <= fImpactParamCut[2]) SETBIT(survivedCut, kMaxImpactParamR);
288 if(TESTBIT(fRequirements, kMaxImpactParameterRpar)) {
289 // HFE Impact parameter cut
290 GetMaxImpactParameterCutR(track,maximpactRcut);
291 if(TMath::Abs(impactR) >= maximpactRcut) SETBIT(fRequirements, kMaxImpactParameterRpar);
293 if(TESTBIT(fRequirements, kMaxImpactParamZ)){
294 // cut on max. Impact Parameter in Z direction
295 if(TMath::Abs(impactZ) <= fImpactParamCut[3]) SETBIT(survivedCut, kMaxImpactParamZ);
297 if(TESTBIT(fRequirements, kMinHFEImpactParamR)){
298 // cut on min. HFE Impact Parameter in Radial direction
299 if(TMath::Abs(hfeimpactR) >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamR);
301 if(TESTBIT(fRequirements, kMinHFEImpactParamRcharge)){
302 // cut on min. HFE Impact Parameter in Radial direction, multiplied by particle charge
303 Double_t charge = (Double_t)track->Charge();
305 if(fMagField < 0)charge *= -1.;//the IP distribution side to be chosen depends on magnetic field polarization
306 if(fOppSideIPcut)charge*=-1.;//in case we choose to select electrons from the side of the photon peak
307 Double_t hfeimpactRcharge = hfeimpactR*charge;//switch selected side of the peak
308 if(hfeimpactRcharge >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamRcharge);
310 if(TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
311 // cut on max. HFE Impact Parameter n sigma in Radial direction
312 // if(fAbsHFEImpactParamNsigmaR) {
313 // //if((TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) && (TMath::Abs(hfeimpactnsigmaR) < 8)) { //mj debug
314 // if(TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) {
315 // SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
316 // // printf("0: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
320 if(hfeimpactnsigmaRcut>0){
321 if(hfeimpactnsigmaR >= hfeimpactnsigmaRcut) {
322 SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
323 //printf("1: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
327 if(hfeimpactnsigmaR <= hfeimpactnsigmaRcut) {
328 SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
329 //printf("2: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
334 if(TESTBIT(fRequirements, kClusterRatioTPC)){
335 // cut on min ratio of found TPC clusters vs findable TPC clusters
336 if(ratioTPC >= fClusterRatioTPC) SETBIT(survivedCut, kClusterRatioTPC);
338 if(TESTBIT(fRequirements, kMinTrackletsTRD)){
339 // cut on minimum number of TRD tracklets
340 AliDebug(1, Form("%s: Min TRD cut: [%d|%d], Require exact number [%s]\n", GetName(), fMinTrackletsTRD, trdTracklets, fTRDtrackletsExact ? "Yes" : "No"));
341 if(fTRDtrackletsExact){
342 if(trdTracklets == fMinTrackletsTRD) {
343 SETBIT(survivedCut, kMinTrackletsTRD);
344 AliDebug(1, Form("%s: Track Selected", GetName()));
347 if(trdTracklets >= fMinTrackletsTRD){
348 SETBIT(survivedCut, kMinTrackletsTRD);
349 AliDebug(1, Form("%s: Track Selected", GetName()));
351 //printf("Min number of tracklets %d\n",fMinTrackletsTRD);
355 if(TESTBIT(fRequirements, kMaxTRDChi2)){
357 AliDebug(1, Form("%s: Cut on TRD chi2: [%f|%f]\n", GetName(),fMaxChi2TRD, trdchi2));
358 if(trdchi2 < fMaxChi2TRD) {
359 SETBIT(survivedCut, kMaxTRDChi2);
360 AliDebug(1,Form("%s: survived %f\n",GetName(),trdchi2));
364 if(TESTBIT(fRequirements, kMinNbITScls)){
365 // cut on minimum number of ITS clusters
366 //printf(Form("Min ITS clusters: [%d|%d]\n", (Int_t)fMinNbITScls, nclsITS));
367 AliDebug(1, Form("%s: Min ITS clusters: [%d|%d]\n", GetName(), fMinNbITScls, nclsITS));
368 if(nclsITS >= ((Int_t)fMinNbITScls)) SETBIT(survivedCut, kMinNbITScls);
371 if(TESTBIT(fRequirements, kMinNClustersTPC)){
372 // cut on minimum number of TPC tracklets
373 //printf(Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
374 AliDebug(1, Form("%s: Min TPC cut: [%d|%d]\n", GetName(), fMinNClustersTPC, nclsTPC));
375 if(nclsTPC >= fMinNClustersTPC) SETBIT(survivedCut, kMinNClustersTPC);
377 if(TESTBIT(fRequirements, kMinNClustersTPCPID)){
378 AliDebug(1, Form("%s: Min TPC PID cut: [%d|%d]\n", GetName(), fMinNClustersTPCPID, nclsTPCPID));
379 if(nclsTPCPID >= fMinNClustersTPCPID) SETBIT(survivedCut, kMinNClustersTPCPID);
381 if(TESTBIT(fRequirements, kDriftITS)){
382 //printf("Require drift\n");
385 if(TESTBIT(itsPixel, 2)) SETBIT(survivedCut, kDriftITS);
388 SETBIT(survivedCut, kDriftITS);
392 if(TESTBIT(fRequirements, kPixelITS)){
393 // cut on ITS pixel layers
394 AliDebug(1, Form("%s: ITS cluster Map: ", GetName()));
395 //PrintBitMap(itsPixel);
398 AliDebug(2, "First");
400 if(TESTBIT(itsPixel, 0))
401 SETBIT(survivedCut, kPixelITS);
403 if(fCheck && !statusL0)
404 SETBIT(survivedCut, kPixelITS);
407 //printf("Second\n");
408 AliDebug(2, "Second");
409 if(TESTBIT(itsPixel, 1))
410 SETBIT(survivedCut, kPixelITS);
412 if(fCheck && !statusL1)
413 SETBIT(survivedCut, kPixelITS);
418 if(TESTBIT(itsPixel,0) && TESTBIT(itsPixel,1))
419 SETBIT(survivedCut, kPixelITS);
421 if(fCheck && !(statusL0 && statusL1))
422 SETBIT(survivedCut, kPixelITS);
427 if(TESTBIT(itsPixel, 0) || TESTBIT(itsPixel, 1))
428 SETBIT(survivedCut, kPixelITS);
430 if(fCheck && !(statusL0 || statusL1))
431 SETBIT(survivedCut, kPixelITS);
433 case kExclusiveSecond:
434 AliDebug(2, "Exlusive second");
435 if(fCheck){ // Cut out tracks which pass a dead ITS Layer 0
436 if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0) && statusL0)
437 SETBIT(survivedCut, kPixelITS);
439 if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0))
440 SETBIT(survivedCut, kPixelITS);
444 // No cut applied, set as survived
445 SETBIT(survivedCut, kPixelITS);
448 // default, defined as no cut applied
450 SETBIT(survivedCut, kPixelITS);
453 AliDebug(1, Form("%s: Survived Cut: %s\n", GetName(), TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"));
456 if(TESTBIT(fRequirements, kTOFPID)){
458 if(track->GetStatus() & AliESDtrack::kTOFpid) SETBIT(survivedCut, kTOFPID);
461 if(TESTBIT(fRequirements, kTOFmismatch)){
462 // cut on TOF mismatch
463 if(!(track->GetStatus() & AliESDtrack::kTOFmismatch)) SETBIT(survivedCut, kTOFmismatch);
466 if(TESTBIT(fRequirements, kTOFlabel)){
467 if(MatchTOFlabel(track)) SETBIT(survivedCut, kTOFlabel);
470 if(TESTBIT(fRequirements, kTPCPIDCleanUp)){
471 // cut on TPC PID cleanup
472 Bool_t fBitsAboveThreshold=GetTPCCountSharedMapBitsAboveThreshold(track);
473 if((fBitsAboveThreshold==0)&&(nclsTPCPID>80)) SETBIT(survivedCut, kTPCPIDCleanUp);
476 if(TESTBIT(fRequirements, kEMCALmatch)){
477 if(track->GetStatus() && AliESDtrack::kEMCALmatch) SETBIT(survivedCut, kEMCALmatch);
480 if(TESTBIT(fRequirements, kRejectKinkDaughter)){
481 //printf("test daughter\n");
482 if(!IsKinkDaughter(track)) SETBIT(survivedCut, kRejectKinkDaughter);
483 //else printf("Found kink daughter\n");
486 if(TESTBIT(fRequirements, kRejectKinkMother)){
487 //printf("test mother\n");
488 if(!IsKinkMother(track)) SETBIT(survivedCut, kRejectKinkMother);
489 //else printf("Found kink mother\n");
491 if(TESTBIT(fRequirements, kTOFsignalDxy)){
492 // cut on TOF matching cluster
493 if((TMath::Abs(tofsignalDx) <= fTOFsignalDx) && (TMath::Abs(tofsignalDz) <= fTOFsignalDz)) SETBIT(survivedCut, kTOFsignalDxy);
495 if(TESTBIT(fRequirements, kITSpattern)){
496 // cut on ITS pattern (every layer with a working ITS module must have an ITS cluster)
497 if(CheckITSpattern(track)) SETBIT(survivedCut, kITSpattern);
499 if(TESTBIT(fRequirements, kAODFilterBit)){
500 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
502 Int_t aodfilter = 1 << fAODFilterBit;
503 if(aodtrack->TestFilterBit(aodfilter)) SETBIT(survivedCut, kAODFilterBit);
507 if(fRequirements == survivedCut){
511 AliDebug(2, Form("%s: Track Survived cuts\n", GetName()));
512 if(IsQAOn()) FillQAhistosRec(track, kAfterCuts);
515 AliDebug(2, Form("%s: Track cut", GetName()));
516 if(IsQAOn()) FillCutCorrelation(survivedCut);
520 //______________________________________________________
521 Bool_t AliHFEextraCuts::CheckMCCuts(AliVParticle */*track*/) const {
523 // Checks cuts on Monte Carlo tracks
524 // returns true if track is selected
525 // QA histograms are filled before track selection and for
526 // selected tracks after track selection
528 return kTRUE; // not yet implemented
531 //______________________________________________________
532 void AliHFEextraCuts::FillQAhistosRec(AliVTrack *track, UInt_t when){
534 // Fill the QA histograms for ESD tracks
535 // Function can be called before cuts or after cut application (second argument)
537 Float_t impactR, impactZ;
538 GetImpactParameters(track, impactR, impactZ);
540 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(0 + when * fgkNQAhistos)))) htmp->Fill(impactR);
541 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(1 + when * fgkNQAhistos)))) htmp->Fill(impactZ);
542 // printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
543 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(2 + when * fgkNQAhistos)))) htmp->Fill(GetTPCclusterRatio(track));
544 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(3 + when * fgkNQAhistos)))) htmp->Fill(track->GetTRDntrackletsPID());
545 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(4 + when * fgkNQAhistos)))) htmp->Fill(GetTPCncls(track));
546 UChar_t itsPixel = track->GetITSClusterMap();
547 TH1 *pixelHist = dynamic_cast<TH1F *>(fQAlist->At(5 + when * fgkNQAhistos));
548 //Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
550 Double_t firstEntry = 0.5;
551 if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
552 pixelHist->Fill(firstEntry + 3);
554 if(itsPixel & BIT(0)){
555 pixelHist->Fill(firstEntry);
556 if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
557 else pixelHist->Fill(firstEntry + 4);
559 if(itsPixel & BIT(1)){
560 pixelHist->Fill(firstEntry + 1);
561 if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
565 // Fill histogram with the status bits
566 TH1 *hStatusBits = dynamic_cast<TH1 *>(fQAlist->At(6 + when * fgkNQAhistos));
568 hStatusBits->Fill(0); // Fill first bin with all tracks
569 if(track->GetStatus() && AliESDtrack::kTOFpid) hStatusBits->Fill(1);
570 if(!(track->GetStatus() && AliESDtrack::kTOFmismatch)) hStatusBits->Fill(2);
571 if(track->GetStatus() && AliESDtrack::kEMCALmatch) hStatusBits->Fill(3);
572 if(GetTPCCountSharedMapBitsAboveThreshold(track)==0) hStatusBits->Fill(4);
574 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(7 + when * fgkNQAhistos)))) htmp->Fill(track->GetTPCsignalN());
576 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(8 + when * fgkNQAhistos)))) htmp->Fill(GetTRDchi(track));
578 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(9 + when * fgkNQAhistos)))) htmp->Fill(GetITSNbOfcls(track));
581 // //______________________________________________________
582 // void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
584 // // Fill the QA histograms for MC tracks
585 // // Function can be called before cuts or after cut application (second argument)
586 // // Not yet implemented
590 //______________________________________________________
591 void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
593 // Fill cut correlation histograms for tracks that didn't pass cuts
595 TH2 *correlation = dynamic_cast<TH2F *>(fQAlist->At(2 * fgkNQAhistos));
597 for(Int_t icut = 0; icut < kNcuts; icut++){
598 if(!TESTBIT(fRequirements, icut)) continue;
599 for(Int_t jcut = icut; jcut < kNcuts; jcut++){
600 if(!TESTBIT(fRequirements, jcut)) continue;
601 if(TESTBIT(survivedCut, icut) && TESTBIT(survivedCut, jcut))
602 correlation->Fill(icut, jcut);
608 //______________________________________________________
609 void AliHFEextraCuts::AddQAHistograms(TList *qaList){
612 // For each cut a histogram before and after track cut is created
613 // Histos before respectively after cut are stored in different lists
614 // Additionally a histogram with the cut correlation is created and stored
615 // in the top directory
620 TString cutstr[2] = {"before", "after"};
622 if(!fQAlist) fQAlist = new TList; // for internal representation, not owner
623 for(Int_t icond = 0; icond < 2; icond++){
624 qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamR%s",GetName(),cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0 + icond * fgkNQAhistos);
625 fQAlist->AddAt(histo1D, 0 + icond * fgkNQAhistos);
626 histo1D->GetXaxis()->SetTitle("Impact Parameter");
627 histo1D->GetYaxis()->SetTitle("Number of Tracks");
628 qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamZ%s",GetName(),cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1 + icond * fgkNQAhistos);
629 fQAlist->AddAt(histo1D, 1 + icond * fgkNQAhistos);
630 histo1D->GetXaxis()->SetTitle("Impact Parameter");
631 histo1D->GetYaxis()->SetTitle("Number of Tracks");
632 qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 100, 0, 1.3)), 2 + icond * fgkNQAhistos);
633 fQAlist->AddAt(histo1D, 2 + icond * fgkNQAhistos);
634 histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
635 histo1D->GetYaxis()->SetTitle("Number of Tracks");
636 qaList->AddAt((histo1D = new TH1F(Form("%s_trdTracklets%s",GetName(),cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3 + icond * fgkNQAhistos);
637 fQAlist->AddAt(histo1D, 3 + icond * fgkNQAhistos);
638 histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
639 histo1D->GetYaxis()->SetTitle("Number of Tracks");
640 qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClusters%s",GetName(),cutstr[icond].Data()), "Number of TPC clusters", 161, 0, 160)), 4 + icond * fgkNQAhistos);
641 fQAlist->AddAt(histo1D, 4 + icond * fgkNQAhistos);
642 histo1D->GetXaxis()->SetTitle("Number of TPC clusters");
643 histo1D->GetYaxis()->SetTitle("Number of Tracks");
644 qaList->AddAt((histo1D = new TH1F(Form("%s_itsPixel%s",GetName(),cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 5 + icond * fgkNQAhistos);
645 fQAlist->AddAt(histo1D, 5 + icond * fgkNQAhistos);
646 histo1D->GetXaxis()->SetTitle("ITS Pixel");
647 histo1D->GetYaxis()->SetTitle("Number of Tracks");
648 Int_t first = histo1D->GetXaxis()->GetFirst();
649 TString binNames[6] = { "First", "Second", "Both", "None", "Exclusive First", "Exclusive Second"};
650 for(Int_t ilabel = 0; ilabel < 6; ilabel++)
651 histo1D->GetXaxis()->SetBinLabel(first + ilabel, binNames[ilabel].Data());
652 qaList->AddAt((histo1D = new TH1F(Form("%s_trackStatus%s",GetName(),cutstr[icond].Data()), "Track Status", 5, 0, 5)), 6 + icond * fgkNQAhistos);
653 fQAlist->AddAt(histo1D, 6 + icond * fgkNQAhistos);
654 histo1D->GetXaxis()->SetTitle("Track Status Bit");
655 histo1D->GetYaxis()->SetTitle("Number of tracks");
656 TString tsnames[5] = {"All", "TOFPID", "No TOFmismatch", "EMCALmatch","No TPC shared bit"};
657 for(Int_t istat = 0; istat < 5; istat++) histo1D->GetXaxis()->SetBinLabel(istat + 1, tsnames[istat].Data());
658 qaList->AddAt((histo1D = new TH1F(Form("%s_tpcdEdxClusters%s",GetName(),cutstr[icond].Data()), "Number of TPC clusters for dEdx calculation", 161, 0, 160)), 7 + icond * fgkNQAhistos);
659 fQAlist->AddAt(histo1D, 7 + icond * fgkNQAhistos);
660 histo1D->GetXaxis()->SetTitle("Number of TPC clusters for dEdx calculation");
661 histo1D->GetYaxis()->SetTitle("counts");
662 qaList->AddAt((histo1D = new TH1F(Form("%s_trdchi2perTracklet%s",GetName(),cutstr[icond].Data()), "chi2 per TRD tracklet", 100, 0, 10)), 8 + icond * fgkNQAhistos);
663 fQAlist->AddAt(histo1D, 8 + icond * fgkNQAhistos);
664 histo1D->GetXaxis()->SetTitle("Chi2 per TRD Tracklet");
665 histo1D->GetYaxis()->SetTitle("Number of Tracks");
666 qaList->AddAt((histo1D = new TH1F(Form("%s_ITScls%s",GetName(),cutstr[icond].Data()), "ITScls", 6, 0, 6)), 9 + icond * fgkNQAhistos);
667 fQAlist->AddAt(histo1D, 9 + icond * fgkNQAhistos);
668 histo1D->GetXaxis()->SetTitle("ITScls");
669 histo1D->GetYaxis()->SetTitle("Number of ITS cls");
672 // Add cut correlation
673 qaList->AddAt((histo2D = new TH2F(Form("%s_cutcorrelation",GetName()), "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2 * fgkNQAhistos);
674 fQAlist->AddAt(histo2D, 2 * fgkNQAhistos);
675 TString labels[kNcuts] = {"MinImpactParamR", "MaxImpactParamR", "MinImpactParamZ", "MaxImpactParamZ", "ClusterRatioTPC", "MinTrackletsTRD", "ITSpixel", "kMinHFEImpactParamR", "kMinHFEImpactParamNsigmaR", "TPC Number of clusters" "Fraction Shared TPC clusters", "TOFPID", "No TOFmismatch", "EMCALmatch", "ImpactParam"};
676 Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
677 for(Int_t icut = 0; icut < kNcuts; icut++){
678 histo2D->GetXaxis()->SetBinLabel(firstx + icut, labels[icut].Data());
679 histo2D->GetYaxis()->SetBinLabel(firsty + icut, labels[icut].Data());
683 //______________________________________________________
684 void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
685 for(Int_t ibit = 32; ibit--; )
686 printf("%d", bitmap & BIT(ibit) ? 1 : 0);
690 //______________________________________________________
691 Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
693 // Check whether ITS area is dead
697 case 2: status = kFALSE; break;
698 case 3: status = kFALSE; break;
699 case 7: status = kFALSE; break;
700 default: status = kTRUE;
705 //______________________________________________________
706 Int_t AliHFEextraCuts::GetITSstatus(const AliVTrack * const track, Int_t layer) const {
708 // Check ITS layer status
711 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
714 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
715 if(esdtrack) esdtrack->GetITSModuleIndexInfo(layer, det, status, xloc, zloc);
721 //______________________________________________________
722 Bool_t AliHFEextraCuts::GetTPCCountSharedMapBitsAboveThreshold(AliVTrack *track){
724 // Checks if number of shared bits is above 1
728 if(track->IsA() == AliESDtrack::Class())
729 shared = &((static_cast<AliESDtrack *>(track))->GetTPCSharedMap());
730 else if(track->IsA() == AliAODTrack::Class())
731 shared = &((static_cast<AliAODTrack *>(track))->GetTPCSharedMap());
735 if(shared->CountBits() >= 2) nsharebit=1;
741 //______________________________________________________
742 UInt_t AliHFEextraCuts::GetTPCncls(AliVTrack *track){
744 // Get Number of findable clusters in the TPC
746 Int_t nClusters = 0; // in case no Information available consider all clusters findable
747 TClass *type = track->IsA();
748 if(TESTBIT(fTPCclusterDef, kFoundIter1)){
749 if(type == AliESDtrack::Class()){
750 AliDebug(2, ("Using def kFoundIter1"));
751 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
752 nClusters = esdtrack->GetTPCNclsIter1();
754 AliDebug(2, ("Number of clusters in iteration 1 not available for AOD tracks"));
756 } else if(TESTBIT(fTPCclusterDef, kCrossedRows)){
757 AliDebug(2, ("Using def kCrossedRows"));
758 nClusters = static_cast<UInt_t>(track->GetTPCClusterInfo(2,1));
759 } else if(TESTBIT(fTPCclusterDef, kFound)){
760 AliDebug(2, ("Using def kFound"));
761 nClusters = track->GetTPCNcls();
764 AliDebug(2, ("Using def kFoundAll"));
765 if(type == AliESDtrack::Class()){
766 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
767 const TBits &clusterTPC = esdtrack->GetTPCClusterMap();
768 nClusters = clusterTPC.CountBits();
771 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
772 const TBits &clusterTPC = aodtrack->GetTPCClusterMap();
773 nClusters = clusterTPC.CountBits();
779 //______________________________________________________
780 Float_t AliHFEextraCuts::GetTPCsharedClustersRatio(AliVTrack *track){
782 // Get fraction of shared TPC clusters
784 Float_t fracClustersTPCShared = 0.0;
785 Int_t nClustersTPC = track->GetTPCNcls();
786 Int_t nClustersTPCShared = 0;
787 TClass *type = track->IsA();
788 if(type == AliESDtrack::Class()){
789 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
790 nClustersTPCShared = esdtrack->GetTPCnclsS();
791 } else if(type == AliAODTrack::Class()){
792 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
793 const TBits &shared = aodtrack->GetTPCSharedMap();
794 nClustersTPCShared = shared.CountBits(0) - shared.CountBits(159);
796 if (nClustersTPC!=0) {
797 fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
799 return fracClustersTPCShared;
802 //______________________________________________________
803 Double_t AliHFEextraCuts::GetTPCclusterRatio(AliVTrack *track){
805 // Get Ratio of found / findable clusters for different definitions
806 // Only implemented for ESD tracks
808 Double_t clusterRatio = 1.; // in case no Information available consider all clusters findable
809 TClass *type = track->IsA();
810 if(TESTBIT(fTPCclusterRatioDef, kCROverFindable)){
811 AliDebug(2, "Using ratio def kCROverFindable");
812 clusterRatio = track->GetTPCNclsF() ? track->GetTPCClusterInfo(2,1)/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // crossed rows/findable
813 } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverCR)){
814 AliDebug(2, "Using ratio def kFoundOverCR");
815 clusterRatio = track->GetTPCClusterInfo(2,0); // found/crossed rows
816 } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindableIter1)){
817 if(type == AliESDtrack::Class()){
818 AliDebug(2, "Using ratio def kFoundOverFindableIter1");
819 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
820 clusterRatio = esdtrack->GetTPCNclsFIter1() ? static_cast<Double_t>(esdtrack->GetTPCNclsIter1())/static_cast<Double_t>(esdtrack->GetTPCNclsFIter1()) : 1.; // crossed
822 AliDebug(2, "Cluster ratio after iteration 1 not possible for AOD tracks");
825 } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindable)) {
826 AliDebug(2, "Using ratio def kFoundOverFindable");
827 clusterRatio = track->GetTPCNclsF() ? static_cast<Double_t>(track->GetTPCNcls())/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // found/findable
830 AliDebug(2, "Using ratio def kFoundAllOverFindable");
831 Int_t allclusters = 0;
832 if(type == AliESDtrack::Class()){
833 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
834 const TBits &clusterTPC = esdtrack->GetTPCClusterMap();
835 allclusters = clusterTPC.CountBits();
838 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
839 const TBits &clusterTPC = aodtrack->GetTPCClusterMap();
840 allclusters = clusterTPC.CountBits();
842 clusterRatio = track->GetTPCNclsF() ? static_cast<Double_t>(allclusters)/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // foundall/findable
847 //___________________________________________________
848 void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z){
850 // Get impact parameter
853 const Double_t kBeampiperadius=3.;
854 Double_t dcaD[2]={-999.,-999.},
855 covD[3]={-999.,-999.,-999.};
856 TClass *type = track->IsA();
857 if(type == AliESDtrack::Class()){
858 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
859 esdtrack->GetImpactParameters(radial, z);
861 else if(type == AliAODTrack::Class()){
864 AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
866 AliDebug(1, "No aod event available\n");
870 //Case ESD track: take copy constructor
871 AliAODTrack *aodtrack = NULL;
872 AliAODTrack *tmptrack = dynamic_cast<AliAODTrack *>(track);
873 if(tmptrack) aodtrack = new AliAODTrack(*tmptrack);
875 AliAODVertex *vtxAODSkip = aodevent->GetPrimaryVertex();
876 if(!vtxAODSkip) return;
877 AliExternalTrackParam etp; etp.CopyFromVTrack(aodtrack);
878 if(etp.PropagateToDCA(vtxAODSkip, aodevent->GetMagneticField(), kBeampiperadius, dcaD, covD)) {
882 //if(vtxAODSkip) delete vtxAODSkip;
883 if(aodtrack) delete aodtrack;
886 //______________________________________________________
887 Bool_t AliHFEextraCuts::IsKinkDaughter(AliVTrack *track){
891 TClass *type = track->IsA();
892 if(type == AliESDtrack::Class()){
893 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
894 if(!esdtrack) return kFALSE;
895 if(esdtrack->GetKinkIndex(0)>0) return kTRUE;
899 else if(type == AliAODTrack::Class()){
900 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
902 AliAODVertex *aodvertex = aodtrack->GetProdVertex();
903 if(!aodvertex) return kFALSE;
904 if(aodvertex->GetType()==AliAODVertex::kKink) return kTRUE;
912 //______________________________________________________
913 Bool_t AliHFEextraCuts::IsKinkMother(AliVTrack *track){
918 TClass *type = track->IsA();
919 if(type == AliESDtrack::Class()){
920 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
921 if(!esdtrack) return kFALSE;
922 if(esdtrack->GetKinkIndex(0)!=0) return kTRUE;
924 } else if(type == AliAODTrack::Class()){
925 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
926 for(int ikink = 0; ikink < fNumberKinkMothers; ikink++){
927 if(fListKinkMothers[ikink] == aodtrack->GetID()) return kTRUE;
935 //______________________________________________________
936 Float_t AliHFEextraCuts::GetTRDchi(AliVTrack *track){
941 TClass *type = track->IsA();
942 if(type == AliESDtrack::Class()){
943 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
944 ntls = esdtrack->GetTRDntracklets();
945 return ntls ? esdtrack->GetTRDchi2()/ntls : -999;
947 else if(type == AliAODTrack::Class()){
948 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
958 //______________________________________________________
959 Int_t AliHFEextraCuts::GetITSNbOfcls(AliVTrack *track){
961 // Get ITS nb of clusters
963 TClass *type = track->IsA();
964 if(type == AliESDtrack::Class()){
965 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
966 return esdtrack->GetITSclusters(0);
969 else if(type == AliAODTrack::Class()){
970 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
972 return aodtrack->GetITSNcls();
978 //______________________________________________________
979 void AliHFEextraCuts::GetHFEImpactParameters(const AliVTrack * const track, Double_t &dcaxy, Double_t &dcansigmaxy){
981 // Get HFE impact parameter (with recalculated primary vertex)
986 AliDebug(1, "No Input event available\n");
989 TString type = track->IsA()->GetName();
990 const Double_t kBeampiperadius=3.;
991 Double_t dcaD[2]={-999.,-999.},
992 covD[3]={-999.,-999.,-999.};
993 Bool_t isRecalcVertex(kFALSE);
995 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
997 AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
999 AliDebug(1, "No esd event available\n");
1003 const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
1004 if(!vtxESDSkip) return;
1006 //case ESD track: take copy constructor
1007 const AliESDtrack *tmptrack = dynamic_cast<const AliESDtrack *>(track);
1010 if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
1011 vtxESDSkip = RemoveDaughtersFromPrimaryVtx(esdevent, track);
1012 isRecalcVertex = kTRUE;
1016 AliESDtrack esdtrack(*tmptrack);
1017 fMagField = fEvent->GetMagneticField();
1018 if(esdtrack.PropagateToDCA(vtxESDSkip, fMagField, kBeampiperadius, dcaD, covD)){
1020 if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
1022 if(isRecalcVertex) delete vtxESDSkip;
1027 //case of AOD tracks
1028 AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
1030 AliDebug(1, "No aod event available\n");
1034 AliAODVertex *vtxAODSkip = aodevent->GetPrimaryVertex();
1035 if(!vtxAODSkip) return;
1037 //Case ESD track: take copy constructor
1038 const AliAODTrack *tmptrack = dynamic_cast<const AliAODTrack *>(track);
1041 if(vtxAODSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
1042 vtxAODSkip = RemoveDaughtersFromPrimaryVtx(aodevent, track);
1043 isRecalcVertex = kTRUE;
1046 AliAODTrack aodtrack(*tmptrack);
1047 AliExternalTrackParam etp; etp.CopyFromVTrack(&aodtrack);
1048 fMagField = aodevent->GetMagneticField();
1049 if(etp.PropagateToDCA(vtxAODSkip,fMagField, kBeampiperadius, dcaD, covD)) {
1051 if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
1053 if(isRecalcVertex) delete vtxAODSkip;
1059 //______________________________________________________
1060 void AliHFEextraCuts::GetHFEImpactParameters(const AliVTrack * const track, Double_t dcaD[2], Double_t covD[3]){
1062 // Get HFE impact parameter (with recalculated primary vertex)
1065 AliDebug(1, "No Input event available\n");
1068 const Double_t kBeampiperadius=3.;
1069 TString type = track->IsA()->GetName();
1070 Bool_t isRecalcVertex(kFALSE);
1072 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
1073 //case of ESD tracks
1074 AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
1076 AliDebug(1, "No esd event available\n");
1080 // Check whether primary vertex is available
1081 const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
1082 if(!vtxESDSkip) return;
1084 const AliESDtrack *tmptrack = dynamic_cast<const AliESDtrack *>(track);
1086 //case ESD track: take copy constructor
1088 if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
1089 vtxESDSkip = RemoveDaughtersFromPrimaryVtx(esdevent, track);
1090 isRecalcVertex = kTRUE;
1093 AliESDtrack esdtrack(*tmptrack);
1094 fMagField = fEvent->GetMagneticField();
1095 esdtrack.PropagateToDCA(vtxESDSkip, fMagField, kBeampiperadius, dcaD, covD);
1097 if(isRecalcVertex) delete vtxESDSkip;
1102 //case of AOD tracks
1103 AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
1105 AliDebug(1, "No aod event available\n");
1109 // Check whether primary vertex is available
1110 AliAODVertex *vtxAODSkip = aodevent->GetPrimaryVertex();
1111 if(!vtxAODSkip) return;
1113 //Case ESD track: take copy constructor
1114 const AliAODTrack *tmptrack = dynamic_cast<const AliAODTrack *>(track);
1117 if(vtxAODSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
1118 vtxAODSkip = RemoveDaughtersFromPrimaryVtx(aodevent, track);
1119 isRecalcVertex = kTRUE;
1122 AliAODTrack aodtrack(*tmptrack);
1123 AliExternalTrackParam etp; etp.CopyFromVTrack(&aodtrack);
1124 fMagField = aodevent->GetMagneticField();
1125 etp.PropagateToDCA(vtxAODSkip, fMagField, kBeampiperadius, dcaD, covD);
1127 if(isRecalcVertex) delete vtxAODSkip;
1133 //______________________________________________________
1134 void AliHFEextraCuts::GetHFEImpactParameterCuts(const AliVTrack * const track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
1136 // Get HFE impact parameter cut(pt dependent)
1139 Double_t pt = track->Pt();
1140 hfeimpactRcut = fIPcutParam[0]+fIPcutParam[1]*exp(fIPcutParam[2]*pt); // abs R cut
1141 hfeimpactnsigmaRcut = fIPcutParam[3]; // sigma cut
1143 //______________________________________________________
1144 void AliHFEextraCuts::GetMaxImpactParameterCutR(const AliVTrack * const track, Double_t &maximpactRcut){
1146 // Get max impact parameter cut r (pt dependent)
1149 Double_t pt = track->Pt();
1151 maximpactRcut = 0.0182 + 0.035/TMath::Power(pt,1.01); // abs R cut
1153 else maximpactRcut = 9999999999.0;
1155 //______________________________________________________
1156 void AliHFEextraCuts::GetTOFsignalDxDz(const AliVTrack * const track, Double_t &tofsignalDx, Double_t &tofsignalDz){
1161 TString type = track->IsA()->GetName();
1162 if(!type.CompareTo("AliESDtrack")){
1163 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
1164 if(!esdtrack) return;
1165 tofsignalDx = esdtrack->GetTOFsignalDx();
1166 tofsignalDz = esdtrack->GetTOFsignalDz();
1171 //______________________________________________________
1172 Bool_t AliHFEextraCuts::MatchTOFlabel(const AliVTrack *const track) const {
1174 // Check whether the TOF label is the same as the track label
1176 const AliESDtrack *esdtrk(NULL);
1177 const AliAODTrack *aodtrk(NULL);
1178 int trklabel(99999), toflabel[3] = {99999,99999,99999};
1179 if((esdtrk = dynamic_cast<const AliESDtrack *>(track))){
1180 trklabel = esdtrk->GetLabel();
1181 esdtrk->GetTOFLabel(toflabel);
1182 } else if((aodtrk = dynamic_cast<const AliAODTrack *>(track))){
1183 trklabel = aodtrk->GetLabel();
1184 aodtrk->GetTOFLabel(toflabel);
1185 } else return kFALSE;
1186 if(TMath::Abs(trklabel) == TMath::Abs(toflabel[0])) return kTRUE;
1189 //______________________________________________________
1190 Bool_t AliHFEextraCuts::CheckITSpattern(const AliVTrack *const track) const {
1192 // Check if every ITS layer, which has a module which is alive, also
1193 // has an ITS cluster
1195 Bool_t patternOK(kTRUE);
1197 for(Int_t ily = 0; ily < 6; ily++){
1198 status = GetITSstatus(track, ily);
1199 if(CheckITSstatus(status)){
1200 // pixel alive, check whether layer has a cluster
1201 if(!TESTBIT(track->GetITSClusterMap(),ily)){
1202 // No cluster even though pixel is alive - reject track
1211 //---------------------------------------------------------------------------
1212 const AliVVertex* AliHFEextraCuts::RemoveDaughtersFromPrimaryVtx(const AliESDEvent * const esdevent, const AliVTrack * const track) {
1214 // This method returns a primary vertex without the daughter tracks of the
1215 // candidate and it recalculates the impact parameters and errors for ESD tracks.
1217 // The output vertex is created with "new".
1220 //const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
1222 AliVertexerTracks vertexer(fEvent->GetMagneticField());
1223 vertexer.SetITSMode();
1224 vertexer.SetMinClusters(4);
1226 skipped[0] = track->GetID();
1227 vertexer.SetSkipTracks(1,skipped);
1229 //diamond constraint
1230 vertexer.SetConstraintOn();
1231 Float_t diamondcovxy[3];
1232 esdevent->GetDiamondCovXY(diamondcovxy);
1233 Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
1234 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1235 AliESDVertex diamond(pos,cov,1.,1);
1236 vertexer.SetVtxStart(&diamond);
1238 const AliVVertex *vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
1243 //---------------------------------------------------------------------------
1244 AliAODVertex* AliHFEextraCuts::RemoveDaughtersFromPrimaryVtx(const AliAODEvent * const aod, const AliVTrack * const track) {
1246 // This method returns a primary vertex without the daughter tracks of the
1247 // candidate and it recalculates the impact parameters and errors for AOD tracks.
1248 // The output vertex is created with "new".
1251 AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
1252 if(!vtxAOD) return 0;
1253 TString title=vtxAOD->GetTitle();
1254 if(!title.Contains("VertexerTracks")) return 0;
1256 AliVertexerTracks vertexer(aod->GetMagneticField());
1258 vertexer.SetITSMode();
1259 vertexer.SetMinClusters(3);
1260 vertexer.SetConstraintOff();
1262 if(title.Contains("WithConstraint")) {
1263 Float_t diamondcovxy[3];
1264 aod->GetDiamondCovXY(diamondcovxy);
1265 Double_t pos[3]={aod->GetDiamondX(),aod->GetDiamondY(),0.};
1266 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1267 AliESDVertex diamond(pos,cov,1.,1);
1268 vertexer.SetVtxStart(&diamond);
1270 Int_t skipped[2]; for(Int_t i=0;i<2;i++) skipped[i]=-1;
1271 Int_t id = (Int_t)track->GetID();
1272 if(!(id<0)) skipped[0] = id;
1274 /*// exclude tracks with global constrained parameters
1275 Int_t nTracks=aod->GetNumberOfTracks();
1276 for(Int_t i=0; i<nTracks; i++){
1277 t = aod->GetTrack(i);
1278 if(t->TestFilterMask(512)){
1279 id = (Int_t)t->GetID();
1280 skipped[nTrksToSkip++] = id;
1284 vertexer.SetSkipTracks(1,skipped);
1285 AliESDVertex *vtxESDNew = vertexer.FindPrimaryVertex(aod);
1287 if(!vtxESDNew) return 0;
1288 if(vtxESDNew->GetNContributors()<=0) {
1289 delete vtxESDNew; vtxESDNew=NULL;
1293 // convert to AliAODVertex
1294 Double_t pos[3],cov[6],chi2perNDF;
1295 vtxESDNew->GetXYZ(pos); // position
1296 vtxESDNew->GetCovMatrix(cov); //covariance matrix
1297 chi2perNDF = vtxESDNew->GetChi2toNDF();
1298 delete vtxESDNew; vtxESDNew=NULL;
1300 AliAODVertex *vtxAODNew = new AliAODVertex(pos,cov,chi2perNDF);