]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/hfe/AliHFEextraCuts.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
from Ante Bilandzic:
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEextraCuts.cxx
CommitLineData
809a4336 1/**************************************************************************
2* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
5* Contributors are mentioned in the code where appropriate. *
6* *
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**************************************************************************/
50685501 15//
16// Extra cuts implemented by the ALICE Heavy Flavour Electron Group
17// Cuts stored here:
18// - ITS pixels
19// - TPC cluster ratio
20// - TRD tracklets
21//
22// Authors:
23// Markus Fasel <M.Fasel@gsi.de>
24//
3a72645a 25#include <TBits.h>
809a4336 26#include <TClass.h>
27#include <TH1F.h>
28#include <TH2F.h>
29#include <TList.h>
30#include <TString.h>
75d81601 31#include <TMath.h>
809a4336 32
959ea9d8 33#include "AliAODEvent.h"
3a72645a 34#include "AliAODTrack.h"
35#include "AliAODPid.h"
11ff28c5 36#include "AliAODVertex.h"
6400afdd 37#include "AliAODTrack.h"
3a72645a 38#include "AliESDEvent.h"
39#include "AliESDVertex.h"
809a4336 40#include "AliESDtrack.h"
70da6c5a 41#include "AliLog.h"
809a4336 42#include "AliMCParticle.h"
3a72645a 43#include "AliVEvent.h"
44#include "AliVTrack.h"
45#include "AliVParticle.h"
46#include "AliVertexerTracks.h"
47#include "AliVVertex.h"
959ea9d8 48#include "AliExternalTrackParam.h"
809a4336 49
50#include "AliHFEextraCuts.h"
51
52ClassImp(AliHFEextraCuts)
53
76d0b522 54const Int_t AliHFEextraCuts::fgkNQAhistos = 10;
8c1c76e9 55
809a4336 56//______________________________________________________
57AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
58 AliCFCutBase(name, title),
3a72645a 59 fEvent(NULL),
809a4336 60 fCutCorrelation(0),
61 fRequirements(0),
3a72645a 62 fMinNClustersTPC(0),
e17c1f86 63 fMinNClustersTPCPID(0),
809a4336 64 fClusterRatioTPC(0.),
65 fMinTrackletsTRD(0),
cedf0381 66 fMaxChi2TRD(5.0),
11ff28c5 67 fMinNbITScls(0),
aa1417ee 68 fTRDtrackletsExact(0),
809a4336 69 fPixelITS(0),
5cd679b7 70 fDriftITS(0),
e3ae862b 71 fTPCclusterDef(0),
72 fTPCclusterRatioDef(0),
8c1c76e9 73 fFractionTPCShared(-1.0),
63bdf450 74 fOppSideIPcut(kFALSE),
1c051dd4 75 fTOFsignalDx(1.0),
76 fTOFsignalDz(1.0),
63bdf450 77 fMagField(-10000),
96167a04 78 fAODFilterBit(0),
70da6c5a 79 fCheck(kFALSE),
80 fQAlist(0x0) ,
dbe3abbe 81 fDebugLevel(0)
809a4336 82{
83 //
84 // Default Constructor
85 //
11ff28c5 86 //printf("Set the number of min ITS clusters %d\n",(Int_t)fMinNbITScls);
809a4336 87 memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
c2690925 88 memset(fIPcutParam, 0, sizeof(Float_t) * 4);
809a4336 89}
90
91//______________________________________________________
92AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
93 AliCFCutBase(c),
3a72645a 94 fEvent(c.fEvent),
809a4336 95 fCutCorrelation(c.fCutCorrelation),
96 fRequirements(c.fRequirements),
3a72645a 97 fMinNClustersTPC(c.fMinNClustersTPC),
e17c1f86 98 fMinNClustersTPCPID(c.fMinNClustersTPCPID),
809a4336 99 fClusterRatioTPC(c.fClusterRatioTPC),
100 fMinTrackletsTRD(c.fMinTrackletsTRD),
cedf0381 101 fMaxChi2TRD(c.fMaxChi2TRD),
11ff28c5 102 fMinNbITScls(c.fMinNbITScls),
aa1417ee 103 fTRDtrackletsExact(c.fTRDtrackletsExact),
809a4336 104 fPixelITS(c.fPixelITS),
5cd679b7 105 fDriftITS(c.fDriftITS),
e3ae862b 106 fTPCclusterDef(c.fTPCclusterDef),
107 fTPCclusterRatioDef(c.fTPCclusterRatioDef),
8c1c76e9 108 fFractionTPCShared(c.fFractionTPCShared),
63bdf450 109 fOppSideIPcut(c.fOppSideIPcut),
1c051dd4 110 fTOFsignalDx(c.fTOFsignalDx),
111 fTOFsignalDz(c.fTOFsignalDz),
63bdf450 112 fMagField(c.fMagField),
96167a04 113 fAODFilterBit(c.fAODFilterBit),
75d81601 114 fCheck(c.fCheck),
dbe3abbe 115 fQAlist(0x0),
70da6c5a 116 fDebugLevel(0)
809a4336 117{
118 //
119 // Copy constructor
120 // Performs a deep copy
121 //
122 memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
c2690925 123 memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
809a4336 124 if(IsQAOn()){
125 fIsQAOn = kTRUE;
126 fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
bf892a6a 127 if(fQAlist) fQAlist->SetOwner();
809a4336 128 }
129}
130
131//______________________________________________________
132AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
133 //
134 // Assignment operator
135 //
136 if(this != &c){
137 AliCFCutBase::operator=(c);
3a72645a 138 fEvent = c.fEvent;
809a4336 139 fCutCorrelation = c.fCutCorrelation;
140 fRequirements = c.fRequirements;
141 fClusterRatioTPC = c.fClusterRatioTPC;
3a72645a 142 fMinNClustersTPC = c.fMinNClustersTPC;
e17c1f86 143 fMinNClustersTPCPID = c.fMinNClustersTPCPID;
809a4336 144 fMinTrackletsTRD = c.fMinTrackletsTRD;
cedf0381 145 fMaxChi2TRD = c.fMaxChi2TRD;
11ff28c5 146 fMinNbITScls = c.fMinNbITScls;
aa1417ee 147 fTRDtrackletsExact = c.fTRDtrackletsExact;
809a4336 148 fPixelITS = c.fPixelITS;
5cd679b7 149 fDriftITS = c.fDriftITS;
e3ae862b 150 fTPCclusterDef = c.fTPCclusterDef;
151 fTPCclusterRatioDef = c.fTPCclusterRatioDef;
8c1c76e9 152 fFractionTPCShared = c.fFractionTPCShared;
63bdf450 153 fOppSideIPcut = c.fOppSideIPcut;
1c051dd4 154 fTOFsignalDx = c.fTOFsignalDx;
155 fTOFsignalDz = c.fTOFsignalDz;
63bdf450 156 fMagField = c.fMagField;
96167a04 157 fAODFilterBit = c.fAODFilterBit;
75d81601 158 fCheck = c.fCheck;
dbe3abbe 159 fDebugLevel = c.fDebugLevel;
809a4336 160 memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
c2690925 161 memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
809a4336 162 if(IsQAOn()){
163 fIsQAOn = kTRUE;
164 fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
bf892a6a 165 if(fQAlist) fQAlist->SetOwner();
809a4336 166 }else fQAlist = 0x0;
167 }
168 return *this;
169}
170
171//______________________________________________________
172AliHFEextraCuts::~AliHFEextraCuts(){
173 //
174 // Destructor
175 //
3a72645a 176 if(fQAlist) delete fQAlist;
177}
178
179//______________________________________________________
180void AliHFEextraCuts::SetRecEventInfo(const TObject *event){
181 //
182 // Set Virtual event an make a copy
183 //
184 if (!event) {
185 AliError("Pointer to AliVEvent !");
186 return;
187 }
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 !");
191 return ;
809a4336 192 }
3a72645a 193 fEvent = (AliVEvent*) event;
194
809a4336 195}
196
197//______________________________________________________
9eeae5d5 198Bool_t AliHFEextraCuts::IsSelected(TObject *o){
809a4336 199 //
200 // Steering function for the track selection
201 //
11ff28c5 202 TClass *type = o->IsA();
203 AliDebug(2, Form("Object type %s", type->GetName()));
204 if((type == AliESDtrack::Class()) || (type == AliAODTrack::Class())){
3a72645a 205 return CheckRecCuts(dynamic_cast<AliVTrack *>(o));
809a4336 206 }
3a72645a 207 return CheckMCCuts(dynamic_cast<AliVParticle *>(o));
809a4336 208}
209
210//______________________________________________________
3a72645a 211Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
809a4336 212 //
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
217 //
cedf0381 218 AliDebug(1, Form("%s: Called", GetName()));
809a4336 219 ULong64_t survivedCut = 0; // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
3a72645a 220 if(IsQAOn()) FillQAhistosRec(track, kBeforeCuts);
809a4336 221 // Apply cuts
7b4c9210 222 Float_t impactR = -999.;
223 Float_t impactZ = -999.;
3a72645a 224 Double_t hfeimpactR, hfeimpactnsigmaR;
225 Double_t hfeimpactRcut, hfeimpactnsigmaRcut;
9250ffbf 226 Double_t maximpactRcut;
3a72645a 227 GetImpactParameters(track, impactR, impactZ);
63bdf450 228 if(TESTBIT(fRequirements, kMinHFEImpactParamR) || TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)||TESTBIT(fRequirements, kMinHFEImpactParamRcharge)){
3a72645a 229 // Protection for PbPb
230 GetHFEImpactParameterCuts(track, hfeimpactRcut, hfeimpactnsigmaRcut);
231 GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
232 }
11ff28c5 233 Int_t nclsITS = GetITSNbOfcls(track);
e3ae862b 234 UInt_t nclsTPC = GetTPCncls(track);
11ff28c5 235 UInt_t nclsTPCPID = track->GetTPCsignalN();
809a4336 236 // printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
9250ffbf 237 Float_t fractionSharedClustersTPC = GetTPCsharedClustersRatio(track);
e3ae862b 238 Double_t ratioTPC = GetTPCclusterRatio(track);
809a4336 239 UChar_t trdTracklets;
8c1c76e9 240 trdTracklets = track->GetTRDntrackletsPID();
cedf0381 241 Float_t trdchi2=-999.;
242 trdchi2=GetTRDchi(track);
809a4336 243 UChar_t itsPixel = track->GetITSClusterMap();
3a72645a 244 Int_t status1 = GetITSstatus(track, 0);
245 Int_t status2 = GetITSstatus(track, 1);
70da6c5a 246 Bool_t statusL0 = CheckITSstatus(status1);
247 Bool_t statusL1 = CheckITSstatus(status2);
1c051dd4 248 Double_t tofsignalDx = 0.0;
249 Double_t tofsignalDz = 0.0;
250 GetTOFsignalDxDz(track,tofsignalDx,tofsignalDz);
251
8c1c76e9 252 if(TESTBIT(fRequirements, kTPCfractionShared)) {
e17c1f86 253 // cut on max fraction of shared TPC clusters
254 if(TMath::Abs(fractionSharedClustersTPC) < fFractionTPCShared) SETBIT(survivedCut, kTPCfractionShared);
9250ffbf 255 }
809a4336 256 if(TESTBIT(fRequirements, kMinImpactParamR)){
257 // cut on min. Impact Parameter in Radial direction
75d81601 258 if(TMath::Abs(impactR) >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
809a4336 259 }
260 if(TESTBIT(fRequirements, kMinImpactParamZ)){
261 // cut on min. Impact Parameter in Z direction
75d81601 262 if(TMath::Abs(impactZ) >= fImpactParamCut[1]) SETBIT(survivedCut, kMinImpactParamZ);
809a4336 263 }
264 if(TESTBIT(fRequirements, kMaxImpactParamR)){
265 // cut on max. Impact Parameter in Radial direction
75d81601 266 if(TMath::Abs(impactR) <= fImpactParamCut[2]) SETBIT(survivedCut, kMaxImpactParamR);
809a4336 267 }
8c1c76e9 268 if(TESTBIT(fRequirements, kMaxImpactParameterRpar)) {
269 // HFE Impact parameter cut
270 GetMaxImpactParameterCutR(track,maximpactRcut);
271 if(TMath::Abs(impactR) >= maximpactRcut) SETBIT(fRequirements, kMaxImpactParameterRpar);
9250ffbf 272 }
809a4336 273 if(TESTBIT(fRequirements, kMaxImpactParamZ)){
274 // cut on max. Impact Parameter in Z direction
75d81601 275 if(TMath::Abs(impactZ) <= fImpactParamCut[3]) SETBIT(survivedCut, kMaxImpactParamZ);
809a4336 276 }
3a72645a 277 if(TESTBIT(fRequirements, kMinHFEImpactParamR)){
278 // cut on min. HFE Impact Parameter in Radial direction
279 if(TMath::Abs(hfeimpactR) >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamR);
280 }
63bdf450 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();
284
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);
289 }
3a72645a 290 if(TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
291 // cut on max. HFE Impact Parameter n sigma in Radial direction
63bdf450 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);
297 // }
298 // }
299 // else {
300 if(hfeimpactnsigmaRcut>0){
301 if(hfeimpactnsigmaR >= hfeimpactnsigmaRcut) {
11ff28c5 302 SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
63bdf450 303 //printf("1: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
11ff28c5 304 }
305 }
63bdf450 306 else{
307 if(hfeimpactnsigmaR <= hfeimpactnsigmaRcut) {
308 SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
309 //printf("2: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
11ff28c5 310 }
311 }
63bdf450 312 //}
3a72645a 313 }
809a4336 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);
317 }
318 if(TESTBIT(fRequirements, kMinTrackletsTRD)){
319 // cut on minimum number of TRD tracklets
cedf0381 320 AliDebug(1, Form("%s: Min TRD cut: [%d|%d], Require exact number [%s]\n", GetName(), fMinTrackletsTRD, trdTracklets, fTRDtrackletsExact ? "Yes" : "No"));
aa1417ee 321 if(fTRDtrackletsExact){
e17c1f86 322 if(trdTracklets == fMinTrackletsTRD) {
323 SETBIT(survivedCut, kMinTrackletsTRD);
cedf0381 324 AliDebug(1, Form("%s: Track Selected", GetName()));
e17c1f86 325 }
aa1417ee 326 }else{
e17c1f86 327 if(trdTracklets >= fMinTrackletsTRD){
328 SETBIT(survivedCut, kMinTrackletsTRD);
cedf0381 329 AliDebug(1, Form("%s: Track Selected", GetName()));
e17c1f86 330 }
aa1417ee 331 //printf("Min number of tracklets %d\n",fMinTrackletsTRD);
332 }
809a4336 333 }
cedf0381 334
335 if(TESTBIT(fRequirements, kMaxTRDChi2)){
336 // cut on TRD chi2
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));
341 }
342 }
343
11ff28c5 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));
cedf0381 347 AliDebug(1, Form("%s: Min ITS clusters: [%d|%d]\n", GetName(), fMinNbITScls, nclsITS));
11ff28c5 348 if(nclsITS >= ((Int_t)fMinNbITScls)) SETBIT(survivedCut, kMinNbITScls);
349 }
350
3a72645a 351 if(TESTBIT(fRequirements, kMinNClustersTPC)){
11ff28c5 352 // cut on minimum number of TPC tracklets
353 //printf(Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
cedf0381 354 AliDebug(1, Form("%s: Min TPC cut: [%d|%d]\n", GetName(), fMinNClustersTPC, nclsTPC));
3a72645a 355 if(nclsTPC >= fMinNClustersTPC) SETBIT(survivedCut, kMinNClustersTPC);
356 }
e17c1f86 357 if(TESTBIT(fRequirements, kMinNClustersTPCPID)){
cedf0381 358 AliDebug(1, Form("%s: Min TPC PID cut: [%d|%d]\n", GetName(), fMinNClustersTPCPID, nclsTPCPID));
e17c1f86 359 if(nclsTPCPID >= fMinNClustersTPCPID) SETBIT(survivedCut, kMinNClustersTPCPID);
360 }
5cd679b7 361 if(TESTBIT(fRequirements, kDriftITS)){
362 //printf("Require drift\n");
363 switch(fDriftITS){
364 case kFirstD:
365 if(TESTBIT(itsPixel, 2)) SETBIT(survivedCut, kDriftITS);
366 break;
367 default:
368 SETBIT(survivedCut, kDriftITS);
369 break;
370 }
371 }
809a4336 372 if(TESTBIT(fRequirements, kPixelITS)){
373 // cut on ITS pixel layers
cedf0381 374 AliDebug(1, Form("%s: ITS cluster Map: ", GetName()));
70da6c5a 375 //PrintBitMap(itsPixel);
809a4336 376 switch(fPixelITS){
70da6c5a 377 case kFirst:
378 AliDebug(2, "First");
11ff28c5 379 //printf("First\n");
70da6c5a 380 if(TESTBIT(itsPixel, 0))
381 SETBIT(survivedCut, kPixelITS);
382 else
383 if(fCheck && !statusL0)
384 SETBIT(survivedCut, kPixelITS);
dbe3abbe 385 break;
386 case kSecond:
11ff28c5 387 //printf("Second\n");
70da6c5a 388 AliDebug(2, "Second");
389 if(TESTBIT(itsPixel, 1))
390 SETBIT(survivedCut, kPixelITS);
391 else
392 if(fCheck && !statusL1)
393 SETBIT(survivedCut, kPixelITS);
dbe3abbe 394 break;
395 case kBoth:
11ff28c5 396 //printf("Both\n");
70da6c5a 397 AliDebug(2, "Both");
398 if(TESTBIT(itsPixel,0) && TESTBIT(itsPixel,1))
399 SETBIT(survivedCut, kPixelITS);
400 else
401 if(fCheck && !(statusL0 && statusL1))
402 SETBIT(survivedCut, kPixelITS);
403 break;
dbe3abbe 404 case kAny:
11ff28c5 405 //printf("Any\n");
70da6c5a 406 AliDebug(2, "Any");
407 if(TESTBIT(itsPixel, 0) || TESTBIT(itsPixel, 1))
408 SETBIT(survivedCut, kPixelITS);
409 else
410 if(fCheck && !(statusL0 || statusL1))
411 SETBIT(survivedCut, kPixelITS);
dbe3abbe 412 break;
e156c3bb 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);
418 } else {
419 if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0))
420 SETBIT(survivedCut, kPixelITS);
421 }
422 break;
11ff28c5 423 case kNone:
424 // No cut applied, set as survived
425 SETBIT(survivedCut, kPixelITS);
426 break;
70da6c5a 427 default:
11ff28c5 428 // default, defined as no cut applied
70da6c5a 429 AliDebug(2, "None");
11ff28c5 430 SETBIT(survivedCut, kPixelITS);
70da6c5a 431 break;
809a4336 432 }
cedf0381 433 AliDebug(1, Form("%s: Survived Cut: %s\n", GetName(), TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"));
809a4336 434 }
8c1c76e9 435
436 if(TESTBIT(fRequirements, kTOFPID)){
e3ae862b 437 // cut on TOF pid
8c1c76e9 438 if(track->GetStatus() & AliESDtrack::kTOFpid) SETBIT(survivedCut, kTOFPID);
e3ae862b 439 }
c2690925 440
8c1c76e9 441 if(TESTBIT(fRequirements, kTOFmismatch)){
c2690925 442 // cut on TOF mismatch
8c1c76e9 443 if(!(track->GetStatus() & AliESDtrack::kTOFmismatch)) SETBIT(survivedCut, kTOFmismatch);
c2690925 444 }
445
8c1c76e9 446 if(TESTBIT(fRequirements, kTPCPIDCleanUp)){
447 // cut on TPC PID cleanup
8c1c76e9 448 Bool_t fBitsAboveThreshold=GetTPCCountSharedMapBitsAboveThreshold(track);
11ff28c5 449 if((fBitsAboveThreshold==0)&&(nclsTPCPID>80)) SETBIT(survivedCut, kTPCPIDCleanUp);
8c1c76e9 450 }
451
452 if(TESTBIT(fRequirements, kEMCALmatch)){
453 if(track->GetStatus() && AliESDtrack::kEMCALmatch) SETBIT(survivedCut, kEMCALmatch);
454 }
c2690925 455
11ff28c5 456 if(TESTBIT(fRequirements, kRejectKinkDaughter)){
457 //printf("test daughter\n");
458 if(!IsKinkDaughter(track)) SETBIT(survivedCut, kRejectKinkDaughter);
459 //else printf("Found kink daughter\n");
460 }
461
462 if(TESTBIT(fRequirements, kRejectKinkMother)){
463 //printf("test mother\n");
464 if(!IsKinkMother(track)) SETBIT(survivedCut, kRejectKinkMother);
465 //else printf("Found kink mother\n");
466 }
1c051dd4 467 if(TESTBIT(fRequirements, kTOFsignalDxy)){
468 // cut on TOF matching cluster
469 if((TMath::Abs(tofsignalDx) <= fTOFsignalDx) && (TMath::Abs(tofsignalDz) <= fTOFsignalDz)) SETBIT(survivedCut, kTOFsignalDxy);
470 }
cedf0381 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);
474 }
96167a04 475 if(TESTBIT(fRequirements, kAODFilterBit)){
476 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
477 if(aodtrack){
478 Int_t aodfilter = 1 << fAODFilterBit;
479 if(aodtrack->TestFilterBit(aodfilter)) SETBIT(survivedCut, kAODFilterBit);
480 }
481 }
11ff28c5 482
8c1c76e9 483 if(fRequirements == survivedCut){
cedf0381 484 //
809a4336 485 // Track selected
486 //
cedf0381 487 AliDebug(2, Form("%s: Track Survived cuts\n", GetName()));
3a72645a 488 if(IsQAOn()) FillQAhistosRec(track, kAfterCuts);
809a4336 489 return kTRUE;
490 }
cedf0381 491 AliDebug(2, Form("%s: Track cut", GetName()));
809a4336 492 if(IsQAOn()) FillCutCorrelation(survivedCut);
493 return kFALSE;
494}
495
496//______________________________________________________
3a72645a 497Bool_t AliHFEextraCuts::CheckMCCuts(AliVParticle */*track*/) const {
809a4336 498 //
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
503 //
504 return kTRUE; // not yet implemented
505}
506
507//______________________________________________________
3a72645a 508void AliHFEextraCuts::FillQAhistosRec(AliVTrack *track, UInt_t when){
809a4336 509 //
510 // Fill the QA histograms for ESD tracks
511 // Function can be called before cuts or after cut application (second argument)
512 //
75d81601 513 Float_t impactR, impactZ;
3a72645a 514 GetImpactParameters(track, impactR, impactZ);
bf892a6a 515 TH1 *htmp = NULL;
8c1c76e9 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);
809a4336 518 // printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
8c1c76e9 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));
809a4336 522 UChar_t itsPixel = track->GetITSClusterMap();
8c1c76e9 523 TH1 *pixelHist = dynamic_cast<TH1F *>(fQAlist->At(5 + when * fgkNQAhistos));
e3fc062d 524 //Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
bf892a6a 525 if(pixelHist){
526 Double_t firstEntry = 0.5;
527 if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
528 pixelHist->Fill(firstEntry + 3);
529 else{
530 if(itsPixel & BIT(0)){
531 pixelHist->Fill(firstEntry);
532 if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
533 else pixelHist->Fill(firstEntry + 4);
534 }
535 if(itsPixel & BIT(1)){
536 pixelHist->Fill(firstEntry + 1);
537 if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
538 }
809a4336 539 }
540 }
8c1c76e9 541 // Fill histogram with the status bits
542 TH1 *hStatusBits = dynamic_cast<TH1 *>(fQAlist->At(6 + when * fgkNQAhistos));
aa1417ee 543 if(hStatusBits) {
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);
549 }
11ff28c5 550 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(7 + when * fgkNQAhistos)))) htmp->Fill(track->GetTPCsignalN());
cedf0381 551
552 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(8 + when * fgkNQAhistos)))) htmp->Fill(GetTRDchi(track));
76d0b522 553
554 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(9 + when * fgkNQAhistos)))) htmp->Fill(GetITSNbOfcls(track));
809a4336 555}
556
557// //______________________________________________________
558// void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
559// //
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
563// //
564// }
565
566//______________________________________________________
567void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
568 //
569 // Fill cut correlation histograms for tracks that didn't pass cuts
570 //
8c1c76e9 571 TH2 *correlation = dynamic_cast<TH2F *>(fQAlist->At(2 * fgkNQAhistos));
bf892a6a 572 if(correlation){
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);
579 }
809a4336 580 }
581 }
582}
583
584//______________________________________________________
585void AliHFEextraCuts::AddQAHistograms(TList *qaList){
586 //
587 // Add QA histograms
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
592 //
3a72645a 593
809a4336 594 TH1 *histo1D = 0x0;
595 TH2 *histo2D = 0x0;
809a4336 596 TString cutstr[2] = {"before", "after"};
3a72645a 597
598 if(!fQAlist) fQAlist = new TList; // for internal representation, not owner
809a4336 599 for(Int_t icond = 0; icond < 2; icond++){
8c1c76e9 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);
809a4336 602 histo1D->GetXaxis()->SetTitle("Impact Parameter");
603 histo1D->GetYaxis()->SetTitle("Number of Tracks");
8c1c76e9 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);
809a4336 606 histo1D->GetXaxis()->SetTitle("Impact Parameter");
607 histo1D->GetYaxis()->SetTitle("Number of Tracks");
76d0b522 608 qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 100, 0, 1.3)), 2 + icond * fgkNQAhistos);
8c1c76e9 609 fQAlist->AddAt(histo1D, 2 + icond * fgkNQAhistos);
809a4336 610 histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
611 histo1D->GetYaxis()->SetTitle("Number of Tracks");
8c1c76e9 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);
809a4336 614 histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
615 histo1D->GetYaxis()->SetTitle("Number of Tracks");
8c1c76e9 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);
3a72645a 618 histo1D->GetXaxis()->SetTitle("Number of TPC clusters");
619 histo1D->GetYaxis()->SetTitle("Number of Tracks");
8c1c76e9 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);
809a4336 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());
8c1c76e9 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");
cedf0381 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");
76d0b522 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");
cedf0381 646
809a4336 647 }
809a4336 648 // Add cut correlation
8c1c76e9 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"};
809a4336 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());
656 }
809a4336 657}
dbe3abbe 658
659//______________________________________________________
660void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
661 for(Int_t ibit = 32; ibit--; )
662 printf("%d", bitmap & BIT(ibit) ? 1 : 0);
663 printf("\n");
664}
665
666//______________________________________________________
75d81601 667Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
dbe3abbe 668 //
669 // Check whether ITS area is dead
670 //
671 Bool_t status;
672 switch(itsStatus){
673 case 2: status = kFALSE; break;
674 case 3: status = kFALSE; break;
675 case 7: status = kFALSE; break;
676 default: status = kTRUE;
677 }
678 return status;
679}
3a72645a 680
3a72645a 681//______________________________________________________
cedf0381 682Int_t AliHFEextraCuts::GetITSstatus(const AliVTrack * const track, Int_t layer) const {
3a72645a 683 //
684 // Check ITS layer status
685 //
686 Int_t status = 0;
687 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
688 Int_t det;
689 Float_t xloc, zloc;
cedf0381 690 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
bf892a6a 691 if(esdtrack) esdtrack->GetITSModuleIndexInfo(layer, det, status, xloc, zloc);
3a72645a 692 }
693 return status;
694}
695
8c1c76e9 696
697//______________________________________________________
698Bool_t AliHFEextraCuts::GetTPCCountSharedMapBitsAboveThreshold(AliVTrack *track){
699 //
700 // Checks if number of shared bits is above 1
701 //
702 Int_t nsharebit = 1;
11ff28c5 703 const TBits *shared;
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());
708 else
709 return 0;
710
711 if(shared->CountBits() >= 2) nsharebit=1;
712 else nsharebit=0;
713
8c1c76e9 714 return nsharebit;
715}
716
3a72645a 717//______________________________________________________
e3ae862b 718UInt_t AliHFEextraCuts::GetTPCncls(AliVTrack *track){
76d0b522 719 //
720 // Get Number of findable clusters in the TPC
721 //
722 Int_t nClusters = 0; // in case no Information available consider all clusters findable
11ff28c5 723 TClass *type = track->IsA();
724 if(TESTBIT(fTPCclusterDef, kFoundIter1)){
76d0b522 725 if(type == AliESDtrack::Class()){
11ff28c5 726 AliDebug(2, ("Using def kFoundIter1"));
727 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
728 nClusters = esdtrack->GetTPCNclsIter1();
729 } else {
730 AliDebug(2, ("Number of clusters in iteration 1 not available for AOD tracks"));
bf892a6a 731 }
11ff28c5 732 } else if(TESTBIT(fTPCclusterDef, kCrossedRows)){
733 AliDebug(2, ("Using def kCrossedRows"));
734 nClusters = static_cast<UInt_t>(track->GetTPCClusterInfo(2,1));
76d0b522 735 } else if(TESTBIT(fTPCclusterDef, kFound)){
11ff28c5 736 AliDebug(2, ("Using def kFound"));
76d0b522 737 nClusters = track->GetTPCNcls();
738 }
739 else {
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();
745 }
746 else {
747 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
748 const TBits &clusterTPC = aodtrack->GetTPCClusterMap();
749 nClusters = clusterTPC.CountBits();
750 }
3a72645a 751 }
76d0b522 752 return nClusters;
3a72645a 753}
754
9250ffbf 755//______________________________________________________
756Float_t AliHFEextraCuts::GetTPCsharedClustersRatio(AliVTrack *track){
757 //
758 // Get fraction of shared TPC clusters
759 //
760 Float_t fracClustersTPCShared = 0.0;
11ff28c5 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();
76d0b522 767 } else if(type == AliAODTrack::Class()){
11ff28c5 768 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
769 const TBits &shared = aodtrack->GetTPCSharedMap();
770 nClustersTPCShared = shared.CountBits(0) - shared.CountBits(159);
771 }
772 if (nClustersTPC!=0) {
76d0b522 773 fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
9250ffbf 774 }
775 return fracClustersTPCShared;
776}
777
e3ae862b 778//______________________________________________________
779Double_t AliHFEextraCuts::GetTPCclusterRatio(AliVTrack *track){
780 //
781 // Get Ratio of found / findable clusters for different definitions
782 // Only implemented for ESD tracks
783 //
784 Double_t clusterRatio = 1.; // in case no Information available consider all clusters findable
11ff28c5 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
797 } else {
798 AliDebug(2, "Cluster ratio after iteration 1 not possible for AOD tracks");
799 clusterRatio = 1.;
e3ae862b 800 }
76d0b522 801 } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindable)) {
11ff28c5 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
e3ae862b 804 }
76d0b522 805 else {
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();
812 }
813 else {
814 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
815 const TBits &clusterTPC = aodtrack->GetTPCClusterMap();
816 allclusters = clusterTPC.CountBits();
817 }
818 clusterRatio = track->GetTPCNclsF() ? static_cast<Double_t>(allclusters)/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // foundall/findable
819 }
e3ae862b 820 return clusterRatio;
e3ae862b 821
76d0b522 822}
823//___________________________________________________
3a72645a 824void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z){
11ff28c5 825 //
826 // Get impact parameter
827 //
959ea9d8 828
829 const Double_t kBeampiperadius=3.;
830 Double_t dcaD[2]={-999.,-999.},
831 covD[3]={-999.,-999.,-999.};
11ff28c5 832 TClass *type = track->IsA();
833 if(type == AliESDtrack::Class()){
834 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
835 esdtrack->GetImpactParameters(radial, z);
836 }
837 else if(type == AliAODTrack::Class()){
959ea9d8 838
839 //case of AOD tracks
840 AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
841 if(!aodevent) {
842 AliDebug(1, "No aod event available\n");
843 return;
bf892a6a 844 }
959ea9d8 845
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);
850
851 AliAODVertex *vtxAODSkip = aodevent->GetPrimaryVertex();
852 if(!vtxAODSkip) return;
853 AliExternalTrackParam etp; etp.CopyFromVTrack(aodtrack);
63bdf450 854 if(etp.PropagateToDCA(vtxAODSkip, aodevent->GetMagneticField(), kBeampiperadius, dcaD, covD)) {
959ea9d8 855 radial = dcaD[0];
856 z = dcaD[1];
857 }
858 //if(vtxAODSkip) delete vtxAODSkip;
859 if(aodtrack) delete aodtrack;
11ff28c5 860 }
861}
862//______________________________________________________
863Bool_t AliHFEextraCuts::IsKinkDaughter(AliVTrack *track){
864 //
865 // Is kink Daughter
866 //
867 TClass *type = track->IsA();
868 if(type == AliESDtrack::Class()){
869 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
870 if(!esdtrack) return kFALSE;
76d0b522 871 if(esdtrack->GetKinkIndex(0)>0) return kTRUE;
872 else return kFALSE;
11ff28c5 873
874 }
875 else if(type == AliAODTrack::Class()){
876 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
877 if(aodtrack){
878 AliAODVertex *aodvertex = aodtrack->GetProdVertex();
879 if(!aodvertex) return kFALSE;
880 if(aodvertex->GetType()==AliAODVertex::kKink) return kTRUE;
881 else return kFALSE;
882 }
883 else return kFALSE;
884 }
885
886 return kFALSE;
887}
888//______________________________________________________
889Bool_t AliHFEextraCuts::IsKinkMother(AliVTrack *track){
890 //
0e30407a 891 // Is kink Mother: only for ESD since need to loop over vertices for AOD
11ff28c5 892 //
0e30407a 893 //
894
11ff28c5 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;
900 else return kFALSE;
901 }
11ff28c5 902
903 return kFALSE;
904
3a72645a 905}
cedf0381 906
907//______________________________________________________
908Float_t AliHFEextraCuts::GetTRDchi(AliVTrack *track){
909 //
910 // Get TRDchi2
911 //
912 Int_t ntls(0);
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;
918 }
919 else if(type == AliAODTrack::Class()){
920 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
921 if(aodtrack){
922 return 999.;
923 }
924 }
925
926 return 999.;
927
928}
929
11ff28c5 930//______________________________________________________
931Int_t AliHFEextraCuts::GetITSNbOfcls(AliVTrack *track){
932 //
933 // Get ITS nb of clusters
934 //
935 TClass *type = track->IsA();
936 if(type == AliESDtrack::Class()){
937 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
938 return esdtrack->GetITSclusters(0);
3a72645a 939
11ff28c5 940 }
941 else if(type == AliAODTrack::Class()){
942 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
943 if(aodtrack){
944 return aodtrack->GetITSNcls();
945 }
946 }
947
948 return 0;
949}
3a72645a 950//______________________________________________________
63bdf450 951void AliHFEextraCuts::GetHFEImpactParameters(const AliVTrack * const track, Double_t &dcaxy, Double_t &dcansigmaxy){
959ea9d8 952 //
953 // Get HFE impact parameter (with recalculated primary vertex)
954 //
955 dcaxy=0;
956 dcansigmaxy=0;
3a72645a 957 if(!fEvent){
958 AliDebug(1, "No Input event available\n");
959 return;
960 }
3a72645a 961 TString type = track->IsA()->GetName();
959ea9d8 962 const Double_t kBeampiperadius=3.;
963 Double_t dcaD[2]={-999.,-999.},
964 covD[3]={-999.,-999.,-999.};
63bdf450 965 Bool_t isRecalcVertex(kFALSE);
959ea9d8 966
967 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
968 //case of ESD tracks
969 AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
970 if(!esdevent) {
971 AliDebug(1, "No esd event available\n");
972 return;
973 }
974
959ea9d8 975 const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
976 if(!vtxESDSkip) return;
959ea9d8 977
63bdf450 978 //case ESD track: take copy constructor
979 const AliESDtrack *tmptrack = dynamic_cast<const AliESDtrack *>(track);
980 if(tmptrack){
981
982 if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
983 vtxESDSkip = RemoveDaughtersFromPrimaryVtx(esdevent, track);
984 isRecalcVertex = kTRUE;
985 }
986
987 if(vtxESDSkip){
988 AliESDtrack esdtrack(*tmptrack);
989 fMagField = fEvent->GetMagneticField();
990 if(esdtrack.PropagateToDCA(vtxESDSkip, fMagField, kBeampiperadius, dcaD, covD)){
991 dcaxy = dcaD[0];
992 if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
993 }
994 if(isRecalcVertex) delete vtxESDSkip;
995 }
959ea9d8 996 }
959ea9d8 997 }
998 else {
999 //case of AOD tracks
1000 AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
1001 if(!aodevent) {
1002 AliDebug(1, "No aod event available\n");
1003 return;
1004 }
1005
959ea9d8 1006 AliAODVertex *vtxAODSkip = aodevent->GetPrimaryVertex();
1007 if(!vtxAODSkip) return;
63bdf450 1008
1009 //Case ESD track: take copy constructor
1010 const AliAODTrack *tmptrack = dynamic_cast<const AliAODTrack *>(track);
1011 if(tmptrack){
1012
1013 if(vtxAODSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
1014 vtxAODSkip = RemoveDaughtersFromPrimaryVtx(aodevent, track);
1015 isRecalcVertex = kTRUE;
1016 }
1017 if(vtxAODSkip){
1018 AliAODTrack aodtrack(*tmptrack);
1019 AliExternalTrackParam etp; etp.CopyFromVTrack(&aodtrack);
1020 fMagField = aodevent->GetMagneticField();
1021 if(etp.PropagateToDCA(vtxAODSkip,fMagField, kBeampiperadius, dcaD, covD)) {
1022 dcaxy = dcaD[0];
1023 if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
1024 }
1025 if(isRecalcVertex) delete vtxAODSkip;
1026 }
959ea9d8 1027 }
3a72645a 1028 }
3a72645a 1029}
1030
11ff28c5 1031//______________________________________________________
63bdf450 1032void AliHFEextraCuts::GetHFEImpactParameters(const AliVTrack * const track, Double_t dcaD[2], Double_t covD[3]){
11ff28c5 1033 //
1034 // Get HFE impact parameter (with recalculated primary vertex)
1035 //
1036 if(!fEvent){
1037 AliDebug(1, "No Input event available\n");
1038 return;
1039 }
1040 const Double_t kBeampiperadius=3.;
1041 TString type = track->IsA()->GetName();
63bdf450 1042 Bool_t isRecalcVertex(kFALSE);
11ff28c5 1043
959ea9d8 1044 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
1045 //case of ESD tracks
1046 AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
1047 if(!esdevent) {
1048 AliDebug(1, "No esd event available\n");
1049 return;
1050 }
1051
63bdf450 1052 // Check whether primary vertex is available
959ea9d8 1053 const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
1054 if(!vtxESDSkip) return;
959ea9d8 1055
63bdf450 1056 const AliESDtrack *tmptrack = dynamic_cast<const AliESDtrack *>(track);
1057 if(tmptrack){
1058 //case ESD track: take copy constructor
1059
1060 if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
1061 vtxESDSkip = RemoveDaughtersFromPrimaryVtx(esdevent, track);
1062 isRecalcVertex = kTRUE;
1063 }
1064 if(vtxESDSkip){
1065 AliESDtrack esdtrack(*tmptrack);
1066 fMagField = fEvent->GetMagneticField();
1067 esdtrack.PropagateToDCA(vtxESDSkip, fMagField, kBeampiperadius, dcaD, covD);
959ea9d8 1068
63bdf450 1069 if(isRecalcVertex) delete vtxESDSkip;
1070 }
1071 }
959ea9d8 1072 }
1073 else {
1074 //case of AOD tracks
1075 AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
1076 if(!aodevent) {
1077 AliDebug(1, "No aod event available\n");
1078 return;
1079 }
1080
63bdf450 1081 // Check whether primary vertex is available
959ea9d8 1082 AliAODVertex *vtxAODSkip = aodevent->GetPrimaryVertex();
1083 if(!vtxAODSkip) return;
959ea9d8 1084
63bdf450 1085 //Case ESD track: take copy constructor
1086 const AliAODTrack *tmptrack = dynamic_cast<const AliAODTrack *>(track);
1087 if(tmptrack){
1088
1089 if(vtxAODSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
1090 vtxAODSkip = RemoveDaughtersFromPrimaryVtx(aodevent, track);
1091 isRecalcVertex = kTRUE;
1092 }
1093 if(vtxAODSkip){
1094 AliAODTrack aodtrack(*tmptrack);
1095 AliExternalTrackParam etp; etp.CopyFromVTrack(&aodtrack);
1096 fMagField = aodevent->GetMagneticField();
1097 etp.PropagateToDCA(vtxAODSkip, fMagField, kBeampiperadius, dcaD, covD);
1098
1099 if(isRecalcVertex) delete vtxAODSkip;
1100 }
1101 }
11ff28c5 1102 }
1103}
3a72645a 1104
1105//______________________________________________________
63bdf450 1106void AliHFEextraCuts::GetHFEImpactParameterCuts(const AliVTrack * const track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
3a72645a 1107 //
1108 // Get HFE impact parameter cut(pt dependent)
1109 //
1110
63bdf450 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
3a72645a 1114}
9250ffbf 1115//______________________________________________________
63bdf450 1116void AliHFEextraCuts::GetMaxImpactParameterCutR(const AliVTrack * const track, Double_t &maximpactRcut){
9250ffbf 1117 //
1118 // Get max impact parameter cut r (pt dependent)
1119 //
1120
63bdf450 1121 Double_t pt = track->Pt();
1122 if(pt > 0.15) {
1123 maximpactRcut = 0.0182 + 0.035/TMath::Power(pt,1.01); // abs R cut
9250ffbf 1124 }
63bdf450 1125 else maximpactRcut = 9999999999.0;
9250ffbf 1126}
1c051dd4 1127//______________________________________________________
63bdf450 1128void AliHFEextraCuts::GetTOFsignalDxDz(const AliVTrack * const track, Double_t &tofsignalDx, Double_t &tofsignalDz){
1c051dd4 1129 //
1130 // TOF matching
1131 //
1132
1133 TString type = track->IsA()->GetName();
1134 if(!type.CompareTo("AliESDtrack")){
63bdf450 1135 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
1c051dd4 1136 if(!esdtrack) return;
1137 tofsignalDx = esdtrack->GetTOFsignalDx();
1138 tofsignalDz = esdtrack->GetTOFsignalDz();
1139 }
1140
1141}
cedf0381 1142
1143//______________________________________________________
1144Bool_t AliHFEextraCuts::CheckITSpattern(const AliVTrack *const track) const {
1145 //
1146 // Check if every ITS layer, which has a module which is alive, also
1147 // has an ITS cluster
1148 //
1149 Bool_t patternOK(kTRUE);
1150 Int_t status(0);
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
1157 patternOK = kFALSE;
1158 break;
1159 }
1160 }
1161 }
1162 return patternOK;
1163}
959ea9d8 1164
1165//---------------------------------------------------------------------------
63bdf450 1166const AliVVertex* AliHFEextraCuts::RemoveDaughtersFromPrimaryVtx(const AliESDEvent * const esdevent, const AliVTrack * const track) {
959ea9d8 1167 //
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.
1170 //
1171 // The output vertex is created with "new".
1172 //
1173
91e50e2b 1174 //const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
959ea9d8 1175
1176 AliVertexerTracks vertexer(fEvent->GetMagneticField());
1177 vertexer.SetITSMode();
1178 vertexer.SetMinClusters(4);
1179 Int_t skipped[2];
1180 skipped[0] = track->GetID();
1181 vertexer.SetSkipTracks(1,skipped);
1182
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.};
63bdf450 1189 AliESDVertex diamond(pos,cov,1.,1);
1190 vertexer.SetVtxStart(&diamond);
959ea9d8 1191
63bdf450 1192 const AliVVertex *vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
959ea9d8 1193
1194 return vtxESDSkip;
1195}
1196
1197//---------------------------------------------------------------------------
63bdf450 1198AliAODVertex* AliHFEextraCuts::RemoveDaughtersFromPrimaryVtx(const AliAODEvent * const aod, const AliVTrack * const track) {
959ea9d8 1199 //
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".
1203 //
1204
1205 AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
1206 if(!vtxAOD) return 0;
1207 TString title=vtxAOD->GetTitle();
1208 if(!title.Contains("VertexerTracks")) return 0;
1209
63bdf450 1210 AliVertexerTracks vertexer(aod->GetMagneticField());
959ea9d8 1211
63bdf450 1212 vertexer.SetITSMode();
1213 vertexer.SetMinClusters(3);
1214 vertexer.SetConstraintOff();
959ea9d8 1215
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.};
63bdf450 1221 AliESDVertex diamond(pos,cov,1.,1);
1222 vertexer.SetVtxStart(&diamond);
959ea9d8 1223 }
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;
1227
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;
1235 }
1236 }*/
1237
63bdf450 1238 vertexer.SetSkipTracks(1,skipped);
1239 AliESDVertex *vtxESDNew = vertexer.FindPrimaryVertex(aod);
959ea9d8 1240
1241 if(!vtxESDNew) return 0;
1242 if(vtxESDNew->GetNContributors()<=0) {
1243 delete vtxESDNew; vtxESDNew=NULL;
1244 return 0;
1245 }
1246
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;
1253
1254 AliAODVertex *vtxAODNew = new AliAODVertex(pos,cov,chi2perNDF);
1255
1256 return vtxAODNew;
1257}