]> 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
updated for systematic study
[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),
bbeea05b 79 fListKinkMothers(1000),
80 fNumberKinkMothers(0),
70da6c5a 81 fCheck(kFALSE),
82 fQAlist(0x0) ,
dbe3abbe 83 fDebugLevel(0)
809a4336 84{
85 //
86 // Default Constructor
87 //
11ff28c5 88 //printf("Set the number of min ITS clusters %d\n",(Int_t)fMinNbITScls);
809a4336 89 memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
c2690925 90 memset(fIPcutParam, 0, sizeof(Float_t) * 4);
809a4336 91}
92
93//______________________________________________________
94AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
95 AliCFCutBase(c),
3a72645a 96 fEvent(c.fEvent),
809a4336 97 fCutCorrelation(c.fCutCorrelation),
98 fRequirements(c.fRequirements),
3a72645a 99 fMinNClustersTPC(c.fMinNClustersTPC),
e17c1f86 100 fMinNClustersTPCPID(c.fMinNClustersTPCPID),
809a4336 101 fClusterRatioTPC(c.fClusterRatioTPC),
102 fMinTrackletsTRD(c.fMinTrackletsTRD),
cedf0381 103 fMaxChi2TRD(c.fMaxChi2TRD),
11ff28c5 104 fMinNbITScls(c.fMinNbITScls),
aa1417ee 105 fTRDtrackletsExact(c.fTRDtrackletsExact),
809a4336 106 fPixelITS(c.fPixelITS),
5cd679b7 107 fDriftITS(c.fDriftITS),
e3ae862b 108 fTPCclusterDef(c.fTPCclusterDef),
109 fTPCclusterRatioDef(c.fTPCclusterRatioDef),
8c1c76e9 110 fFractionTPCShared(c.fFractionTPCShared),
63bdf450 111 fOppSideIPcut(c.fOppSideIPcut),
1c051dd4 112 fTOFsignalDx(c.fTOFsignalDx),
113 fTOFsignalDz(c.fTOFsignalDz),
63bdf450 114 fMagField(c.fMagField),
96167a04 115 fAODFilterBit(c.fAODFilterBit),
bbeea05b 116 fListKinkMothers(c.fListKinkMothers),
117 fNumberKinkMothers(c.fNumberKinkMothers),
75d81601 118 fCheck(c.fCheck),
dbe3abbe 119 fQAlist(0x0),
70da6c5a 120 fDebugLevel(0)
809a4336 121{
122 //
123 // Copy constructor
124 // Performs a deep copy
125 //
126 memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
c2690925 127 memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
809a4336 128 if(IsQAOn()){
129 fIsQAOn = kTRUE;
130 fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
bf892a6a 131 if(fQAlist) fQAlist->SetOwner();
809a4336 132 }
133}
134
135//______________________________________________________
136AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
137 //
138 // Assignment operator
139 //
140 if(this != &c){
141 AliCFCutBase::operator=(c);
3a72645a 142 fEvent = c.fEvent;
809a4336 143 fCutCorrelation = c.fCutCorrelation;
144 fRequirements = c.fRequirements;
145 fClusterRatioTPC = c.fClusterRatioTPC;
3a72645a 146 fMinNClustersTPC = c.fMinNClustersTPC;
e17c1f86 147 fMinNClustersTPCPID = c.fMinNClustersTPCPID;
809a4336 148 fMinTrackletsTRD = c.fMinTrackletsTRD;
cedf0381 149 fMaxChi2TRD = c.fMaxChi2TRD;
11ff28c5 150 fMinNbITScls = c.fMinNbITScls;
aa1417ee 151 fTRDtrackletsExact = c.fTRDtrackletsExact;
809a4336 152 fPixelITS = c.fPixelITS;
5cd679b7 153 fDriftITS = c.fDriftITS;
e3ae862b 154 fTPCclusterDef = c.fTPCclusterDef;
155 fTPCclusterRatioDef = c.fTPCclusterRatioDef;
8c1c76e9 156 fFractionTPCShared = c.fFractionTPCShared;
63bdf450 157 fOppSideIPcut = c.fOppSideIPcut;
1c051dd4 158 fTOFsignalDx = c.fTOFsignalDx;
159 fTOFsignalDz = c.fTOFsignalDz;
63bdf450 160 fMagField = c.fMagField;
96167a04 161 fAODFilterBit = c.fAODFilterBit;
bbeea05b 162 fListKinkMothers = c.fListKinkMothers;
163 fNumberKinkMothers = c.fNumberKinkMothers;
75d81601 164 fCheck = c.fCheck;
dbe3abbe 165 fDebugLevel = c.fDebugLevel;
809a4336 166 memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
c2690925 167 memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
809a4336 168 if(IsQAOn()){
169 fIsQAOn = kTRUE;
170 fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
bf892a6a 171 if(fQAlist) fQAlist->SetOwner();
809a4336 172 }else fQAlist = 0x0;
173 }
174 return *this;
175}
176
177//______________________________________________________
178AliHFEextraCuts::~AliHFEextraCuts(){
179 //
180 // Destructor
181 //
3a72645a 182 if(fQAlist) delete fQAlist;
183}
184
185//______________________________________________________
186void AliHFEextraCuts::SetRecEventInfo(const TObject *event){
187 //
188 // Set Virtual event an make a copy
189 //
190 if (!event) {
191 AliError("Pointer to AliVEvent !");
192 return;
193 }
194 TString className(event->ClassName());
195 if (! (className.CompareTo("AliESDEvent")==0 || className.CompareTo("AliAODEvent")==0)) {
196 AliError("argument must point to an AliESDEvent or AliAODEvent !");
197 return ;
809a4336 198 }
3a72645a 199 fEvent = (AliVEvent*) event;
bbeea05b 200
201 AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
202 if(aodevent){
203 // Initialize lookup for kink mother rejection
204 if(aodevent->GetNumberOfVertices() > fListKinkMothers.GetSize()) fListKinkMothers.Set(aodevent->GetNumberOfVertices());
205 fNumberKinkMothers = 0;
206 for(Int_t ivtx = 0; ivtx < aodevent->GetNumberOfVertices(); ivtx++){
207 AliAODVertex *aodvtx = aodevent->GetVertex(ivtx);
208 if(aodvtx->GetType() != AliAODVertex::kKink) continue;
209 AliAODTrack *mother = (AliAODTrack *) aodvtx->GetParent();
210 if(!mother) continue;
211 Int_t idmother = mother->GetID();
212 fListKinkMothers[fNumberKinkMothers++] = idmother;
213 }
214 }
809a4336 215}
216
217//______________________________________________________
9eeae5d5 218Bool_t AliHFEextraCuts::IsSelected(TObject *o){
809a4336 219 //
220 // Steering function for the track selection
221 //
11ff28c5 222 TClass *type = o->IsA();
223 AliDebug(2, Form("Object type %s", type->GetName()));
224 if((type == AliESDtrack::Class()) || (type == AliAODTrack::Class())){
3a72645a 225 return CheckRecCuts(dynamic_cast<AliVTrack *>(o));
809a4336 226 }
3a72645a 227 return CheckMCCuts(dynamic_cast<AliVParticle *>(o));
809a4336 228}
229
230//______________________________________________________
3a72645a 231Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
809a4336 232 //
233 // Checks cuts on reconstructed tracks
234 // returns true if track is selected
235 // QA histograms are filled before track selection and for
236 // selected tracks after track selection
237 //
cedf0381 238 AliDebug(1, Form("%s: Called", GetName()));
809a4336 239 ULong64_t survivedCut = 0; // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
3a72645a 240 if(IsQAOn()) FillQAhistosRec(track, kBeforeCuts);
809a4336 241 // Apply cuts
7b4c9210 242 Float_t impactR = -999.;
243 Float_t impactZ = -999.;
3a72645a 244 Double_t hfeimpactR, hfeimpactnsigmaR;
245 Double_t hfeimpactRcut, hfeimpactnsigmaRcut;
9250ffbf 246 Double_t maximpactRcut;
3a72645a 247 GetImpactParameters(track, impactR, impactZ);
63bdf450 248 if(TESTBIT(fRequirements, kMinHFEImpactParamR) || TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)||TESTBIT(fRequirements, kMinHFEImpactParamRcharge)){
3a72645a 249 // Protection for PbPb
250 GetHFEImpactParameterCuts(track, hfeimpactRcut, hfeimpactnsigmaRcut);
251 GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
252 }
11ff28c5 253 Int_t nclsITS = GetITSNbOfcls(track);
e3ae862b 254 UInt_t nclsTPC = GetTPCncls(track);
11ff28c5 255 UInt_t nclsTPCPID = track->GetTPCsignalN();
809a4336 256 // printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
9250ffbf 257 Float_t fractionSharedClustersTPC = GetTPCsharedClustersRatio(track);
e3ae862b 258 Double_t ratioTPC = GetTPCclusterRatio(track);
809a4336 259 UChar_t trdTracklets;
8c1c76e9 260 trdTracklets = track->GetTRDntrackletsPID();
cedf0381 261 Float_t trdchi2=-999.;
262 trdchi2=GetTRDchi(track);
809a4336 263 UChar_t itsPixel = track->GetITSClusterMap();
3a72645a 264 Int_t status1 = GetITSstatus(track, 0);
265 Int_t status2 = GetITSstatus(track, 1);
70da6c5a 266 Bool_t statusL0 = CheckITSstatus(status1);
267 Bool_t statusL1 = CheckITSstatus(status2);
1c051dd4 268 Double_t tofsignalDx = 0.0;
269 Double_t tofsignalDz = 0.0;
270 GetTOFsignalDxDz(track,tofsignalDx,tofsignalDz);
271
8c1c76e9 272 if(TESTBIT(fRequirements, kTPCfractionShared)) {
e17c1f86 273 // cut on max fraction of shared TPC clusters
274 if(TMath::Abs(fractionSharedClustersTPC) < fFractionTPCShared) SETBIT(survivedCut, kTPCfractionShared);
9250ffbf 275 }
809a4336 276 if(TESTBIT(fRequirements, kMinImpactParamR)){
277 // cut on min. Impact Parameter in Radial direction
75d81601 278 if(TMath::Abs(impactR) >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
809a4336 279 }
280 if(TESTBIT(fRequirements, kMinImpactParamZ)){
281 // cut on min. Impact Parameter in Z direction
75d81601 282 if(TMath::Abs(impactZ) >= fImpactParamCut[1]) SETBIT(survivedCut, kMinImpactParamZ);
809a4336 283 }
284 if(TESTBIT(fRequirements, kMaxImpactParamR)){
285 // cut on max. Impact Parameter in Radial direction
75d81601 286 if(TMath::Abs(impactR) <= fImpactParamCut[2]) SETBIT(survivedCut, kMaxImpactParamR);
809a4336 287 }
8c1c76e9 288 if(TESTBIT(fRequirements, kMaxImpactParameterRpar)) {
289 // HFE Impact parameter cut
290 GetMaxImpactParameterCutR(track,maximpactRcut);
291 if(TMath::Abs(impactR) >= maximpactRcut) SETBIT(fRequirements, kMaxImpactParameterRpar);
9250ffbf 292 }
809a4336 293 if(TESTBIT(fRequirements, kMaxImpactParamZ)){
294 // cut on max. Impact Parameter in Z direction
75d81601 295 if(TMath::Abs(impactZ) <= fImpactParamCut[3]) SETBIT(survivedCut, kMaxImpactParamZ);
809a4336 296 }
3a72645a 297 if(TESTBIT(fRequirements, kMinHFEImpactParamR)){
298 // cut on min. HFE Impact Parameter in Radial direction
299 if(TMath::Abs(hfeimpactR) >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamR);
300 }
63bdf450 301 if(TESTBIT(fRequirements, kMinHFEImpactParamRcharge)){
302 // cut on min. HFE Impact Parameter in Radial direction, multiplied by particle charge
303 Double_t charge = (Double_t)track->Charge();
304
305 if(fMagField < 0)charge *= -1.;//the IP distribution side to be chosen depends on magnetic field polarization
306 if(fOppSideIPcut)charge*=-1.;//in case we choose to select electrons from the side of the photon peak
307 Double_t hfeimpactRcharge = hfeimpactR*charge;//switch selected side of the peak
308 if(hfeimpactRcharge >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamRcharge);
309 }
3a72645a 310 if(TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
311 // cut on max. HFE Impact Parameter n sigma in Radial direction
63bdf450 312 // if(fAbsHFEImpactParamNsigmaR) {
313 // //if((TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) && (TMath::Abs(hfeimpactnsigmaR) < 8)) { //mj debug
314 // if(TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) {
315 // SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
316 // // printf("0: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
317 // }
318 // }
319 // else {
320 if(hfeimpactnsigmaRcut>0){
321 if(hfeimpactnsigmaR >= hfeimpactnsigmaRcut) {
11ff28c5 322 SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
63bdf450 323 //printf("1: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
11ff28c5 324 }
325 }
63bdf450 326 else{
327 if(hfeimpactnsigmaR <= hfeimpactnsigmaRcut) {
328 SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
329 //printf("2: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
11ff28c5 330 }
331 }
63bdf450 332 //}
3a72645a 333 }
809a4336 334 if(TESTBIT(fRequirements, kClusterRatioTPC)){
335 // cut on min ratio of found TPC clusters vs findable TPC clusters
336 if(ratioTPC >= fClusterRatioTPC) SETBIT(survivedCut, kClusterRatioTPC);
337 }
338 if(TESTBIT(fRequirements, kMinTrackletsTRD)){
339 // cut on minimum number of TRD tracklets
cedf0381 340 AliDebug(1, Form("%s: Min TRD cut: [%d|%d], Require exact number [%s]\n", GetName(), fMinTrackletsTRD, trdTracklets, fTRDtrackletsExact ? "Yes" : "No"));
aa1417ee 341 if(fTRDtrackletsExact){
e17c1f86 342 if(trdTracklets == fMinTrackletsTRD) {
343 SETBIT(survivedCut, kMinTrackletsTRD);
cedf0381 344 AliDebug(1, Form("%s: Track Selected", GetName()));
e17c1f86 345 }
aa1417ee 346 }else{
e17c1f86 347 if(trdTracklets >= fMinTrackletsTRD){
348 SETBIT(survivedCut, kMinTrackletsTRD);
cedf0381 349 AliDebug(1, Form("%s: Track Selected", GetName()));
e17c1f86 350 }
aa1417ee 351 //printf("Min number of tracklets %d\n",fMinTrackletsTRD);
352 }
809a4336 353 }
cedf0381 354
355 if(TESTBIT(fRequirements, kMaxTRDChi2)){
356 // cut on TRD chi2
357 AliDebug(1, Form("%s: Cut on TRD chi2: [%f|%f]\n", GetName(),fMaxChi2TRD, trdchi2));
358 if(trdchi2 < fMaxChi2TRD) {
359 SETBIT(survivedCut, kMaxTRDChi2);
360 AliDebug(1,Form("%s: survived %f\n",GetName(),trdchi2));
361 }
362 }
363
11ff28c5 364 if(TESTBIT(fRequirements, kMinNbITScls)){
365 // cut on minimum number of ITS clusters
366 //printf(Form("Min ITS clusters: [%d|%d]\n", (Int_t)fMinNbITScls, nclsITS));
cedf0381 367 AliDebug(1, Form("%s: Min ITS clusters: [%d|%d]\n", GetName(), fMinNbITScls, nclsITS));
11ff28c5 368 if(nclsITS >= ((Int_t)fMinNbITScls)) SETBIT(survivedCut, kMinNbITScls);
369 }
370
3a72645a 371 if(TESTBIT(fRequirements, kMinNClustersTPC)){
11ff28c5 372 // cut on minimum number of TPC tracklets
373 //printf(Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
cedf0381 374 AliDebug(1, Form("%s: Min TPC cut: [%d|%d]\n", GetName(), fMinNClustersTPC, nclsTPC));
3a72645a 375 if(nclsTPC >= fMinNClustersTPC) SETBIT(survivedCut, kMinNClustersTPC);
376 }
e17c1f86 377 if(TESTBIT(fRequirements, kMinNClustersTPCPID)){
cedf0381 378 AliDebug(1, Form("%s: Min TPC PID cut: [%d|%d]\n", GetName(), fMinNClustersTPCPID, nclsTPCPID));
e17c1f86 379 if(nclsTPCPID >= fMinNClustersTPCPID) SETBIT(survivedCut, kMinNClustersTPCPID);
380 }
5cd679b7 381 if(TESTBIT(fRequirements, kDriftITS)){
382 //printf("Require drift\n");
383 switch(fDriftITS){
384 case kFirstD:
385 if(TESTBIT(itsPixel, 2)) SETBIT(survivedCut, kDriftITS);
386 break;
387 default:
388 SETBIT(survivedCut, kDriftITS);
389 break;
390 }
391 }
809a4336 392 if(TESTBIT(fRequirements, kPixelITS)){
393 // cut on ITS pixel layers
cedf0381 394 AliDebug(1, Form("%s: ITS cluster Map: ", GetName()));
70da6c5a 395 //PrintBitMap(itsPixel);
809a4336 396 switch(fPixelITS){
70da6c5a 397 case kFirst:
398 AliDebug(2, "First");
11ff28c5 399 //printf("First\n");
70da6c5a 400 if(TESTBIT(itsPixel, 0))
401 SETBIT(survivedCut, kPixelITS);
402 else
403 if(fCheck && !statusL0)
404 SETBIT(survivedCut, kPixelITS);
dbe3abbe 405 break;
406 case kSecond:
11ff28c5 407 //printf("Second\n");
70da6c5a 408 AliDebug(2, "Second");
409 if(TESTBIT(itsPixel, 1))
410 SETBIT(survivedCut, kPixelITS);
411 else
412 if(fCheck && !statusL1)
413 SETBIT(survivedCut, kPixelITS);
dbe3abbe 414 break;
415 case kBoth:
11ff28c5 416 //printf("Both\n");
70da6c5a 417 AliDebug(2, "Both");
418 if(TESTBIT(itsPixel,0) && TESTBIT(itsPixel,1))
419 SETBIT(survivedCut, kPixelITS);
420 else
421 if(fCheck && !(statusL0 && statusL1))
422 SETBIT(survivedCut, kPixelITS);
423 break;
dbe3abbe 424 case kAny:
11ff28c5 425 //printf("Any\n");
70da6c5a 426 AliDebug(2, "Any");
427 if(TESTBIT(itsPixel, 0) || TESTBIT(itsPixel, 1))
428 SETBIT(survivedCut, kPixelITS);
429 else
430 if(fCheck && !(statusL0 || statusL1))
431 SETBIT(survivedCut, kPixelITS);
dbe3abbe 432 break;
e156c3bb 433 case kExclusiveSecond:
434 AliDebug(2, "Exlusive second");
435 if(fCheck){ // Cut out tracks which pass a dead ITS Layer 0
436 if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0) && statusL0)
437 SETBIT(survivedCut, kPixelITS);
438 } else {
439 if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0))
440 SETBIT(survivedCut, kPixelITS);
441 }
442 break;
11ff28c5 443 case kNone:
444 // No cut applied, set as survived
445 SETBIT(survivedCut, kPixelITS);
446 break;
70da6c5a 447 default:
11ff28c5 448 // default, defined as no cut applied
70da6c5a 449 AliDebug(2, "None");
11ff28c5 450 SETBIT(survivedCut, kPixelITS);
70da6c5a 451 break;
809a4336 452 }
cedf0381 453 AliDebug(1, Form("%s: Survived Cut: %s\n", GetName(), TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"));
809a4336 454 }
8c1c76e9 455
456 if(TESTBIT(fRequirements, kTOFPID)){
e3ae862b 457 // cut on TOF pid
8c1c76e9 458 if(track->GetStatus() & AliESDtrack::kTOFpid) SETBIT(survivedCut, kTOFPID);
e3ae862b 459 }
c2690925 460
8c1c76e9 461 if(TESTBIT(fRequirements, kTOFmismatch)){
c2690925 462 // cut on TOF mismatch
8c1c76e9 463 if(!(track->GetStatus() & AliESDtrack::kTOFmismatch)) SETBIT(survivedCut, kTOFmismatch);
c2690925 464 }
465
2e4ed823 466 if(TESTBIT(fRequirements, kTOFlabel)){
467 if(MatchTOFlabel(track)) SETBIT(survivedCut, kTOFlabel);
468 }
469
8c1c76e9 470 if(TESTBIT(fRequirements, kTPCPIDCleanUp)){
471 // cut on TPC PID cleanup
8c1c76e9 472 Bool_t fBitsAboveThreshold=GetTPCCountSharedMapBitsAboveThreshold(track);
11ff28c5 473 if((fBitsAboveThreshold==0)&&(nclsTPCPID>80)) SETBIT(survivedCut, kTPCPIDCleanUp);
8c1c76e9 474 }
475
476 if(TESTBIT(fRequirements, kEMCALmatch)){
477 if(track->GetStatus() && AliESDtrack::kEMCALmatch) SETBIT(survivedCut, kEMCALmatch);
478 }
c2690925 479
11ff28c5 480 if(TESTBIT(fRequirements, kRejectKinkDaughter)){
481 //printf("test daughter\n");
482 if(!IsKinkDaughter(track)) SETBIT(survivedCut, kRejectKinkDaughter);
483 //else printf("Found kink daughter\n");
484 }
485
486 if(TESTBIT(fRequirements, kRejectKinkMother)){
487 //printf("test mother\n");
488 if(!IsKinkMother(track)) SETBIT(survivedCut, kRejectKinkMother);
489 //else printf("Found kink mother\n");
bbeea05b 490 }
1c051dd4 491 if(TESTBIT(fRequirements, kTOFsignalDxy)){
492 // cut on TOF matching cluster
493 if((TMath::Abs(tofsignalDx) <= fTOFsignalDx) && (TMath::Abs(tofsignalDz) <= fTOFsignalDz)) SETBIT(survivedCut, kTOFsignalDxy);
494 }
cedf0381 495 if(TESTBIT(fRequirements, kITSpattern)){
496 // cut on ITS pattern (every layer with a working ITS module must have an ITS cluster)
497 if(CheckITSpattern(track)) SETBIT(survivedCut, kITSpattern);
498 }
96167a04 499 if(TESTBIT(fRequirements, kAODFilterBit)){
500 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
501 if(aodtrack){
502 Int_t aodfilter = 1 << fAODFilterBit;
503 if(aodtrack->TestFilterBit(aodfilter)) SETBIT(survivedCut, kAODFilterBit);
504 }
505 }
11ff28c5 506
8c1c76e9 507 if(fRequirements == survivedCut){
cedf0381 508 //
809a4336 509 // Track selected
510 //
cedf0381 511 AliDebug(2, Form("%s: Track Survived cuts\n", GetName()));
3a72645a 512 if(IsQAOn()) FillQAhistosRec(track, kAfterCuts);
809a4336 513 return kTRUE;
514 }
cedf0381 515 AliDebug(2, Form("%s: Track cut", GetName()));
809a4336 516 if(IsQAOn()) FillCutCorrelation(survivedCut);
517 return kFALSE;
518}
519
520//______________________________________________________
3a72645a 521Bool_t AliHFEextraCuts::CheckMCCuts(AliVParticle */*track*/) const {
809a4336 522 //
523 // Checks cuts on Monte Carlo tracks
524 // returns true if track is selected
525 // QA histograms are filled before track selection and for
526 // selected tracks after track selection
527 //
528 return kTRUE; // not yet implemented
529}
530
531//______________________________________________________
3a72645a 532void AliHFEextraCuts::FillQAhistosRec(AliVTrack *track, UInt_t when){
809a4336 533 //
534 // Fill the QA histograms for ESD tracks
535 // Function can be called before cuts or after cut application (second argument)
536 //
75d81601 537 Float_t impactR, impactZ;
3a72645a 538 GetImpactParameters(track, impactR, impactZ);
bf892a6a 539 TH1 *htmp = NULL;
8c1c76e9 540 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(0 + when * fgkNQAhistos)))) htmp->Fill(impactR);
541 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(1 + when * fgkNQAhistos)))) htmp->Fill(impactZ);
809a4336 542 // printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
8c1c76e9 543 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(2 + when * fgkNQAhistos)))) htmp->Fill(GetTPCclusterRatio(track));
544 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(3 + when * fgkNQAhistos)))) htmp->Fill(track->GetTRDntrackletsPID());
545 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(4 + when * fgkNQAhistos)))) htmp->Fill(GetTPCncls(track));
809a4336 546 UChar_t itsPixel = track->GetITSClusterMap();
8c1c76e9 547 TH1 *pixelHist = dynamic_cast<TH1F *>(fQAlist->At(5 + when * fgkNQAhistos));
e3fc062d 548 //Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
bf892a6a 549 if(pixelHist){
550 Double_t firstEntry = 0.5;
551 if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
552 pixelHist->Fill(firstEntry + 3);
553 else{
554 if(itsPixel & BIT(0)){
555 pixelHist->Fill(firstEntry);
556 if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
557 else pixelHist->Fill(firstEntry + 4);
558 }
559 if(itsPixel & BIT(1)){
560 pixelHist->Fill(firstEntry + 1);
561 if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
562 }
809a4336 563 }
564 }
8c1c76e9 565 // Fill histogram with the status bits
566 TH1 *hStatusBits = dynamic_cast<TH1 *>(fQAlist->At(6 + when * fgkNQAhistos));
aa1417ee 567 if(hStatusBits) {
568 hStatusBits->Fill(0); // Fill first bin with all tracks
569 if(track->GetStatus() && AliESDtrack::kTOFpid) hStatusBits->Fill(1);
570 if(!(track->GetStatus() && AliESDtrack::kTOFmismatch)) hStatusBits->Fill(2);
571 if(track->GetStatus() && AliESDtrack::kEMCALmatch) hStatusBits->Fill(3);
572 if(GetTPCCountSharedMapBitsAboveThreshold(track)==0) hStatusBits->Fill(4);
573 }
11ff28c5 574 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(7 + when * fgkNQAhistos)))) htmp->Fill(track->GetTPCsignalN());
cedf0381 575
576 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(8 + when * fgkNQAhistos)))) htmp->Fill(GetTRDchi(track));
76d0b522 577
578 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(9 + when * fgkNQAhistos)))) htmp->Fill(GetITSNbOfcls(track));
809a4336 579}
580
581// //______________________________________________________
582// void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
583// //
584// // Fill the QA histograms for MC tracks
585// // Function can be called before cuts or after cut application (second argument)
586// // Not yet implemented
587// //
588// }
589
590//______________________________________________________
591void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
592 //
593 // Fill cut correlation histograms for tracks that didn't pass cuts
594 //
8c1c76e9 595 TH2 *correlation = dynamic_cast<TH2F *>(fQAlist->At(2 * fgkNQAhistos));
bf892a6a 596 if(correlation){
597 for(Int_t icut = 0; icut < kNcuts; icut++){
598 if(!TESTBIT(fRequirements, icut)) continue;
599 for(Int_t jcut = icut; jcut < kNcuts; jcut++){
600 if(!TESTBIT(fRequirements, jcut)) continue;
601 if(TESTBIT(survivedCut, icut) && TESTBIT(survivedCut, jcut))
602 correlation->Fill(icut, jcut);
603 }
809a4336 604 }
605 }
606}
607
608//______________________________________________________
609void AliHFEextraCuts::AddQAHistograms(TList *qaList){
610 //
611 // Add QA histograms
612 // For each cut a histogram before and after track cut is created
613 // Histos before respectively after cut are stored in different lists
614 // Additionally a histogram with the cut correlation is created and stored
615 // in the top directory
616 //
3a72645a 617
809a4336 618 TH1 *histo1D = 0x0;
619 TH2 *histo2D = 0x0;
809a4336 620 TString cutstr[2] = {"before", "after"};
3a72645a 621
622 if(!fQAlist) fQAlist = new TList; // for internal representation, not owner
809a4336 623 for(Int_t icond = 0; icond < 2; icond++){
8c1c76e9 624 qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamR%s",GetName(),cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0 + icond * fgkNQAhistos);
625 fQAlist->AddAt(histo1D, 0 + icond * fgkNQAhistos);
809a4336 626 histo1D->GetXaxis()->SetTitle("Impact Parameter");
627 histo1D->GetYaxis()->SetTitle("Number of Tracks");
8c1c76e9 628 qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamZ%s",GetName(),cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1 + icond * fgkNQAhistos);
629 fQAlist->AddAt(histo1D, 1 + icond * fgkNQAhistos);
809a4336 630 histo1D->GetXaxis()->SetTitle("Impact Parameter");
631 histo1D->GetYaxis()->SetTitle("Number of Tracks");
76d0b522 632 qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 100, 0, 1.3)), 2 + icond * fgkNQAhistos);
8c1c76e9 633 fQAlist->AddAt(histo1D, 2 + icond * fgkNQAhistos);
809a4336 634 histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
635 histo1D->GetYaxis()->SetTitle("Number of Tracks");
8c1c76e9 636 qaList->AddAt((histo1D = new TH1F(Form("%s_trdTracklets%s",GetName(),cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3 + icond * fgkNQAhistos);
637 fQAlist->AddAt(histo1D, 3 + icond * fgkNQAhistos);
809a4336 638 histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
639 histo1D->GetYaxis()->SetTitle("Number of Tracks");
8c1c76e9 640 qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClusters%s",GetName(),cutstr[icond].Data()), "Number of TPC clusters", 161, 0, 160)), 4 + icond * fgkNQAhistos);
641 fQAlist->AddAt(histo1D, 4 + icond * fgkNQAhistos);
3a72645a 642 histo1D->GetXaxis()->SetTitle("Number of TPC clusters");
643 histo1D->GetYaxis()->SetTitle("Number of Tracks");
8c1c76e9 644 qaList->AddAt((histo1D = new TH1F(Form("%s_itsPixel%s",GetName(),cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 5 + icond * fgkNQAhistos);
645 fQAlist->AddAt(histo1D, 5 + icond * fgkNQAhistos);
809a4336 646 histo1D->GetXaxis()->SetTitle("ITS Pixel");
647 histo1D->GetYaxis()->SetTitle("Number of Tracks");
648 Int_t first = histo1D->GetXaxis()->GetFirst();
649 TString binNames[6] = { "First", "Second", "Both", "None", "Exclusive First", "Exclusive Second"};
650 for(Int_t ilabel = 0; ilabel < 6; ilabel++)
651 histo1D->GetXaxis()->SetBinLabel(first + ilabel, binNames[ilabel].Data());
8c1c76e9 652 qaList->AddAt((histo1D = new TH1F(Form("%s_trackStatus%s",GetName(),cutstr[icond].Data()), "Track Status", 5, 0, 5)), 6 + icond * fgkNQAhistos);
653 fQAlist->AddAt(histo1D, 6 + icond * fgkNQAhistos);
654 histo1D->GetXaxis()->SetTitle("Track Status Bit");
655 histo1D->GetYaxis()->SetTitle("Number of tracks");
656 TString tsnames[5] = {"All", "TOFPID", "No TOFmismatch", "EMCALmatch","No TPC shared bit"};
657 for(Int_t istat = 0; istat < 5; istat++) histo1D->GetXaxis()->SetBinLabel(istat + 1, tsnames[istat].Data());
658 qaList->AddAt((histo1D = new TH1F(Form("%s_tpcdEdxClusters%s",GetName(),cutstr[icond].Data()), "Number of TPC clusters for dEdx calculation", 161, 0, 160)), 7 + icond * fgkNQAhistos);
659 fQAlist->AddAt(histo1D, 7 + icond * fgkNQAhistos);
660 histo1D->GetXaxis()->SetTitle("Number of TPC clusters for dEdx calculation");
661 histo1D->GetYaxis()->SetTitle("counts");
cedf0381 662 qaList->AddAt((histo1D = new TH1F(Form("%s_trdchi2perTracklet%s",GetName(),cutstr[icond].Data()), "chi2 per TRD tracklet", 100, 0, 10)), 8 + icond * fgkNQAhistos);
663 fQAlist->AddAt(histo1D, 8 + icond * fgkNQAhistos);
664 histo1D->GetXaxis()->SetTitle("Chi2 per TRD Tracklet");
665 histo1D->GetYaxis()->SetTitle("Number of Tracks");
76d0b522 666 qaList->AddAt((histo1D = new TH1F(Form("%s_ITScls%s",GetName(),cutstr[icond].Data()), "ITScls", 6, 0, 6)), 9 + icond * fgkNQAhistos);
667 fQAlist->AddAt(histo1D, 9 + icond * fgkNQAhistos);
668 histo1D->GetXaxis()->SetTitle("ITScls");
669 histo1D->GetYaxis()->SetTitle("Number of ITS cls");
cedf0381 670
809a4336 671 }
809a4336 672 // Add cut correlation
8c1c76e9 673 qaList->AddAt((histo2D = new TH2F(Form("%s_cutcorrelation",GetName()), "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2 * fgkNQAhistos);
674 fQAlist->AddAt(histo2D, 2 * fgkNQAhistos);
675 TString labels[kNcuts] = {"MinImpactParamR", "MaxImpactParamR", "MinImpactParamZ", "MaxImpactParamZ", "ClusterRatioTPC", "MinTrackletsTRD", "ITSpixel", "kMinHFEImpactParamR", "kMinHFEImpactParamNsigmaR", "TPC Number of clusters" "Fraction Shared TPC clusters", "TOFPID", "No TOFmismatch", "EMCALmatch", "ImpactParam"};
809a4336 676 Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
677 for(Int_t icut = 0; icut < kNcuts; icut++){
678 histo2D->GetXaxis()->SetBinLabel(firstx + icut, labels[icut].Data());
679 histo2D->GetYaxis()->SetBinLabel(firsty + icut, labels[icut].Data());
680 }
809a4336 681}
dbe3abbe 682
683//______________________________________________________
684void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
685 for(Int_t ibit = 32; ibit--; )
686 printf("%d", bitmap & BIT(ibit) ? 1 : 0);
687 printf("\n");
688}
689
690//______________________________________________________
75d81601 691Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
dbe3abbe 692 //
693 // Check whether ITS area is dead
694 //
695 Bool_t status;
696 switch(itsStatus){
697 case 2: status = kFALSE; break;
698 case 3: status = kFALSE; break;
699 case 7: status = kFALSE; break;
700 default: status = kTRUE;
701 }
702 return status;
703}
3a72645a 704
3a72645a 705//______________________________________________________
cedf0381 706Int_t AliHFEextraCuts::GetITSstatus(const AliVTrack * const track, Int_t layer) const {
3a72645a 707 //
708 // Check ITS layer status
709 //
710 Int_t status = 0;
711 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
712 Int_t det;
713 Float_t xloc, zloc;
cedf0381 714 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
bf892a6a 715 if(esdtrack) esdtrack->GetITSModuleIndexInfo(layer, det, status, xloc, zloc);
3a72645a 716 }
717 return status;
718}
719
8c1c76e9 720
721//______________________________________________________
722Bool_t AliHFEextraCuts::GetTPCCountSharedMapBitsAboveThreshold(AliVTrack *track){
723 //
724 // Checks if number of shared bits is above 1
725 //
726 Int_t nsharebit = 1;
11ff28c5 727 const TBits *shared;
728 if(track->IsA() == AliESDtrack::Class())
729 shared = &((static_cast<AliESDtrack *>(track))->GetTPCSharedMap());
730 else if(track->IsA() == AliAODTrack::Class())
731 shared = &((static_cast<AliAODTrack *>(track))->GetTPCSharedMap());
732 else
733 return 0;
734
735 if(shared->CountBits() >= 2) nsharebit=1;
736 else nsharebit=0;
737
8c1c76e9 738 return nsharebit;
739}
740
3a72645a 741//______________________________________________________
e3ae862b 742UInt_t AliHFEextraCuts::GetTPCncls(AliVTrack *track){
76d0b522 743 //
744 // Get Number of findable clusters in the TPC
745 //
746 Int_t nClusters = 0; // in case no Information available consider all clusters findable
11ff28c5 747 TClass *type = track->IsA();
748 if(TESTBIT(fTPCclusterDef, kFoundIter1)){
76d0b522 749 if(type == AliESDtrack::Class()){
11ff28c5 750 AliDebug(2, ("Using def kFoundIter1"));
751 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
752 nClusters = esdtrack->GetTPCNclsIter1();
753 } else {
754 AliDebug(2, ("Number of clusters in iteration 1 not available for AOD tracks"));
bf892a6a 755 }
11ff28c5 756 } else if(TESTBIT(fTPCclusterDef, kCrossedRows)){
757 AliDebug(2, ("Using def kCrossedRows"));
758 nClusters = static_cast<UInt_t>(track->GetTPCClusterInfo(2,1));
76d0b522 759 } else if(TESTBIT(fTPCclusterDef, kFound)){
11ff28c5 760 AliDebug(2, ("Using def kFound"));
76d0b522 761 nClusters = track->GetTPCNcls();
762 }
763 else {
764 AliDebug(2, ("Using def kFoundAll"));
765 if(type == AliESDtrack::Class()){
766 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
767 const TBits &clusterTPC = esdtrack->GetTPCClusterMap();
768 nClusters = clusterTPC.CountBits();
769 }
770 else {
771 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
772 const TBits &clusterTPC = aodtrack->GetTPCClusterMap();
773 nClusters = clusterTPC.CountBits();
774 }
3a72645a 775 }
76d0b522 776 return nClusters;
3a72645a 777}
778
9250ffbf 779//______________________________________________________
780Float_t AliHFEextraCuts::GetTPCsharedClustersRatio(AliVTrack *track){
781 //
782 // Get fraction of shared TPC clusters
783 //
784 Float_t fracClustersTPCShared = 0.0;
11ff28c5 785 Int_t nClustersTPC = track->GetTPCNcls();
786 Int_t nClustersTPCShared = 0;
787 TClass *type = track->IsA();
788 if(type == AliESDtrack::Class()){
789 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
790 nClustersTPCShared = esdtrack->GetTPCnclsS();
76d0b522 791 } else if(type == AliAODTrack::Class()){
11ff28c5 792 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
793 const TBits &shared = aodtrack->GetTPCSharedMap();
794 nClustersTPCShared = shared.CountBits(0) - shared.CountBits(159);
795 }
796 if (nClustersTPC!=0) {
76d0b522 797 fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
9250ffbf 798 }
799 return fracClustersTPCShared;
800}
801
e3ae862b 802//______________________________________________________
803Double_t AliHFEextraCuts::GetTPCclusterRatio(AliVTrack *track){
804 //
805 // Get Ratio of found / findable clusters for different definitions
806 // Only implemented for ESD tracks
807 //
808 Double_t clusterRatio = 1.; // in case no Information available consider all clusters findable
11ff28c5 809 TClass *type = track->IsA();
810 if(TESTBIT(fTPCclusterRatioDef, kCROverFindable)){
811 AliDebug(2, "Using ratio def kCROverFindable");
812 clusterRatio = track->GetTPCNclsF() ? track->GetTPCClusterInfo(2,1)/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // crossed rows/findable
813 } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverCR)){
814 AliDebug(2, "Using ratio def kFoundOverCR");
815 clusterRatio = track->GetTPCClusterInfo(2,0); // found/crossed rows
816 } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindableIter1)){
817 if(type == AliESDtrack::Class()){
818 AliDebug(2, "Using ratio def kFoundOverFindableIter1");
819 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
820 clusterRatio = esdtrack->GetTPCNclsFIter1() ? static_cast<Double_t>(esdtrack->GetTPCNclsIter1())/static_cast<Double_t>(esdtrack->GetTPCNclsFIter1()) : 1.; // crossed
821 } else {
822 AliDebug(2, "Cluster ratio after iteration 1 not possible for AOD tracks");
823 clusterRatio = 1.;
e3ae862b 824 }
76d0b522 825 } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindable)) {
11ff28c5 826 AliDebug(2, "Using ratio def kFoundOverFindable");
827 clusterRatio = track->GetTPCNclsF() ? static_cast<Double_t>(track->GetTPCNcls())/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // found/findable
e3ae862b 828 }
76d0b522 829 else {
830 AliDebug(2, "Using ratio def kFoundAllOverFindable");
831 Int_t allclusters = 0;
832 if(type == AliESDtrack::Class()){
833 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
834 const TBits &clusterTPC = esdtrack->GetTPCClusterMap();
835 allclusters = clusterTPC.CountBits();
836 }
837 else {
838 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
839 const TBits &clusterTPC = aodtrack->GetTPCClusterMap();
840 allclusters = clusterTPC.CountBits();
841 }
842 clusterRatio = track->GetTPCNclsF() ? static_cast<Double_t>(allclusters)/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // foundall/findable
843 }
e3ae862b 844 return clusterRatio;
e3ae862b 845
76d0b522 846}
847//___________________________________________________
3a72645a 848void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z){
11ff28c5 849 //
850 // Get impact parameter
851 //
959ea9d8 852
853 const Double_t kBeampiperadius=3.;
854 Double_t dcaD[2]={-999.,-999.},
855 covD[3]={-999.,-999.,-999.};
11ff28c5 856 TClass *type = track->IsA();
857 if(type == AliESDtrack::Class()){
858 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
859 esdtrack->GetImpactParameters(radial, z);
860 }
861 else if(type == AliAODTrack::Class()){
959ea9d8 862
863 //case of AOD tracks
864 AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
865 if(!aodevent) {
866 AliDebug(1, "No aod event available\n");
867 return;
bf892a6a 868 }
959ea9d8 869
870 //Case ESD track: take copy constructor
871 AliAODTrack *aodtrack = NULL;
872 AliAODTrack *tmptrack = dynamic_cast<AliAODTrack *>(track);
873 if(tmptrack) aodtrack = new AliAODTrack(*tmptrack);
874
875 AliAODVertex *vtxAODSkip = aodevent->GetPrimaryVertex();
876 if(!vtxAODSkip) return;
877 AliExternalTrackParam etp; etp.CopyFromVTrack(aodtrack);
63bdf450 878 if(etp.PropagateToDCA(vtxAODSkip, aodevent->GetMagneticField(), kBeampiperadius, dcaD, covD)) {
959ea9d8 879 radial = dcaD[0];
880 z = dcaD[1];
881 }
882 //if(vtxAODSkip) delete vtxAODSkip;
883 if(aodtrack) delete aodtrack;
11ff28c5 884 }
885}
886//______________________________________________________
887Bool_t AliHFEextraCuts::IsKinkDaughter(AliVTrack *track){
888 //
889 // Is kink Daughter
890 //
891 TClass *type = track->IsA();
892 if(type == AliESDtrack::Class()){
893 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
894 if(!esdtrack) return kFALSE;
76d0b522 895 if(esdtrack->GetKinkIndex(0)>0) return kTRUE;
896 else return kFALSE;
11ff28c5 897
898 }
899 else if(type == AliAODTrack::Class()){
900 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
901 if(aodtrack){
902 AliAODVertex *aodvertex = aodtrack->GetProdVertex();
903 if(!aodvertex) return kFALSE;
904 if(aodvertex->GetType()==AliAODVertex::kKink) return kTRUE;
905 else return kFALSE;
906 }
907 else return kFALSE;
908 }
909
910 return kFALSE;
911}
912//______________________________________________________
913Bool_t AliHFEextraCuts::IsKinkMother(AliVTrack *track){
914 //
bbeea05b 915 // Is kink Mother
0e30407a 916 //
917
11ff28c5 918 TClass *type = track->IsA();
919 if(type == AliESDtrack::Class()){
920 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
921 if(!esdtrack) return kFALSE;
922 if(esdtrack->GetKinkIndex(0)!=0) return kTRUE;
923 else return kFALSE;
bbeea05b 924 } else if(type == AliAODTrack::Class()){
925 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
926 for(int ikink = 0; ikink < fNumberKinkMothers; ikink++){
927 if(fListKinkMothers[ikink] == aodtrack->GetID()) return kTRUE;
928 }
11ff28c5 929 }
11ff28c5 930
931 return kFALSE;
932
3a72645a 933}
cedf0381 934
935//______________________________________________________
936Float_t AliHFEextraCuts::GetTRDchi(AliVTrack *track){
937 //
938 // Get TRDchi2
939 //
940 Int_t ntls(0);
941 TClass *type = track->IsA();
942 if(type == AliESDtrack::Class()){
943 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
944 ntls = esdtrack->GetTRDntracklets();
945 return ntls ? esdtrack->GetTRDchi2()/ntls : -999;
946 }
947 else if(type == AliAODTrack::Class()){
948 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
949 if(aodtrack){
950 return 999.;
951 }
952 }
953
954 return 999.;
955
956}
957
11ff28c5 958//______________________________________________________
959Int_t AliHFEextraCuts::GetITSNbOfcls(AliVTrack *track){
960 //
961 // Get ITS nb of clusters
962 //
963 TClass *type = track->IsA();
964 if(type == AliESDtrack::Class()){
965 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
966 return esdtrack->GetITSclusters(0);
3a72645a 967
11ff28c5 968 }
969 else if(type == AliAODTrack::Class()){
970 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
971 if(aodtrack){
972 return aodtrack->GetITSNcls();
973 }
974 }
975
976 return 0;
977}
3a72645a 978//______________________________________________________
63bdf450 979void AliHFEextraCuts::GetHFEImpactParameters(const AliVTrack * const track, Double_t &dcaxy, Double_t &dcansigmaxy){
959ea9d8 980 //
981 // Get HFE impact parameter (with recalculated primary vertex)
982 //
983 dcaxy=0;
984 dcansigmaxy=0;
3a72645a 985 if(!fEvent){
986 AliDebug(1, "No Input event available\n");
987 return;
988 }
3a72645a 989 TString type = track->IsA()->GetName();
959ea9d8 990 const Double_t kBeampiperadius=3.;
991 Double_t dcaD[2]={-999.,-999.},
992 covD[3]={-999.,-999.,-999.};
63bdf450 993 Bool_t isRecalcVertex(kFALSE);
959ea9d8 994
995 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
996 //case of ESD tracks
997 AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
998 if(!esdevent) {
999 AliDebug(1, "No esd event available\n");
1000 return;
1001 }
1002
959ea9d8 1003 const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
1004 if(!vtxESDSkip) return;
959ea9d8 1005
63bdf450 1006 //case ESD track: take copy constructor
1007 const AliESDtrack *tmptrack = dynamic_cast<const AliESDtrack *>(track);
1008 if(tmptrack){
1009
1010 if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
1011 vtxESDSkip = RemoveDaughtersFromPrimaryVtx(esdevent, track);
1012 isRecalcVertex = kTRUE;
1013 }
1014
1015 if(vtxESDSkip){
1016 AliESDtrack esdtrack(*tmptrack);
1017 fMagField = fEvent->GetMagneticField();
1018 if(esdtrack.PropagateToDCA(vtxESDSkip, fMagField, kBeampiperadius, dcaD, covD)){
1019 dcaxy = dcaD[0];
1020 if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
1021 }
1022 if(isRecalcVertex) delete vtxESDSkip;
1023 }
959ea9d8 1024 }
959ea9d8 1025 }
1026 else {
1027 //case of AOD tracks
1028 AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
1029 if(!aodevent) {
1030 AliDebug(1, "No aod event available\n");
1031 return;
1032 }
1033
959ea9d8 1034 AliAODVertex *vtxAODSkip = aodevent->GetPrimaryVertex();
1035 if(!vtxAODSkip) return;
63bdf450 1036
1037 //Case ESD track: take copy constructor
1038 const AliAODTrack *tmptrack = dynamic_cast<const AliAODTrack *>(track);
1039 if(tmptrack){
1040
1041 if(vtxAODSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
1042 vtxAODSkip = RemoveDaughtersFromPrimaryVtx(aodevent, track);
1043 isRecalcVertex = kTRUE;
1044 }
1045 if(vtxAODSkip){
1046 AliAODTrack aodtrack(*tmptrack);
1047 AliExternalTrackParam etp; etp.CopyFromVTrack(&aodtrack);
1048 fMagField = aodevent->GetMagneticField();
1049 if(etp.PropagateToDCA(vtxAODSkip,fMagField, kBeampiperadius, dcaD, covD)) {
1050 dcaxy = dcaD[0];
1051 if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
1052 }
1053 if(isRecalcVertex) delete vtxAODSkip;
1054 }
959ea9d8 1055 }
3a72645a 1056 }
3a72645a 1057}
1058
11ff28c5 1059//______________________________________________________
63bdf450 1060void AliHFEextraCuts::GetHFEImpactParameters(const AliVTrack * const track, Double_t dcaD[2], Double_t covD[3]){
11ff28c5 1061 //
1062 // Get HFE impact parameter (with recalculated primary vertex)
1063 //
1064 if(!fEvent){
1065 AliDebug(1, "No Input event available\n");
1066 return;
1067 }
1068 const Double_t kBeampiperadius=3.;
1069 TString type = track->IsA()->GetName();
63bdf450 1070 Bool_t isRecalcVertex(kFALSE);
11ff28c5 1071
959ea9d8 1072 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
1073 //case of ESD tracks
1074 AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
1075 if(!esdevent) {
1076 AliDebug(1, "No esd event available\n");
1077 return;
1078 }
1079
63bdf450 1080 // Check whether primary vertex is available
959ea9d8 1081 const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
1082 if(!vtxESDSkip) return;
959ea9d8 1083
63bdf450 1084 const AliESDtrack *tmptrack = dynamic_cast<const AliESDtrack *>(track);
1085 if(tmptrack){
1086 //case ESD track: take copy constructor
1087
1088 if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
1089 vtxESDSkip = RemoveDaughtersFromPrimaryVtx(esdevent, track);
1090 isRecalcVertex = kTRUE;
1091 }
1092 if(vtxESDSkip){
1093 AliESDtrack esdtrack(*tmptrack);
1094 fMagField = fEvent->GetMagneticField();
1095 esdtrack.PropagateToDCA(vtxESDSkip, fMagField, kBeampiperadius, dcaD, covD);
959ea9d8 1096
63bdf450 1097 if(isRecalcVertex) delete vtxESDSkip;
1098 }
1099 }
959ea9d8 1100 }
1101 else {
1102 //case of AOD tracks
1103 AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
1104 if(!aodevent) {
1105 AliDebug(1, "No aod event available\n");
1106 return;
1107 }
1108
63bdf450 1109 // Check whether primary vertex is available
959ea9d8 1110 AliAODVertex *vtxAODSkip = aodevent->GetPrimaryVertex();
1111 if(!vtxAODSkip) return;
959ea9d8 1112
63bdf450 1113 //Case ESD track: take copy constructor
1114 const AliAODTrack *tmptrack = dynamic_cast<const AliAODTrack *>(track);
1115 if(tmptrack){
1116
1117 if(vtxAODSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
1118 vtxAODSkip = RemoveDaughtersFromPrimaryVtx(aodevent, track);
1119 isRecalcVertex = kTRUE;
1120 }
1121 if(vtxAODSkip){
1122 AliAODTrack aodtrack(*tmptrack);
1123 AliExternalTrackParam etp; etp.CopyFromVTrack(&aodtrack);
1124 fMagField = aodevent->GetMagneticField();
1125 etp.PropagateToDCA(vtxAODSkip, fMagField, kBeampiperadius, dcaD, covD);
1126
1127 if(isRecalcVertex) delete vtxAODSkip;
1128 }
1129 }
11ff28c5 1130 }
1131}
3a72645a 1132
1133//______________________________________________________
63bdf450 1134void AliHFEextraCuts::GetHFEImpactParameterCuts(const AliVTrack * const track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
3a72645a 1135 //
1136 // Get HFE impact parameter cut(pt dependent)
1137 //
1138
63bdf450 1139 Double_t pt = track->Pt();
1140 hfeimpactRcut = fIPcutParam[0]+fIPcutParam[1]*exp(fIPcutParam[2]*pt); // abs R cut
1141 hfeimpactnsigmaRcut = fIPcutParam[3]; // sigma cut
3a72645a 1142}
9250ffbf 1143//______________________________________________________
63bdf450 1144void AliHFEextraCuts::GetMaxImpactParameterCutR(const AliVTrack * const track, Double_t &maximpactRcut){
9250ffbf 1145 //
1146 // Get max impact parameter cut r (pt dependent)
1147 //
1148
63bdf450 1149 Double_t pt = track->Pt();
1150 if(pt > 0.15) {
1151 maximpactRcut = 0.0182 + 0.035/TMath::Power(pt,1.01); // abs R cut
9250ffbf 1152 }
63bdf450 1153 else maximpactRcut = 9999999999.0;
9250ffbf 1154}
1c051dd4 1155//______________________________________________________
63bdf450 1156void AliHFEextraCuts::GetTOFsignalDxDz(const AliVTrack * const track, Double_t &tofsignalDx, Double_t &tofsignalDz){
1c051dd4 1157 //
1158 // TOF matching
1159 //
1160
1161 TString type = track->IsA()->GetName();
1162 if(!type.CompareTo("AliESDtrack")){
63bdf450 1163 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
1c051dd4 1164 if(!esdtrack) return;
1165 tofsignalDx = esdtrack->GetTOFsignalDx();
1166 tofsignalDz = esdtrack->GetTOFsignalDz();
1167 }
1168
1169}
cedf0381 1170
2e4ed823 1171//______________________________________________________
1172Bool_t AliHFEextraCuts::MatchTOFlabel(const AliVTrack *const track) const {
1173 //
1174 // Check whether the TOF label is the same as the track label
1175 //
1176 const AliESDtrack *esdtrk(NULL);
1177 const AliAODTrack *aodtrk(NULL);
1178 int trklabel(99999), toflabel[3] = {99999,99999,99999};
1179 if((esdtrk = dynamic_cast<const AliESDtrack *>(track))){
1180 trklabel = esdtrk->GetLabel();
1181 esdtrk->GetTOFLabel(toflabel);
1182 } else if((aodtrk = dynamic_cast<const AliAODTrack *>(track))){
2b66effc 1183 trklabel = aodtrk->GetLabel();
1184 aodtrk->GetTOFLabel(toflabel);
2e4ed823 1185 } else return kFALSE;
1186 if(TMath::Abs(trklabel) == TMath::Abs(toflabel[0])) return kTRUE;
1187 return kFALSE;
1188}
cedf0381 1189//______________________________________________________
1190Bool_t AliHFEextraCuts::CheckITSpattern(const AliVTrack *const track) const {
1191 //
1192 // Check if every ITS layer, which has a module which is alive, also
1193 // has an ITS cluster
1194 //
1195 Bool_t patternOK(kTRUE);
1196 Int_t status(0);
1197 for(Int_t ily = 0; ily < 6; ily++){
1198 status = GetITSstatus(track, ily);
1199 if(CheckITSstatus(status)){
1200 // pixel alive, check whether layer has a cluster
1201 if(!TESTBIT(track->GetITSClusterMap(),ily)){
1202 // No cluster even though pixel is alive - reject track
1203 patternOK = kFALSE;
1204 break;
1205 }
1206 }
1207 }
1208 return patternOK;
1209}
959ea9d8 1210
1211//---------------------------------------------------------------------------
63bdf450 1212const AliVVertex* AliHFEextraCuts::RemoveDaughtersFromPrimaryVtx(const AliESDEvent * const esdevent, const AliVTrack * const track) {
959ea9d8 1213 //
1214 // This method returns a primary vertex without the daughter tracks of the
1215 // candidate and it recalculates the impact parameters and errors for ESD tracks.
1216 //
1217 // The output vertex is created with "new".
1218 //
1219
91e50e2b 1220 //const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
959ea9d8 1221
1222 AliVertexerTracks vertexer(fEvent->GetMagneticField());
1223 vertexer.SetITSMode();
1224 vertexer.SetMinClusters(4);
1225 Int_t skipped[2];
1226 skipped[0] = track->GetID();
1227 vertexer.SetSkipTracks(1,skipped);
1228
1229 //diamond constraint
1230 vertexer.SetConstraintOn();
1231 Float_t diamondcovxy[3];
1232 esdevent->GetDiamondCovXY(diamondcovxy);
1233 Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
1234 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
63bdf450 1235 AliESDVertex diamond(pos,cov,1.,1);
1236 vertexer.SetVtxStart(&diamond);
959ea9d8 1237
63bdf450 1238 const AliVVertex *vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
959ea9d8 1239
1240 return vtxESDSkip;
1241}
1242
1243//---------------------------------------------------------------------------
63bdf450 1244AliAODVertex* AliHFEextraCuts::RemoveDaughtersFromPrimaryVtx(const AliAODEvent * const aod, const AliVTrack * const track) {
959ea9d8 1245 //
1246 // This method returns a primary vertex without the daughter tracks of the
1247 // candidate and it recalculates the impact parameters and errors for AOD tracks.
1248 // The output vertex is created with "new".
1249 //
1250
1251 AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
1252 if(!vtxAOD) return 0;
1253 TString title=vtxAOD->GetTitle();
1254 if(!title.Contains("VertexerTracks")) return 0;
1255
63bdf450 1256 AliVertexerTracks vertexer(aod->GetMagneticField());
959ea9d8 1257
63bdf450 1258 vertexer.SetITSMode();
1259 vertexer.SetMinClusters(3);
1260 vertexer.SetConstraintOff();
959ea9d8 1261
1262 if(title.Contains("WithConstraint")) {
1263 Float_t diamondcovxy[3];
1264 aod->GetDiamondCovXY(diamondcovxy);
1265 Double_t pos[3]={aod->GetDiamondX(),aod->GetDiamondY(),0.};
1266 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
63bdf450 1267 AliESDVertex diamond(pos,cov,1.,1);
1268 vertexer.SetVtxStart(&diamond);
959ea9d8 1269 }
1270 Int_t skipped[2]; for(Int_t i=0;i<2;i++) skipped[i]=-1;
1271 Int_t id = (Int_t)track->GetID();
1272 if(!(id<0)) skipped[0] = id;
1273
1274 /*// exclude tracks with global constrained parameters
1275 Int_t nTracks=aod->GetNumberOfTracks();
1276 for(Int_t i=0; i<nTracks; i++){
1277 t = aod->GetTrack(i);
1278 if(t->TestFilterMask(512)){
1279 id = (Int_t)t->GetID();
1280 skipped[nTrksToSkip++] = id;
1281 }
1282 }*/
1283
63bdf450 1284 vertexer.SetSkipTracks(1,skipped);
1285 AliESDVertex *vtxESDNew = vertexer.FindPrimaryVertex(aod);
959ea9d8 1286
1287 if(!vtxESDNew) return 0;
1288 if(vtxESDNew->GetNContributors()<=0) {
1289 delete vtxESDNew; vtxESDNew=NULL;
1290 return 0;
1291 }
1292
1293 // convert to AliAODVertex
1294 Double_t pos[3],cov[6],chi2perNDF;
1295 vtxESDNew->GetXYZ(pos); // position
1296 vtxESDNew->GetCovMatrix(cov); //covariance matrix
1297 chi2perNDF = vtxESDNew->GetChi2toNDF();
1298 delete vtxESDNew; vtxESDNew=NULL;
1299
1300 AliAODVertex *vtxAODNew = new AliAODVertex(pos,cov,chi2perNDF);
1301
1302 return vtxAODNew;
1303}