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),
84 // Default Constructor
86 //printf("Set the number of min ITS clusters %d\n",(Int_t)fMinNbITScls);
87 memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
88 memset(fIPcutParam, 0, sizeof(Float_t) * 4);
91 //______________________________________________________
92 AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
95 fCutCorrelation(c.fCutCorrelation),
96 fRequirements(c.fRequirements),
97 fMinNClustersTPC(c.fMinNClustersTPC),
98 fMinNClustersTPCPID(c.fMinNClustersTPCPID),
99 fClusterRatioTPC(c.fClusterRatioTPC),
100 fMinTrackletsTRD(c.fMinTrackletsTRD),
101 fMaxChi2TRD(c.fMaxChi2TRD),
102 fMinNbITScls(c.fMinNbITScls),
103 fTRDtrackletsExact(c.fTRDtrackletsExact),
104 fPixelITS(c.fPixelITS),
105 fDriftITS(c.fDriftITS),
106 fTPCclusterDef(c.fTPCclusterDef),
107 fTPCclusterRatioDef(c.fTPCclusterRatioDef),
108 fFractionTPCShared(c.fFractionTPCShared),
109 fOppSideIPcut(c.fOppSideIPcut),
110 fTOFsignalDx(c.fTOFsignalDx),
111 fTOFsignalDz(c.fTOFsignalDz),
112 fMagField(c.fMagField),
113 fAODFilterBit(c.fAODFilterBit),
120 // Performs a deep copy
122 memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
123 memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
126 fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
127 if(fQAlist) fQAlist->SetOwner();
131 //______________________________________________________
132 AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
134 // Assignment operator
137 AliCFCutBase::operator=(c);
139 fCutCorrelation = c.fCutCorrelation;
140 fRequirements = c.fRequirements;
141 fClusterRatioTPC = c.fClusterRatioTPC;
142 fMinNClustersTPC = c.fMinNClustersTPC;
143 fMinNClustersTPCPID = c.fMinNClustersTPCPID;
144 fMinTrackletsTRD = c.fMinTrackletsTRD;
145 fMaxChi2TRD = c.fMaxChi2TRD;
146 fMinNbITScls = c.fMinNbITScls;
147 fTRDtrackletsExact = c.fTRDtrackletsExact;
148 fPixelITS = c.fPixelITS;
149 fDriftITS = c.fDriftITS;
150 fTPCclusterDef = c.fTPCclusterDef;
151 fTPCclusterRatioDef = c.fTPCclusterRatioDef;
152 fFractionTPCShared = c.fFractionTPCShared;
153 fOppSideIPcut = c.fOppSideIPcut;
154 fTOFsignalDx = c.fTOFsignalDx;
155 fTOFsignalDz = c.fTOFsignalDz;
156 fMagField = c.fMagField;
157 fAODFilterBit = c.fAODFilterBit;
159 fDebugLevel = c.fDebugLevel;
160 memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
161 memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
164 fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
165 if(fQAlist) fQAlist->SetOwner();
171 //______________________________________________________
172 AliHFEextraCuts::~AliHFEextraCuts(){
176 if(fQAlist) delete fQAlist;
179 //______________________________________________________
180 void AliHFEextraCuts::SetRecEventInfo(const TObject *event){
182 // Set Virtual event an make a copy
185 AliError("Pointer to AliVEvent !");
188 TString className(event->ClassName());
189 if (! (className.CompareTo("AliESDEvent")==0 || className.CompareTo("AliAODEvent")==0)) {
190 AliError("argument must point to an AliESDEvent or AliAODEvent !");
193 fEvent = (AliVEvent*) event;
197 //______________________________________________________
198 Bool_t AliHFEextraCuts::IsSelected(TObject *o){
200 // Steering function for the track selection
202 TClass *type = o->IsA();
203 AliDebug(2, Form("Object type %s", type->GetName()));
204 if((type == AliESDtrack::Class()) || (type == AliAODTrack::Class())){
205 return CheckRecCuts(dynamic_cast<AliVTrack *>(o));
207 return CheckMCCuts(dynamic_cast<AliVParticle *>(o));
210 //______________________________________________________
211 Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
213 // Checks cuts on reconstructed tracks
214 // returns true if track is selected
215 // QA histograms are filled before track selection and for
216 // selected tracks after track selection
218 AliDebug(1, Form("%s: Called", GetName()));
219 ULong64_t survivedCut = 0; // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
220 if(IsQAOn()) FillQAhistosRec(track, kBeforeCuts);
222 Float_t impactR = -999.;
223 Float_t impactZ = -999.;
224 Double_t hfeimpactR, hfeimpactnsigmaR;
225 Double_t hfeimpactRcut, hfeimpactnsigmaRcut;
226 Double_t maximpactRcut;
227 GetImpactParameters(track, impactR, impactZ);
228 if(TESTBIT(fRequirements, kMinHFEImpactParamR) || TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)||TESTBIT(fRequirements, kMinHFEImpactParamRcharge)){
229 // Protection for PbPb
230 GetHFEImpactParameterCuts(track, hfeimpactRcut, hfeimpactnsigmaRcut);
231 GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
233 Int_t nclsITS = GetITSNbOfcls(track);
234 UInt_t nclsTPC = GetTPCncls(track);
235 UInt_t nclsTPCPID = track->GetTPCsignalN();
236 // printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
237 Float_t fractionSharedClustersTPC = GetTPCsharedClustersRatio(track);
238 Double_t ratioTPC = GetTPCclusterRatio(track);
239 UChar_t trdTracklets;
240 trdTracklets = track->GetTRDntrackletsPID();
241 Float_t trdchi2=-999.;
242 trdchi2=GetTRDchi(track);
243 UChar_t itsPixel = track->GetITSClusterMap();
244 Int_t status1 = GetITSstatus(track, 0);
245 Int_t status2 = GetITSstatus(track, 1);
246 Bool_t statusL0 = CheckITSstatus(status1);
247 Bool_t statusL1 = CheckITSstatus(status2);
248 Double_t tofsignalDx = 0.0;
249 Double_t tofsignalDz = 0.0;
250 GetTOFsignalDxDz(track,tofsignalDx,tofsignalDz);
252 if(TESTBIT(fRequirements, kTPCfractionShared)) {
253 // cut on max fraction of shared TPC clusters
254 if(TMath::Abs(fractionSharedClustersTPC) < fFractionTPCShared) SETBIT(survivedCut, kTPCfractionShared);
256 if(TESTBIT(fRequirements, kMinImpactParamR)){
257 // cut on min. Impact Parameter in Radial direction
258 if(TMath::Abs(impactR) >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
260 if(TESTBIT(fRequirements, kMinImpactParamZ)){
261 // cut on min. Impact Parameter in Z direction
262 if(TMath::Abs(impactZ) >= fImpactParamCut[1]) SETBIT(survivedCut, kMinImpactParamZ);
264 if(TESTBIT(fRequirements, kMaxImpactParamR)){
265 // cut on max. Impact Parameter in Radial direction
266 if(TMath::Abs(impactR) <= fImpactParamCut[2]) SETBIT(survivedCut, kMaxImpactParamR);
268 if(TESTBIT(fRequirements, kMaxImpactParameterRpar)) {
269 // HFE Impact parameter cut
270 GetMaxImpactParameterCutR(track,maximpactRcut);
271 if(TMath::Abs(impactR) >= maximpactRcut) SETBIT(fRequirements, kMaxImpactParameterRpar);
273 if(TESTBIT(fRequirements, kMaxImpactParamZ)){
274 // cut on max. Impact Parameter in Z direction
275 if(TMath::Abs(impactZ) <= fImpactParamCut[3]) SETBIT(survivedCut, kMaxImpactParamZ);
277 if(TESTBIT(fRequirements, kMinHFEImpactParamR)){
278 // cut on min. HFE Impact Parameter in Radial direction
279 if(TMath::Abs(hfeimpactR) >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamR);
281 if(TESTBIT(fRequirements, kMinHFEImpactParamRcharge)){
282 // cut on min. HFE Impact Parameter in Radial direction, multiplied by particle charge
283 Double_t charge = (Double_t)track->Charge();
285 if(fMagField < 0)charge *= -1.;//the IP distribution side to be chosen depends on magnetic field polarization
286 if(fOppSideIPcut)charge*=-1.;//in case we choose to select electrons from the side of the photon peak
287 Double_t hfeimpactRcharge = hfeimpactR*charge;//switch selected side of the peak
288 if(hfeimpactRcharge >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamRcharge);
290 if(TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
291 // cut on max. HFE Impact Parameter n sigma in Radial direction
292 // if(fAbsHFEImpactParamNsigmaR) {
293 // //if((TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) && (TMath::Abs(hfeimpactnsigmaR) < 8)) { //mj debug
294 // if(TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) {
295 // SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
296 // // printf("0: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
300 if(hfeimpactnsigmaRcut>0){
301 if(hfeimpactnsigmaR >= hfeimpactnsigmaRcut) {
302 SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
303 //printf("1: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
307 if(hfeimpactnsigmaR <= hfeimpactnsigmaRcut) {
308 SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
309 //printf("2: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
314 if(TESTBIT(fRequirements, kClusterRatioTPC)){
315 // cut on min ratio of found TPC clusters vs findable TPC clusters
316 if(ratioTPC >= fClusterRatioTPC) SETBIT(survivedCut, kClusterRatioTPC);
318 if(TESTBIT(fRequirements, kMinTrackletsTRD)){
319 // cut on minimum number of TRD tracklets
320 AliDebug(1, Form("%s: Min TRD cut: [%d|%d], Require exact number [%s]\n", GetName(), fMinTrackletsTRD, trdTracklets, fTRDtrackletsExact ? "Yes" : "No"));
321 if(fTRDtrackletsExact){
322 if(trdTracklets == fMinTrackletsTRD) {
323 SETBIT(survivedCut, kMinTrackletsTRD);
324 AliDebug(1, Form("%s: Track Selected", GetName()));
327 if(trdTracklets >= fMinTrackletsTRD){
328 SETBIT(survivedCut, kMinTrackletsTRD);
329 AliDebug(1, Form("%s: Track Selected", GetName()));
331 //printf("Min number of tracklets %d\n",fMinTrackletsTRD);
335 if(TESTBIT(fRequirements, kMaxTRDChi2)){
337 AliDebug(1, Form("%s: Cut on TRD chi2: [%f|%f]\n", GetName(),fMaxChi2TRD, trdchi2));
338 if(trdchi2 < fMaxChi2TRD) {
339 SETBIT(survivedCut, kMaxTRDChi2);
340 AliDebug(1,Form("%s: survived %f\n",GetName(),trdchi2));
344 if(TESTBIT(fRequirements, kMinNbITScls)){
345 // cut on minimum number of ITS clusters
346 //printf(Form("Min ITS clusters: [%d|%d]\n", (Int_t)fMinNbITScls, nclsITS));
347 AliDebug(1, Form("%s: Min ITS clusters: [%d|%d]\n", GetName(), fMinNbITScls, nclsITS));
348 if(nclsITS >= ((Int_t)fMinNbITScls)) SETBIT(survivedCut, kMinNbITScls);
351 if(TESTBIT(fRequirements, kMinNClustersTPC)){
352 // cut on minimum number of TPC tracklets
353 //printf(Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
354 AliDebug(1, Form("%s: Min TPC cut: [%d|%d]\n", GetName(), fMinNClustersTPC, nclsTPC));
355 if(nclsTPC >= fMinNClustersTPC) SETBIT(survivedCut, kMinNClustersTPC);
357 if(TESTBIT(fRequirements, kMinNClustersTPCPID)){
358 AliDebug(1, Form("%s: Min TPC PID cut: [%d|%d]\n", GetName(), fMinNClustersTPCPID, nclsTPCPID));
359 if(nclsTPCPID >= fMinNClustersTPCPID) SETBIT(survivedCut, kMinNClustersTPCPID);
361 if(TESTBIT(fRequirements, kDriftITS)){
362 //printf("Require drift\n");
365 if(TESTBIT(itsPixel, 2)) SETBIT(survivedCut, kDriftITS);
368 SETBIT(survivedCut, kDriftITS);
372 if(TESTBIT(fRequirements, kPixelITS)){
373 // cut on ITS pixel layers
374 AliDebug(1, Form("%s: ITS cluster Map: ", GetName()));
375 //PrintBitMap(itsPixel);
378 AliDebug(2, "First");
380 if(TESTBIT(itsPixel, 0))
381 SETBIT(survivedCut, kPixelITS);
383 if(fCheck && !statusL0)
384 SETBIT(survivedCut, kPixelITS);
387 //printf("Second\n");
388 AliDebug(2, "Second");
389 if(TESTBIT(itsPixel, 1))
390 SETBIT(survivedCut, kPixelITS);
392 if(fCheck && !statusL1)
393 SETBIT(survivedCut, kPixelITS);
398 if(TESTBIT(itsPixel,0) && TESTBIT(itsPixel,1))
399 SETBIT(survivedCut, kPixelITS);
401 if(fCheck && !(statusL0 && statusL1))
402 SETBIT(survivedCut, kPixelITS);
407 if(TESTBIT(itsPixel, 0) || TESTBIT(itsPixel, 1))
408 SETBIT(survivedCut, kPixelITS);
410 if(fCheck && !(statusL0 || statusL1))
411 SETBIT(survivedCut, kPixelITS);
413 case kExclusiveSecond:
414 AliDebug(2, "Exlusive second");
415 if(fCheck){ // Cut out tracks which pass a dead ITS Layer 0
416 if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0) && statusL0)
417 SETBIT(survivedCut, kPixelITS);
419 if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0))
420 SETBIT(survivedCut, kPixelITS);
424 // No cut applied, set as survived
425 SETBIT(survivedCut, kPixelITS);
428 // default, defined as no cut applied
430 SETBIT(survivedCut, kPixelITS);
433 AliDebug(1, Form("%s: Survived Cut: %s\n", GetName(), TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"));
436 if(TESTBIT(fRequirements, kTOFPID)){
438 if(track->GetStatus() & AliESDtrack::kTOFpid) SETBIT(survivedCut, kTOFPID);
441 if(TESTBIT(fRequirements, kTOFmismatch)){
442 // cut on TOF mismatch
443 if(!(track->GetStatus() & AliESDtrack::kTOFmismatch)) SETBIT(survivedCut, kTOFmismatch);
446 if(TESTBIT(fRequirements, kTPCPIDCleanUp)){
447 // cut on TPC PID cleanup
448 Bool_t fBitsAboveThreshold=GetTPCCountSharedMapBitsAboveThreshold(track);
449 if((fBitsAboveThreshold==0)&&(nclsTPCPID>80)) SETBIT(survivedCut, kTPCPIDCleanUp);
452 if(TESTBIT(fRequirements, kEMCALmatch)){
453 if(track->GetStatus() && AliESDtrack::kEMCALmatch) SETBIT(survivedCut, kEMCALmatch);
456 if(TESTBIT(fRequirements, kRejectKinkDaughter)){
457 //printf("test daughter\n");
458 if(!IsKinkDaughter(track)) SETBIT(survivedCut, kRejectKinkDaughter);
459 //else printf("Found kink daughter\n");
462 if(TESTBIT(fRequirements, kRejectKinkMother)){
463 //printf("test mother\n");
464 if(!IsKinkMother(track)) SETBIT(survivedCut, kRejectKinkMother);
465 //else printf("Found kink mother\n");
467 if(TESTBIT(fRequirements, kTOFsignalDxy)){
468 // cut on TOF matching cluster
469 if((TMath::Abs(tofsignalDx) <= fTOFsignalDx) && (TMath::Abs(tofsignalDz) <= fTOFsignalDz)) SETBIT(survivedCut, kTOFsignalDxy);
471 if(TESTBIT(fRequirements, kITSpattern)){
472 // cut on ITS pattern (every layer with a working ITS module must have an ITS cluster)
473 if(CheckITSpattern(track)) SETBIT(survivedCut, kITSpattern);
475 if(TESTBIT(fRequirements, kAODFilterBit)){
476 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
478 Int_t aodfilter = 1 << fAODFilterBit;
479 if(aodtrack->TestFilterBit(aodfilter)) SETBIT(survivedCut, kAODFilterBit);
483 if(fRequirements == survivedCut){
487 AliDebug(2, Form("%s: Track Survived cuts\n", GetName()));
488 if(IsQAOn()) FillQAhistosRec(track, kAfterCuts);
491 AliDebug(2, Form("%s: Track cut", GetName()));
492 if(IsQAOn()) FillCutCorrelation(survivedCut);
496 //______________________________________________________
497 Bool_t AliHFEextraCuts::CheckMCCuts(AliVParticle */*track*/) const {
499 // Checks cuts on Monte Carlo tracks
500 // returns true if track is selected
501 // QA histograms are filled before track selection and for
502 // selected tracks after track selection
504 return kTRUE; // not yet implemented
507 //______________________________________________________
508 void AliHFEextraCuts::FillQAhistosRec(AliVTrack *track, UInt_t when){
510 // Fill the QA histograms for ESD tracks
511 // Function can be called before cuts or after cut application (second argument)
513 Float_t impactR, impactZ;
514 GetImpactParameters(track, impactR, impactZ);
516 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(0 + when * fgkNQAhistos)))) htmp->Fill(impactR);
517 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(1 + when * fgkNQAhistos)))) htmp->Fill(impactZ);
518 // printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
519 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(2 + when * fgkNQAhistos)))) htmp->Fill(GetTPCclusterRatio(track));
520 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(3 + when * fgkNQAhistos)))) htmp->Fill(track->GetTRDntrackletsPID());
521 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(4 + when * fgkNQAhistos)))) htmp->Fill(GetTPCncls(track));
522 UChar_t itsPixel = track->GetITSClusterMap();
523 TH1 *pixelHist = dynamic_cast<TH1F *>(fQAlist->At(5 + when * fgkNQAhistos));
524 //Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
526 Double_t firstEntry = 0.5;
527 if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
528 pixelHist->Fill(firstEntry + 3);
530 if(itsPixel & BIT(0)){
531 pixelHist->Fill(firstEntry);
532 if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
533 else pixelHist->Fill(firstEntry + 4);
535 if(itsPixel & BIT(1)){
536 pixelHist->Fill(firstEntry + 1);
537 if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
541 // Fill histogram with the status bits
542 TH1 *hStatusBits = dynamic_cast<TH1 *>(fQAlist->At(6 + when * fgkNQAhistos));
544 hStatusBits->Fill(0); // Fill first bin with all tracks
545 if(track->GetStatus() && AliESDtrack::kTOFpid) hStatusBits->Fill(1);
546 if(!(track->GetStatus() && AliESDtrack::kTOFmismatch)) hStatusBits->Fill(2);
547 if(track->GetStatus() && AliESDtrack::kEMCALmatch) hStatusBits->Fill(3);
548 if(GetTPCCountSharedMapBitsAboveThreshold(track)==0) hStatusBits->Fill(4);
550 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(7 + when * fgkNQAhistos)))) htmp->Fill(track->GetTPCsignalN());
552 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(8 + when * fgkNQAhistos)))) htmp->Fill(GetTRDchi(track));
554 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(9 + when * fgkNQAhistos)))) htmp->Fill(GetITSNbOfcls(track));
557 // //______________________________________________________
558 // void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
560 // // Fill the QA histograms for MC tracks
561 // // Function can be called before cuts or after cut application (second argument)
562 // // Not yet implemented
566 //______________________________________________________
567 void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
569 // Fill cut correlation histograms for tracks that didn't pass cuts
571 TH2 *correlation = dynamic_cast<TH2F *>(fQAlist->At(2 * fgkNQAhistos));
573 for(Int_t icut = 0; icut < kNcuts; icut++){
574 if(!TESTBIT(fRequirements, icut)) continue;
575 for(Int_t jcut = icut; jcut < kNcuts; jcut++){
576 if(!TESTBIT(fRequirements, jcut)) continue;
577 if(TESTBIT(survivedCut, icut) && TESTBIT(survivedCut, jcut))
578 correlation->Fill(icut, jcut);
584 //______________________________________________________
585 void AliHFEextraCuts::AddQAHistograms(TList *qaList){
588 // For each cut a histogram before and after track cut is created
589 // Histos before respectively after cut are stored in different lists
590 // Additionally a histogram with the cut correlation is created and stored
591 // in the top directory
596 TString cutstr[2] = {"before", "after"};
598 if(!fQAlist) fQAlist = new TList; // for internal representation, not owner
599 for(Int_t icond = 0; icond < 2; icond++){
600 qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamR%s",GetName(),cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0 + icond * fgkNQAhistos);
601 fQAlist->AddAt(histo1D, 0 + icond * fgkNQAhistos);
602 histo1D->GetXaxis()->SetTitle("Impact Parameter");
603 histo1D->GetYaxis()->SetTitle("Number of Tracks");
604 qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamZ%s",GetName(),cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1 + icond * fgkNQAhistos);
605 fQAlist->AddAt(histo1D, 1 + icond * fgkNQAhistos);
606 histo1D->GetXaxis()->SetTitle("Impact Parameter");
607 histo1D->GetYaxis()->SetTitle("Number of Tracks");
608 qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 100, 0, 1.3)), 2 + icond * fgkNQAhistos);
609 fQAlist->AddAt(histo1D, 2 + icond * fgkNQAhistos);
610 histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
611 histo1D->GetYaxis()->SetTitle("Number of Tracks");
612 qaList->AddAt((histo1D = new TH1F(Form("%s_trdTracklets%s",GetName(),cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3 + icond * fgkNQAhistos);
613 fQAlist->AddAt(histo1D, 3 + icond * fgkNQAhistos);
614 histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
615 histo1D->GetYaxis()->SetTitle("Number of Tracks");
616 qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClusters%s",GetName(),cutstr[icond].Data()), "Number of TPC clusters", 161, 0, 160)), 4 + icond * fgkNQAhistos);
617 fQAlist->AddAt(histo1D, 4 + icond * fgkNQAhistos);
618 histo1D->GetXaxis()->SetTitle("Number of TPC clusters");
619 histo1D->GetYaxis()->SetTitle("Number of Tracks");
620 qaList->AddAt((histo1D = new TH1F(Form("%s_itsPixel%s",GetName(),cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 5 + icond * fgkNQAhistos);
621 fQAlist->AddAt(histo1D, 5 + icond * fgkNQAhistos);
622 histo1D->GetXaxis()->SetTitle("ITS Pixel");
623 histo1D->GetYaxis()->SetTitle("Number of Tracks");
624 Int_t first = histo1D->GetXaxis()->GetFirst();
625 TString binNames[6] = { "First", "Second", "Both", "None", "Exclusive First", "Exclusive Second"};
626 for(Int_t ilabel = 0; ilabel < 6; ilabel++)
627 histo1D->GetXaxis()->SetBinLabel(first + ilabel, binNames[ilabel].Data());
628 qaList->AddAt((histo1D = new TH1F(Form("%s_trackStatus%s",GetName(),cutstr[icond].Data()), "Track Status", 5, 0, 5)), 6 + icond * fgkNQAhistos);
629 fQAlist->AddAt(histo1D, 6 + icond * fgkNQAhistos);
630 histo1D->GetXaxis()->SetTitle("Track Status Bit");
631 histo1D->GetYaxis()->SetTitle("Number of tracks");
632 TString tsnames[5] = {"All", "TOFPID", "No TOFmismatch", "EMCALmatch","No TPC shared bit"};
633 for(Int_t istat = 0; istat < 5; istat++) histo1D->GetXaxis()->SetBinLabel(istat + 1, tsnames[istat].Data());
634 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);
635 fQAlist->AddAt(histo1D, 7 + icond * fgkNQAhistos);
636 histo1D->GetXaxis()->SetTitle("Number of TPC clusters for dEdx calculation");
637 histo1D->GetYaxis()->SetTitle("counts");
638 qaList->AddAt((histo1D = new TH1F(Form("%s_trdchi2perTracklet%s",GetName(),cutstr[icond].Data()), "chi2 per TRD tracklet", 100, 0, 10)), 8 + icond * fgkNQAhistos);
639 fQAlist->AddAt(histo1D, 8 + icond * fgkNQAhistos);
640 histo1D->GetXaxis()->SetTitle("Chi2 per TRD Tracklet");
641 histo1D->GetYaxis()->SetTitle("Number of Tracks");
642 qaList->AddAt((histo1D = new TH1F(Form("%s_ITScls%s",GetName(),cutstr[icond].Data()), "ITScls", 6, 0, 6)), 9 + icond * fgkNQAhistos);
643 fQAlist->AddAt(histo1D, 9 + icond * fgkNQAhistos);
644 histo1D->GetXaxis()->SetTitle("ITScls");
645 histo1D->GetYaxis()->SetTitle("Number of ITS cls");
648 // Add cut correlation
649 qaList->AddAt((histo2D = new TH2F(Form("%s_cutcorrelation",GetName()), "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2 * fgkNQAhistos);
650 fQAlist->AddAt(histo2D, 2 * fgkNQAhistos);
651 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"};
652 Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
653 for(Int_t icut = 0; icut < kNcuts; icut++){
654 histo2D->GetXaxis()->SetBinLabel(firstx + icut, labels[icut].Data());
655 histo2D->GetYaxis()->SetBinLabel(firsty + icut, labels[icut].Data());
659 //______________________________________________________
660 void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
661 for(Int_t ibit = 32; ibit--; )
662 printf("%d", bitmap & BIT(ibit) ? 1 : 0);
666 //______________________________________________________
667 Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
669 // Check whether ITS area is dead
673 case 2: status = kFALSE; break;
674 case 3: status = kFALSE; break;
675 case 7: status = kFALSE; break;
676 default: status = kTRUE;
681 //______________________________________________________
682 Int_t AliHFEextraCuts::GetITSstatus(const AliVTrack * const track, Int_t layer) const {
684 // Check ITS layer status
687 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
690 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
691 if(esdtrack) esdtrack->GetITSModuleIndexInfo(layer, det, status, xloc, zloc);
697 //______________________________________________________
698 Bool_t AliHFEextraCuts::GetTPCCountSharedMapBitsAboveThreshold(AliVTrack *track){
700 // Checks if number of shared bits is above 1
704 if(track->IsA() == AliESDtrack::Class())
705 shared = &((static_cast<AliESDtrack *>(track))->GetTPCSharedMap());
706 else if(track->IsA() == AliAODTrack::Class())
707 shared = &((static_cast<AliAODTrack *>(track))->GetTPCSharedMap());
711 if(shared->CountBits() >= 2) nsharebit=1;
717 //______________________________________________________
718 UInt_t AliHFEextraCuts::GetTPCncls(AliVTrack *track){
720 // Get Number of findable clusters in the TPC
722 Int_t nClusters = 0; // in case no Information available consider all clusters findable
723 TClass *type = track->IsA();
724 if(TESTBIT(fTPCclusterDef, kFoundIter1)){
725 if(type == AliESDtrack::Class()){
726 AliDebug(2, ("Using def kFoundIter1"));
727 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
728 nClusters = esdtrack->GetTPCNclsIter1();
730 AliDebug(2, ("Number of clusters in iteration 1 not available for AOD tracks"));
732 } else if(TESTBIT(fTPCclusterDef, kCrossedRows)){
733 AliDebug(2, ("Using def kCrossedRows"));
734 nClusters = static_cast<UInt_t>(track->GetTPCClusterInfo(2,1));
735 } else if(TESTBIT(fTPCclusterDef, kFound)){
736 AliDebug(2, ("Using def kFound"));
737 nClusters = track->GetTPCNcls();
740 AliDebug(2, ("Using def kFoundAll"));
741 if(type == AliESDtrack::Class()){
742 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
743 const TBits &clusterTPC = esdtrack->GetTPCClusterMap();
744 nClusters = clusterTPC.CountBits();
747 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
748 const TBits &clusterTPC = aodtrack->GetTPCClusterMap();
749 nClusters = clusterTPC.CountBits();
755 //______________________________________________________
756 Float_t AliHFEextraCuts::GetTPCsharedClustersRatio(AliVTrack *track){
758 // Get fraction of shared TPC clusters
760 Float_t fracClustersTPCShared = 0.0;
761 Int_t nClustersTPC = track->GetTPCNcls();
762 Int_t nClustersTPCShared = 0;
763 TClass *type = track->IsA();
764 if(type == AliESDtrack::Class()){
765 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
766 nClustersTPCShared = esdtrack->GetTPCnclsS();
767 } else if(type == AliAODTrack::Class()){
768 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
769 const TBits &shared = aodtrack->GetTPCSharedMap();
770 nClustersTPCShared = shared.CountBits(0) - shared.CountBits(159);
772 if (nClustersTPC!=0) {
773 fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
775 return fracClustersTPCShared;
778 //______________________________________________________
779 Double_t AliHFEextraCuts::GetTPCclusterRatio(AliVTrack *track){
781 // Get Ratio of found / findable clusters for different definitions
782 // Only implemented for ESD tracks
784 Double_t clusterRatio = 1.; // in case no Information available consider all clusters findable
785 TClass *type = track->IsA();
786 if(TESTBIT(fTPCclusterRatioDef, kCROverFindable)){
787 AliDebug(2, "Using ratio def kCROverFindable");
788 clusterRatio = track->GetTPCNclsF() ? track->GetTPCClusterInfo(2,1)/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // crossed rows/findable
789 } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverCR)){
790 AliDebug(2, "Using ratio def kFoundOverCR");
791 clusterRatio = track->GetTPCClusterInfo(2,0); // found/crossed rows
792 } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindableIter1)){
793 if(type == AliESDtrack::Class()){
794 AliDebug(2, "Using ratio def kFoundOverFindableIter1");
795 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
796 clusterRatio = esdtrack->GetTPCNclsFIter1() ? static_cast<Double_t>(esdtrack->GetTPCNclsIter1())/static_cast<Double_t>(esdtrack->GetTPCNclsFIter1()) : 1.; // crossed
798 AliDebug(2, "Cluster ratio after iteration 1 not possible for AOD tracks");
801 } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindable)) {
802 AliDebug(2, "Using ratio def kFoundOverFindable");
803 clusterRatio = track->GetTPCNclsF() ? static_cast<Double_t>(track->GetTPCNcls())/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // found/findable
806 AliDebug(2, "Using ratio def kFoundAllOverFindable");
807 Int_t allclusters = 0;
808 if(type == AliESDtrack::Class()){
809 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
810 const TBits &clusterTPC = esdtrack->GetTPCClusterMap();
811 allclusters = clusterTPC.CountBits();
814 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
815 const TBits &clusterTPC = aodtrack->GetTPCClusterMap();
816 allclusters = clusterTPC.CountBits();
818 clusterRatio = track->GetTPCNclsF() ? static_cast<Double_t>(allclusters)/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // foundall/findable
823 //___________________________________________________
824 void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z){
826 // Get impact parameter
829 const Double_t kBeampiperadius=3.;
830 Double_t dcaD[2]={-999.,-999.},
831 covD[3]={-999.,-999.,-999.};
832 TClass *type = track->IsA();
833 if(type == AliESDtrack::Class()){
834 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
835 esdtrack->GetImpactParameters(radial, z);
837 else if(type == AliAODTrack::Class()){
840 AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
842 AliDebug(1, "No aod event available\n");
846 //Case ESD track: take copy constructor
847 AliAODTrack *aodtrack = NULL;
848 AliAODTrack *tmptrack = dynamic_cast<AliAODTrack *>(track);
849 if(tmptrack) aodtrack = new AliAODTrack(*tmptrack);
851 AliAODVertex *vtxAODSkip = aodevent->GetPrimaryVertex();
852 if(!vtxAODSkip) return;
853 AliExternalTrackParam etp; etp.CopyFromVTrack(aodtrack);
854 if(etp.PropagateToDCA(vtxAODSkip, aodevent->GetMagneticField(), kBeampiperadius, dcaD, covD)) {
858 //if(vtxAODSkip) delete vtxAODSkip;
859 if(aodtrack) delete aodtrack;
862 //______________________________________________________
863 Bool_t AliHFEextraCuts::IsKinkDaughter(AliVTrack *track){
867 TClass *type = track->IsA();
868 if(type == AliESDtrack::Class()){
869 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
870 if(!esdtrack) return kFALSE;
871 if(esdtrack->GetKinkIndex(0)>0) return kTRUE;
875 else if(type == AliAODTrack::Class()){
876 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
878 AliAODVertex *aodvertex = aodtrack->GetProdVertex();
879 if(!aodvertex) return kFALSE;
880 if(aodvertex->GetType()==AliAODVertex::kKink) return kTRUE;
888 //______________________________________________________
889 Bool_t AliHFEextraCuts::IsKinkMother(AliVTrack *track){
891 // Is kink Mother: only for ESD since need to loop over vertices for AOD
895 TClass *type = track->IsA();
896 if(type == AliESDtrack::Class()){
897 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
898 if(!esdtrack) return kFALSE;
899 if(esdtrack->GetKinkIndex(0)!=0) return kTRUE;
907 //______________________________________________________
908 Float_t AliHFEextraCuts::GetTRDchi(AliVTrack *track){
913 TClass *type = track->IsA();
914 if(type == AliESDtrack::Class()){
915 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
916 ntls = esdtrack->GetTRDntracklets();
917 return ntls ? esdtrack->GetTRDchi2()/ntls : -999;
919 else if(type == AliAODTrack::Class()){
920 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
930 //______________________________________________________
931 Int_t AliHFEextraCuts::GetITSNbOfcls(AliVTrack *track){
933 // Get ITS nb of clusters
935 TClass *type = track->IsA();
936 if(type == AliESDtrack::Class()){
937 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
938 return esdtrack->GetITSclusters(0);
941 else if(type == AliAODTrack::Class()){
942 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
944 return aodtrack->GetITSNcls();
950 //______________________________________________________
951 void AliHFEextraCuts::GetHFEImpactParameters(const AliVTrack * const track, Double_t &dcaxy, Double_t &dcansigmaxy){
953 // Get HFE impact parameter (with recalculated primary vertex)
958 AliDebug(1, "No Input event available\n");
961 TString type = track->IsA()->GetName();
962 const Double_t kBeampiperadius=3.;
963 Double_t dcaD[2]={-999.,-999.},
964 covD[3]={-999.,-999.,-999.};
965 Bool_t isRecalcVertex(kFALSE);
967 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
969 AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
971 AliDebug(1, "No esd event available\n");
975 const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
976 if(!vtxESDSkip) return;
978 //case ESD track: take copy constructor
979 const AliESDtrack *tmptrack = dynamic_cast<const AliESDtrack *>(track);
982 if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
983 vtxESDSkip = RemoveDaughtersFromPrimaryVtx(esdevent, track);
984 isRecalcVertex = kTRUE;
988 AliESDtrack esdtrack(*tmptrack);
989 fMagField = fEvent->GetMagneticField();
990 if(esdtrack.PropagateToDCA(vtxESDSkip, fMagField, kBeampiperadius, dcaD, covD)){
992 if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
994 if(isRecalcVertex) delete vtxESDSkip;
1000 AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
1002 AliDebug(1, "No aod event available\n");
1006 AliAODVertex *vtxAODSkip = aodevent->GetPrimaryVertex();
1007 if(!vtxAODSkip) return;
1009 //Case ESD track: take copy constructor
1010 const AliAODTrack *tmptrack = dynamic_cast<const AliAODTrack *>(track);
1013 if(vtxAODSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
1014 vtxAODSkip = RemoveDaughtersFromPrimaryVtx(aodevent, track);
1015 isRecalcVertex = kTRUE;
1018 AliAODTrack aodtrack(*tmptrack);
1019 AliExternalTrackParam etp; etp.CopyFromVTrack(&aodtrack);
1020 fMagField = aodevent->GetMagneticField();
1021 if(etp.PropagateToDCA(vtxAODSkip,fMagField, kBeampiperadius, dcaD, covD)) {
1023 if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
1025 if(isRecalcVertex) delete vtxAODSkip;
1031 //______________________________________________________
1032 void AliHFEextraCuts::GetHFEImpactParameters(const AliVTrack * const track, Double_t dcaD[2], Double_t covD[3]){
1034 // Get HFE impact parameter (with recalculated primary vertex)
1037 AliDebug(1, "No Input event available\n");
1040 const Double_t kBeampiperadius=3.;
1041 TString type = track->IsA()->GetName();
1042 Bool_t isRecalcVertex(kFALSE);
1044 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
1045 //case of ESD tracks
1046 AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
1048 AliDebug(1, "No esd event available\n");
1052 // Check whether primary vertex is available
1053 const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
1054 if(!vtxESDSkip) return;
1056 const AliESDtrack *tmptrack = dynamic_cast<const AliESDtrack *>(track);
1058 //case ESD track: take copy constructor
1060 if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
1061 vtxESDSkip = RemoveDaughtersFromPrimaryVtx(esdevent, track);
1062 isRecalcVertex = kTRUE;
1065 AliESDtrack esdtrack(*tmptrack);
1066 fMagField = fEvent->GetMagneticField();
1067 esdtrack.PropagateToDCA(vtxESDSkip, fMagField, kBeampiperadius, dcaD, covD);
1069 if(isRecalcVertex) delete vtxESDSkip;
1074 //case of AOD tracks
1075 AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
1077 AliDebug(1, "No aod event available\n");
1081 // Check whether primary vertex is available
1082 AliAODVertex *vtxAODSkip = aodevent->GetPrimaryVertex();
1083 if(!vtxAODSkip) return;
1085 //Case ESD track: take copy constructor
1086 const AliAODTrack *tmptrack = dynamic_cast<const AliAODTrack *>(track);
1089 if(vtxAODSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
1090 vtxAODSkip = RemoveDaughtersFromPrimaryVtx(aodevent, track);
1091 isRecalcVertex = kTRUE;
1094 AliAODTrack aodtrack(*tmptrack);
1095 AliExternalTrackParam etp; etp.CopyFromVTrack(&aodtrack);
1096 fMagField = aodevent->GetMagneticField();
1097 etp.PropagateToDCA(vtxAODSkip, fMagField, kBeampiperadius, dcaD, covD);
1099 if(isRecalcVertex) delete vtxAODSkip;
1105 //______________________________________________________
1106 void AliHFEextraCuts::GetHFEImpactParameterCuts(const AliVTrack * const track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
1108 // Get HFE impact parameter cut(pt dependent)
1111 Double_t pt = track->Pt();
1112 hfeimpactRcut = fIPcutParam[0]+fIPcutParam[1]*exp(fIPcutParam[2]*pt); // abs R cut
1113 hfeimpactnsigmaRcut = fIPcutParam[3]; // sigma cut
1115 //______________________________________________________
1116 void AliHFEextraCuts::GetMaxImpactParameterCutR(const AliVTrack * const track, Double_t &maximpactRcut){
1118 // Get max impact parameter cut r (pt dependent)
1121 Double_t pt = track->Pt();
1123 maximpactRcut = 0.0182 + 0.035/TMath::Power(pt,1.01); // abs R cut
1125 else maximpactRcut = 9999999999.0;
1127 //______________________________________________________
1128 void AliHFEextraCuts::GetTOFsignalDxDz(const AliVTrack * const track, Double_t &tofsignalDx, Double_t &tofsignalDz){
1133 TString type = track->IsA()->GetName();
1134 if(!type.CompareTo("AliESDtrack")){
1135 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
1136 if(!esdtrack) return;
1137 tofsignalDx = esdtrack->GetTOFsignalDx();
1138 tofsignalDz = esdtrack->GetTOFsignalDz();
1143 //______________________________________________________
1144 Bool_t AliHFEextraCuts::CheckITSpattern(const AliVTrack *const track) const {
1146 // Check if every ITS layer, which has a module which is alive, also
1147 // has an ITS cluster
1149 Bool_t patternOK(kTRUE);
1151 for(Int_t ily = 0; ily < 6; ily++){
1152 status = GetITSstatus(track, ily);
1153 if(CheckITSstatus(status)){
1154 // pixel alive, check whether layer has a cluster
1155 if(!TESTBIT(track->GetITSClusterMap(),ily)){
1156 // No cluster even though pixel is alive - reject track
1165 //---------------------------------------------------------------------------
1166 const AliVVertex* AliHFEextraCuts::RemoveDaughtersFromPrimaryVtx(const AliESDEvent * const esdevent, const AliVTrack * const track) {
1168 // This method returns a primary vertex without the daughter tracks of the
1169 // candidate and it recalculates the impact parameters and errors for ESD tracks.
1171 // The output vertex is created with "new".
1174 //const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
1176 AliVertexerTracks vertexer(fEvent->GetMagneticField());
1177 vertexer.SetITSMode();
1178 vertexer.SetMinClusters(4);
1180 skipped[0] = track->GetID();
1181 vertexer.SetSkipTracks(1,skipped);
1183 //diamond constraint
1184 vertexer.SetConstraintOn();
1185 Float_t diamondcovxy[3];
1186 esdevent->GetDiamondCovXY(diamondcovxy);
1187 Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
1188 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1189 AliESDVertex diamond(pos,cov,1.,1);
1190 vertexer.SetVtxStart(&diamond);
1192 const AliVVertex *vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
1197 //---------------------------------------------------------------------------
1198 AliAODVertex* AliHFEextraCuts::RemoveDaughtersFromPrimaryVtx(const AliAODEvent * const aod, const AliVTrack * const track) {
1200 // This method returns a primary vertex without the daughter tracks of the
1201 // candidate and it recalculates the impact parameters and errors for AOD tracks.
1202 // The output vertex is created with "new".
1205 AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
1206 if(!vtxAOD) return 0;
1207 TString title=vtxAOD->GetTitle();
1208 if(!title.Contains("VertexerTracks")) return 0;
1210 AliVertexerTracks vertexer(aod->GetMagneticField());
1212 vertexer.SetITSMode();
1213 vertexer.SetMinClusters(3);
1214 vertexer.SetConstraintOff();
1216 if(title.Contains("WithConstraint")) {
1217 Float_t diamondcovxy[3];
1218 aod->GetDiamondCovXY(diamondcovxy);
1219 Double_t pos[3]={aod->GetDiamondX(),aod->GetDiamondY(),0.};
1220 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
1221 AliESDVertex diamond(pos,cov,1.,1);
1222 vertexer.SetVtxStart(&diamond);
1224 Int_t skipped[2]; for(Int_t i=0;i<2;i++) skipped[i]=-1;
1225 Int_t id = (Int_t)track->GetID();
1226 if(!(id<0)) skipped[0] = id;
1228 /*// exclude tracks with global constrained parameters
1229 Int_t nTracks=aod->GetNumberOfTracks();
1230 for(Int_t i=0; i<nTracks; i++){
1231 t = aod->GetTrack(i);
1232 if(t->TestFilterMask(512)){
1233 id = (Int_t)t->GetID();
1234 skipped[nTrksToSkip++] = id;
1238 vertexer.SetSkipTracks(1,skipped);
1239 AliESDVertex *vtxESDNew = vertexer.FindPrimaryVertex(aod);
1241 if(!vtxESDNew) return 0;
1242 if(vtxESDNew->GetNContributors()<=0) {
1243 delete vtxESDNew; vtxESDNew=NULL;
1247 // convert to AliAODVertex
1248 Double_t pos[3],cov[6],chi2perNDF;
1249 vtxESDNew->GetXYZ(pos); // position
1250 vtxESDNew->GetCovMatrix(cov); //covariance matrix
1251 chi2perNDF = vtxESDNew->GetChi2toNDF();
1252 delete vtxESDNew; vtxESDNew=NULL;
1254 AliAODVertex *vtxAODNew = new AliAODVertex(pos,cov,chi2perNDF);