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