]> 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
Removed hardcoded AliESDs.root in AliESDInputHandler(RP)::Notify()
[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
3a72645a 33#include "AliAODTrack.h"
34#include "AliAODPid.h"
11ff28c5 35#include "AliAODVertex.h"
3a72645a 36#include "AliESDEvent.h"
37#include "AliESDVertex.h"
809a4336 38#include "AliESDtrack.h"
70da6c5a 39#include "AliLog.h"
809a4336 40#include "AliMCParticle.h"
3a72645a 41#include "AliVEvent.h"
42#include "AliVTrack.h"
43#include "AliVParticle.h"
44#include "AliVertexerTracks.h"
45#include "AliVVertex.h"
809a4336 46
47#include "AliHFEextraCuts.h"
48
49ClassImp(AliHFEextraCuts)
50
e17c1f86 51const Int_t AliHFEextraCuts::fgkNQAhistos = 8;
8c1c76e9 52
809a4336 53//______________________________________________________
54AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
55 AliCFCutBase(name, title),
3a72645a 56 fEvent(NULL),
809a4336 57 fCutCorrelation(0),
58 fRequirements(0),
3a72645a 59 fMinNClustersTPC(0),
e17c1f86 60 fMinNClustersTPCPID(0),
809a4336 61 fClusterRatioTPC(0.),
62 fMinTrackletsTRD(0),
11ff28c5 63 fMinNbITScls(0),
aa1417ee 64 fTRDtrackletsExact(0),
809a4336 65 fPixelITS(0),
e3ae862b 66 fTPCclusterDef(0),
67 fTPCclusterRatioDef(0),
8c1c76e9 68 fFractionTPCShared(-1.0),
11ff28c5 69 fAbsHFEImpactParamNsigmaR(kTRUE),
70da6c5a 70 fCheck(kFALSE),
71 fQAlist(0x0) ,
dbe3abbe 72 fDebugLevel(0)
809a4336 73{
74 //
75 // Default Constructor
76 //
11ff28c5 77 //printf("Set the number of min ITS clusters %d\n",(Int_t)fMinNbITScls);
809a4336 78 memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
c2690925 79 memset(fIPcutParam, 0, sizeof(Float_t) * 4);
809a4336 80}
81
82//______________________________________________________
83AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
84 AliCFCutBase(c),
3a72645a 85 fEvent(c.fEvent),
809a4336 86 fCutCorrelation(c.fCutCorrelation),
87 fRequirements(c.fRequirements),
3a72645a 88 fMinNClustersTPC(c.fMinNClustersTPC),
e17c1f86 89 fMinNClustersTPCPID(c.fMinNClustersTPCPID),
809a4336 90 fClusterRatioTPC(c.fClusterRatioTPC),
91 fMinTrackletsTRD(c.fMinTrackletsTRD),
11ff28c5 92 fMinNbITScls(c.fMinNbITScls),
aa1417ee 93 fTRDtrackletsExact(c.fTRDtrackletsExact),
809a4336 94 fPixelITS(c.fPixelITS),
e3ae862b 95 fTPCclusterDef(c.fTPCclusterDef),
96 fTPCclusterRatioDef(c.fTPCclusterRatioDef),
8c1c76e9 97 fFractionTPCShared(c.fFractionTPCShared),
11ff28c5 98 fAbsHFEImpactParamNsigmaR(c.fAbsHFEImpactParamNsigmaR),
75d81601 99 fCheck(c.fCheck),
dbe3abbe 100 fQAlist(0x0),
70da6c5a 101 fDebugLevel(0)
809a4336 102{
103 //
104 // Copy constructor
105 // Performs a deep copy
106 //
107 memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
c2690925 108 memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
809a4336 109 if(IsQAOn()){
110 fIsQAOn = kTRUE;
111 fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
bf892a6a 112 if(fQAlist) fQAlist->SetOwner();
809a4336 113 }
114}
115
116//______________________________________________________
117AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
118 //
119 // Assignment operator
120 //
121 if(this != &c){
122 AliCFCutBase::operator=(c);
3a72645a 123 fEvent = c.fEvent;
809a4336 124 fCutCorrelation = c.fCutCorrelation;
125 fRequirements = c.fRequirements;
126 fClusterRatioTPC = c.fClusterRatioTPC;
3a72645a 127 fMinNClustersTPC = c.fMinNClustersTPC;
e17c1f86 128 fMinNClustersTPCPID = c.fMinNClustersTPCPID;
809a4336 129 fMinTrackletsTRD = c.fMinTrackletsTRD;
11ff28c5 130 fMinNbITScls = c.fMinNbITScls;
aa1417ee 131 fTRDtrackletsExact = c.fTRDtrackletsExact;
809a4336 132 fPixelITS = c.fPixelITS;
e3ae862b 133 fTPCclusterDef = c.fTPCclusterDef;
134 fTPCclusterRatioDef = c.fTPCclusterRatioDef;
8c1c76e9 135 fFractionTPCShared = c.fFractionTPCShared;
11ff28c5 136 fAbsHFEImpactParamNsigmaR = c.fAbsHFEImpactParamNsigmaR;
75d81601 137 fCheck = c.fCheck;
dbe3abbe 138 fDebugLevel = c.fDebugLevel;
809a4336 139 memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
c2690925 140 memcpy(fIPcutParam, c.fIPcutParam, sizeof(Float_t) * 4);
809a4336 141 if(IsQAOn()){
142 fIsQAOn = kTRUE;
143 fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
bf892a6a 144 if(fQAlist) fQAlist->SetOwner();
809a4336 145 }else fQAlist = 0x0;
146 }
147 return *this;
148}
149
150//______________________________________________________
151AliHFEextraCuts::~AliHFEextraCuts(){
152 //
153 // Destructor
154 //
3a72645a 155 if(fQAlist) delete fQAlist;
156}
157
158//______________________________________________________
159void AliHFEextraCuts::SetRecEventInfo(const TObject *event){
160 //
161 // Set Virtual event an make a copy
162 //
163 if (!event) {
164 AliError("Pointer to AliVEvent !");
165 return;
166 }
167 TString className(event->ClassName());
168 if (! (className.CompareTo("AliESDEvent")==0 || className.CompareTo("AliAODEvent")==0)) {
169 AliError("argument must point to an AliESDEvent or AliAODEvent !");
170 return ;
809a4336 171 }
3a72645a 172 fEvent = (AliVEvent*) event;
173
809a4336 174}
175
176//______________________________________________________
9eeae5d5 177Bool_t AliHFEextraCuts::IsSelected(TObject *o){
809a4336 178 //
179 // Steering function for the track selection
180 //
11ff28c5 181 TClass *type = o->IsA();
182 AliDebug(2, Form("Object type %s", type->GetName()));
183 if((type == AliESDtrack::Class()) || (type == AliAODTrack::Class())){
3a72645a 184 return CheckRecCuts(dynamic_cast<AliVTrack *>(o));
809a4336 185 }
3a72645a 186 return CheckMCCuts(dynamic_cast<AliVParticle *>(o));
809a4336 187}
188
189//______________________________________________________
3a72645a 190Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
809a4336 191 //
192 // Checks cuts on reconstructed tracks
193 // returns true if track is selected
194 // QA histograms are filled before track selection and for
195 // selected tracks after track selection
196 //
70da6c5a 197 AliDebug(1, "Called");
809a4336 198 ULong64_t survivedCut = 0; // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
3a72645a 199 if(IsQAOn()) FillQAhistosRec(track, kBeforeCuts);
809a4336 200 // Apply cuts
e3ae862b 201 Float_t impactR, impactZ;
3a72645a 202 Double_t hfeimpactR, hfeimpactnsigmaR;
203 Double_t hfeimpactRcut, hfeimpactnsigmaRcut;
9250ffbf 204 Double_t maximpactRcut;
3a72645a 205 GetImpactParameters(track, impactR, impactZ);
206 if(TESTBIT(fRequirements, kMinHFEImpactParamR) || TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
207 // Protection for PbPb
208 GetHFEImpactParameterCuts(track, hfeimpactRcut, hfeimpactnsigmaRcut);
209 GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
210 }
11ff28c5 211 Int_t nclsITS = GetITSNbOfcls(track);
e3ae862b 212 UInt_t nclsTPC = GetTPCncls(track);
11ff28c5 213 UInt_t nclsTPCPID = track->GetTPCsignalN();
809a4336 214 // printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
9250ffbf 215 Float_t fractionSharedClustersTPC = GetTPCsharedClustersRatio(track);
e3ae862b 216 Double_t ratioTPC = GetTPCclusterRatio(track);
809a4336 217 UChar_t trdTracklets;
8c1c76e9 218 trdTracklets = track->GetTRDntrackletsPID();
809a4336 219 UChar_t itsPixel = track->GetITSClusterMap();
3a72645a 220 Int_t status1 = GetITSstatus(track, 0);
221 Int_t status2 = GetITSstatus(track, 1);
70da6c5a 222 Bool_t statusL0 = CheckITSstatus(status1);
223 Bool_t statusL1 = CheckITSstatus(status2);
8c1c76e9 224 if(TESTBIT(fRequirements, kTPCfractionShared)) {
e17c1f86 225 // cut on max fraction of shared TPC clusters
226 if(TMath::Abs(fractionSharedClustersTPC) < fFractionTPCShared) SETBIT(survivedCut, kTPCfractionShared);
9250ffbf 227 }
809a4336 228 if(TESTBIT(fRequirements, kMinImpactParamR)){
229 // cut on min. Impact Parameter in Radial direction
75d81601 230 if(TMath::Abs(impactR) >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
809a4336 231 }
232 if(TESTBIT(fRequirements, kMinImpactParamZ)){
233 // cut on min. Impact Parameter in Z direction
75d81601 234 if(TMath::Abs(impactZ) >= fImpactParamCut[1]) SETBIT(survivedCut, kMinImpactParamZ);
809a4336 235 }
236 if(TESTBIT(fRequirements, kMaxImpactParamR)){
237 // cut on max. Impact Parameter in Radial direction
75d81601 238 if(TMath::Abs(impactR) <= fImpactParamCut[2]) SETBIT(survivedCut, kMaxImpactParamR);
809a4336 239 }
8c1c76e9 240 if(TESTBIT(fRequirements, kMaxImpactParameterRpar)) {
241 // HFE Impact parameter cut
242 GetMaxImpactParameterCutR(track,maximpactRcut);
243 if(TMath::Abs(impactR) >= maximpactRcut) SETBIT(fRequirements, kMaxImpactParameterRpar);
9250ffbf 244 }
809a4336 245 if(TESTBIT(fRequirements, kMaxImpactParamZ)){
246 // cut on max. Impact Parameter in Z direction
75d81601 247 if(TMath::Abs(impactZ) <= fImpactParamCut[3]) SETBIT(survivedCut, kMaxImpactParamZ);
809a4336 248 }
3a72645a 249 if(TESTBIT(fRequirements, kMinHFEImpactParamR)){
250 // cut on min. HFE Impact Parameter in Radial direction
251 if(TMath::Abs(hfeimpactR) >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamR);
252 }
253 if(TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
254 // cut on max. HFE Impact Parameter n sigma in Radial direction
11ff28c5 255 if(fAbsHFEImpactParamNsigmaR) {
256// if((TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) && (TMath::Abs(hfeimpactnsigmaR) < 8)) { //mj debug
257 if(TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) {
258 SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
259// printf("0: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
260 }
261 }
262 else {
263 if(hfeimpactnsigmaRcut>0){
264 if(hfeimpactnsigmaR >= hfeimpactnsigmaRcut) {
265 SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
266 //printf("1: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
267 }
268 }
269 else{
270 if(hfeimpactnsigmaR <= hfeimpactnsigmaRcut) {
271 SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
272 //printf("2: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
273 }
274 }
275 }
3a72645a 276 }
809a4336 277 if(TESTBIT(fRequirements, kClusterRatioTPC)){
278 // cut on min ratio of found TPC clusters vs findable TPC clusters
279 if(ratioTPC >= fClusterRatioTPC) SETBIT(survivedCut, kClusterRatioTPC);
280 }
281 if(TESTBIT(fRequirements, kMinTrackletsTRD)){
282 // cut on minimum number of TRD tracklets
aa1417ee 283 AliDebug(1, Form("Min TRD cut: [%d|%d], Require exact number [%s]\n", fMinTrackletsTRD, trdTracklets, fTRDtrackletsExact ? "Yes" : "No"));
284 if(fTRDtrackletsExact){
e17c1f86 285 if(trdTracklets == fMinTrackletsTRD) {
286 SETBIT(survivedCut, kMinTrackletsTRD);
287 AliDebug(1, "Track Selected");
288 }
aa1417ee 289 }else{
e17c1f86 290 if(trdTracklets >= fMinTrackletsTRD){
291 SETBIT(survivedCut, kMinTrackletsTRD);
292 AliDebug(1, "Track Selected");
293 }
aa1417ee 294 //printf("Min number of tracklets %d\n",fMinTrackletsTRD);
295 }
809a4336 296 }
11ff28c5 297
298 if(TESTBIT(fRequirements, kMinNbITScls)){
299 // cut on minimum number of ITS clusters
300 //printf(Form("Min ITS clusters: [%d|%d]\n", (Int_t)fMinNbITScls, nclsITS));
301 AliDebug(1, Form("Min ITS clusters: [%d|%d]\n", fMinNbITScls, nclsITS));
302 if(nclsITS >= ((Int_t)fMinNbITScls)) SETBIT(survivedCut, kMinNbITScls);
303 }
304
3a72645a 305 if(TESTBIT(fRequirements, kMinNClustersTPC)){
11ff28c5 306 // cut on minimum number of TPC tracklets
307 //printf(Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
3a72645a 308 AliDebug(1, Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
309 if(nclsTPC >= fMinNClustersTPC) SETBIT(survivedCut, kMinNClustersTPC);
310 }
e17c1f86 311 if(TESTBIT(fRequirements, kMinNClustersTPCPID)){
312 AliDebug(1, Form("Min TPC PID cut: [%d|%d]\n", fMinNClustersTPCPID, nclsTPCPID));
313 if(nclsTPCPID >= fMinNClustersTPCPID) SETBIT(survivedCut, kMinNClustersTPCPID);
314 }
809a4336 315 if(TESTBIT(fRequirements, kPixelITS)){
316 // cut on ITS pixel layers
70da6c5a 317 AliDebug(1, "ITS cluster Map: ");
318 //PrintBitMap(itsPixel);
809a4336 319 switch(fPixelITS){
70da6c5a 320 case kFirst:
321 AliDebug(2, "First");
11ff28c5 322 //printf("First\n");
70da6c5a 323 if(TESTBIT(itsPixel, 0))
324 SETBIT(survivedCut, kPixelITS);
325 else
326 if(fCheck && !statusL0)
327 SETBIT(survivedCut, kPixelITS);
dbe3abbe 328 break;
329 case kSecond:
11ff28c5 330 //printf("Second\n");
70da6c5a 331 AliDebug(2, "Second");
332 if(TESTBIT(itsPixel, 1))
333 SETBIT(survivedCut, kPixelITS);
334 else
335 if(fCheck && !statusL1)
336 SETBIT(survivedCut, kPixelITS);
dbe3abbe 337 break;
338 case kBoth:
11ff28c5 339 //printf("Both\n");
70da6c5a 340 AliDebug(2, "Both");
341 if(TESTBIT(itsPixel,0) && TESTBIT(itsPixel,1))
342 SETBIT(survivedCut, kPixelITS);
343 else
344 if(fCheck && !(statusL0 && statusL1))
345 SETBIT(survivedCut, kPixelITS);
346 break;
dbe3abbe 347 case kAny:
11ff28c5 348 //printf("Any\n");
70da6c5a 349 AliDebug(2, "Any");
350 if(TESTBIT(itsPixel, 0) || TESTBIT(itsPixel, 1))
351 SETBIT(survivedCut, kPixelITS);
352 else
353 if(fCheck && !(statusL0 || statusL1))
354 SETBIT(survivedCut, kPixelITS);
dbe3abbe 355 break;
e156c3bb 356 case kExclusiveSecond:
357 AliDebug(2, "Exlusive second");
358 if(fCheck){ // Cut out tracks which pass a dead ITS Layer 0
359 if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0) && statusL0)
360 SETBIT(survivedCut, kPixelITS);
361 } else {
362 if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0))
363 SETBIT(survivedCut, kPixelITS);
364 }
365 break;
11ff28c5 366 case kNone:
367 // No cut applied, set as survived
368 SETBIT(survivedCut, kPixelITS);
369 break;
70da6c5a 370 default:
11ff28c5 371 // default, defined as no cut applied
70da6c5a 372 AliDebug(2, "None");
11ff28c5 373 SETBIT(survivedCut, kPixelITS);
70da6c5a 374 break;
809a4336 375 }
70da6c5a 376 AliDebug(1, Form("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"));
809a4336 377 }
8c1c76e9 378
379 if(TESTBIT(fRequirements, kTOFPID)){
e3ae862b 380 // cut on TOF pid
8c1c76e9 381 if(track->GetStatus() & AliESDtrack::kTOFpid) SETBIT(survivedCut, kTOFPID);
e3ae862b 382 }
c2690925 383
8c1c76e9 384 if(TESTBIT(fRequirements, kTOFmismatch)){
c2690925 385 // cut on TOF mismatch
8c1c76e9 386 if(!(track->GetStatus() & AliESDtrack::kTOFmismatch)) SETBIT(survivedCut, kTOFmismatch);
c2690925 387 }
388
8c1c76e9 389 if(TESTBIT(fRequirements, kTPCPIDCleanUp)){
390 // cut on TPC PID cleanup
8c1c76e9 391 Bool_t fBitsAboveThreshold=GetTPCCountSharedMapBitsAboveThreshold(track);
11ff28c5 392 if((fBitsAboveThreshold==0)&&(nclsTPCPID>80)) SETBIT(survivedCut, kTPCPIDCleanUp);
8c1c76e9 393 }
394
395 if(TESTBIT(fRequirements, kEMCALmatch)){
396 if(track->GetStatus() && AliESDtrack::kEMCALmatch) SETBIT(survivedCut, kEMCALmatch);
397 }
c2690925 398
11ff28c5 399 if(TESTBIT(fRequirements, kRejectKinkDaughter)){
400 //printf("test daughter\n");
401 if(!IsKinkDaughter(track)) SETBIT(survivedCut, kRejectKinkDaughter);
402 //else printf("Found kink daughter\n");
403 }
404
405 if(TESTBIT(fRequirements, kRejectKinkMother)){
406 //printf("test mother\n");
407 if(!IsKinkMother(track)) SETBIT(survivedCut, kRejectKinkMother);
408 //else printf("Found kink mother\n");
409 }
410
8c1c76e9 411 if(fRequirements == survivedCut){
809a4336 412 //
413 // Track selected
414 //
3a72645a 415 AliDebug(2, "Track Survived cuts\n");
416 if(IsQAOn()) FillQAhistosRec(track, kAfterCuts);
809a4336 417 return kTRUE;
418 }
3a72645a 419 AliDebug(2, "Track cut");
809a4336 420 if(IsQAOn()) FillCutCorrelation(survivedCut);
421 return kFALSE;
422}
423
424//______________________________________________________
3a72645a 425Bool_t AliHFEextraCuts::CheckMCCuts(AliVParticle */*track*/) const {
809a4336 426 //
427 // Checks cuts on Monte Carlo tracks
428 // returns true if track is selected
429 // QA histograms are filled before track selection and for
430 // selected tracks after track selection
431 //
432 return kTRUE; // not yet implemented
433}
434
435//______________________________________________________
3a72645a 436void AliHFEextraCuts::FillQAhistosRec(AliVTrack *track, UInt_t when){
809a4336 437 //
438 // Fill the QA histograms for ESD tracks
439 // Function can be called before cuts or after cut application (second argument)
440 //
75d81601 441 Float_t impactR, impactZ;
3a72645a 442 GetImpactParameters(track, impactR, impactZ);
bf892a6a 443 TH1 *htmp = NULL;
8c1c76e9 444 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(0 + when * fgkNQAhistos)))) htmp->Fill(impactR);
445 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(1 + when * fgkNQAhistos)))) htmp->Fill(impactZ);
809a4336 446 // printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
8c1c76e9 447 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(2 + when * fgkNQAhistos)))) htmp->Fill(GetTPCclusterRatio(track));
448 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(3 + when * fgkNQAhistos)))) htmp->Fill(track->GetTRDntrackletsPID());
449 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(4 + when * fgkNQAhistos)))) htmp->Fill(GetTPCncls(track));
809a4336 450 UChar_t itsPixel = track->GetITSClusterMap();
8c1c76e9 451 TH1 *pixelHist = dynamic_cast<TH1F *>(fQAlist->At(5 + when * fgkNQAhistos));
e3fc062d 452 //Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
bf892a6a 453 if(pixelHist){
454 Double_t firstEntry = 0.5;
455 if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
456 pixelHist->Fill(firstEntry + 3);
457 else{
458 if(itsPixel & BIT(0)){
459 pixelHist->Fill(firstEntry);
460 if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
461 else pixelHist->Fill(firstEntry + 4);
462 }
463 if(itsPixel & BIT(1)){
464 pixelHist->Fill(firstEntry + 1);
465 if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
466 }
809a4336 467 }
468 }
8c1c76e9 469 // Fill histogram with the status bits
470 TH1 *hStatusBits = dynamic_cast<TH1 *>(fQAlist->At(6 + when * fgkNQAhistos));
aa1417ee 471 if(hStatusBits) {
472 hStatusBits->Fill(0); // Fill first bin with all tracks
473 if(track->GetStatus() && AliESDtrack::kTOFpid) hStatusBits->Fill(1);
474 if(!(track->GetStatus() && AliESDtrack::kTOFmismatch)) hStatusBits->Fill(2);
475 if(track->GetStatus() && AliESDtrack::kEMCALmatch) hStatusBits->Fill(3);
476 if(GetTPCCountSharedMapBitsAboveThreshold(track)==0) hStatusBits->Fill(4);
477 }
11ff28c5 478 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(7 + when * fgkNQAhistos)))) htmp->Fill(track->GetTPCsignalN());
809a4336 479}
480
481// //______________________________________________________
482// void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
483// //
484// // Fill the QA histograms for MC tracks
485// // Function can be called before cuts or after cut application (second argument)
486// // Not yet implemented
487// //
488// }
489
490//______________________________________________________
491void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
492 //
493 // Fill cut correlation histograms for tracks that didn't pass cuts
494 //
8c1c76e9 495 TH2 *correlation = dynamic_cast<TH2F *>(fQAlist->At(2 * fgkNQAhistos));
bf892a6a 496 if(correlation){
497 for(Int_t icut = 0; icut < kNcuts; icut++){
498 if(!TESTBIT(fRequirements, icut)) continue;
499 for(Int_t jcut = icut; jcut < kNcuts; jcut++){
500 if(!TESTBIT(fRequirements, jcut)) continue;
501 if(TESTBIT(survivedCut, icut) && TESTBIT(survivedCut, jcut))
502 correlation->Fill(icut, jcut);
503 }
809a4336 504 }
505 }
506}
507
508//______________________________________________________
509void AliHFEextraCuts::AddQAHistograms(TList *qaList){
510 //
511 // Add QA histograms
512 // For each cut a histogram before and after track cut is created
513 // Histos before respectively after cut are stored in different lists
514 // Additionally a histogram with the cut correlation is created and stored
515 // in the top directory
516 //
3a72645a 517
809a4336 518 TH1 *histo1D = 0x0;
519 TH2 *histo2D = 0x0;
809a4336 520 TString cutstr[2] = {"before", "after"};
3a72645a 521
522 if(!fQAlist) fQAlist = new TList; // for internal representation, not owner
809a4336 523 for(Int_t icond = 0; icond < 2; icond++){
8c1c76e9 524 qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamR%s",GetName(),cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0 + icond * fgkNQAhistos);
525 fQAlist->AddAt(histo1D, 0 + icond * fgkNQAhistos);
809a4336 526 histo1D->GetXaxis()->SetTitle("Impact Parameter");
527 histo1D->GetYaxis()->SetTitle("Number of Tracks");
8c1c76e9 528 qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamZ%s",GetName(),cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1 + icond * fgkNQAhistos);
529 fQAlist->AddAt(histo1D, 1 + icond * fgkNQAhistos);
809a4336 530 histo1D->GetXaxis()->SetTitle("Impact Parameter");
531 histo1D->GetYaxis()->SetTitle("Number of Tracks");
8c1c76e9 532 qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 100, 0, 1)), 2 + icond * fgkNQAhistos);
533 fQAlist->AddAt(histo1D, 2 + icond * fgkNQAhistos);
809a4336 534 histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
535 histo1D->GetYaxis()->SetTitle("Number of Tracks");
8c1c76e9 536 qaList->AddAt((histo1D = new TH1F(Form("%s_trdTracklets%s",GetName(),cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3 + icond * fgkNQAhistos);
537 fQAlist->AddAt(histo1D, 3 + icond * fgkNQAhistos);
809a4336 538 histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
539 histo1D->GetYaxis()->SetTitle("Number of Tracks");
8c1c76e9 540 qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClusters%s",GetName(),cutstr[icond].Data()), "Number of TPC clusters", 161, 0, 160)), 4 + icond * fgkNQAhistos);
541 fQAlist->AddAt(histo1D, 4 + icond * fgkNQAhistos);
3a72645a 542 histo1D->GetXaxis()->SetTitle("Number of TPC clusters");
543 histo1D->GetYaxis()->SetTitle("Number of Tracks");
8c1c76e9 544 qaList->AddAt((histo1D = new TH1F(Form("%s_itsPixel%s",GetName(),cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 5 + icond * fgkNQAhistos);
545 fQAlist->AddAt(histo1D, 5 + icond * fgkNQAhistos);
809a4336 546 histo1D->GetXaxis()->SetTitle("ITS Pixel");
547 histo1D->GetYaxis()->SetTitle("Number of Tracks");
548 Int_t first = histo1D->GetXaxis()->GetFirst();
549 TString binNames[6] = { "First", "Second", "Both", "None", "Exclusive First", "Exclusive Second"};
550 for(Int_t ilabel = 0; ilabel < 6; ilabel++)
551 histo1D->GetXaxis()->SetBinLabel(first + ilabel, binNames[ilabel].Data());
8c1c76e9 552 qaList->AddAt((histo1D = new TH1F(Form("%s_trackStatus%s",GetName(),cutstr[icond].Data()), "Track Status", 5, 0, 5)), 6 + icond * fgkNQAhistos);
553 fQAlist->AddAt(histo1D, 6 + icond * fgkNQAhistos);
554 histo1D->GetXaxis()->SetTitle("Track Status Bit");
555 histo1D->GetYaxis()->SetTitle("Number of tracks");
556 TString tsnames[5] = {"All", "TOFPID", "No TOFmismatch", "EMCALmatch","No TPC shared bit"};
557 for(Int_t istat = 0; istat < 5; istat++) histo1D->GetXaxis()->SetBinLabel(istat + 1, tsnames[istat].Data());
558 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);
559 fQAlist->AddAt(histo1D, 7 + icond * fgkNQAhistos);
560 histo1D->GetXaxis()->SetTitle("Number of TPC clusters for dEdx calculation");
561 histo1D->GetYaxis()->SetTitle("counts");
809a4336 562 }
809a4336 563 // Add cut correlation
8c1c76e9 564 qaList->AddAt((histo2D = new TH2F(Form("%s_cutcorrelation",GetName()), "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2 * fgkNQAhistos);
565 fQAlist->AddAt(histo2D, 2 * fgkNQAhistos);
566 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 567 Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
568 for(Int_t icut = 0; icut < kNcuts; icut++){
569 histo2D->GetXaxis()->SetBinLabel(firstx + icut, labels[icut].Data());
570 histo2D->GetYaxis()->SetBinLabel(firsty + icut, labels[icut].Data());
571 }
809a4336 572}
dbe3abbe 573
574//______________________________________________________
575void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
576 for(Int_t ibit = 32; ibit--; )
577 printf("%d", bitmap & BIT(ibit) ? 1 : 0);
578 printf("\n");
579}
580
581//______________________________________________________
75d81601 582Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
dbe3abbe 583 //
584 // Check whether ITS area is dead
585 //
586 Bool_t status;
587 switch(itsStatus){
588 case 2: status = kFALSE; break;
589 case 3: status = kFALSE; break;
590 case 7: status = kFALSE; break;
591 default: status = kTRUE;
592 }
593 return status;
594}
3a72645a 595
3a72645a 596//______________________________________________________
597Int_t AliHFEextraCuts::GetITSstatus(AliVTrack *track, Int_t layer){
598 //
599 // Check ITS layer status
600 //
601 Int_t status = 0;
602 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
603 Int_t det;
604 Float_t xloc, zloc;
605 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
bf892a6a 606 if(esdtrack) esdtrack->GetITSModuleIndexInfo(layer, det, status, xloc, zloc);
3a72645a 607 }
608 return status;
609}
610
8c1c76e9 611
612//______________________________________________________
613Bool_t AliHFEextraCuts::GetTPCCountSharedMapBitsAboveThreshold(AliVTrack *track){
614 //
615 // Checks if number of shared bits is above 1
616 //
617 Int_t nsharebit = 1;
11ff28c5 618 const TBits *shared;
619 if(track->IsA() == AliESDtrack::Class())
620 shared = &((static_cast<AliESDtrack *>(track))->GetTPCSharedMap());
621 else if(track->IsA() == AliAODTrack::Class())
622 shared = &((static_cast<AliAODTrack *>(track))->GetTPCSharedMap());
623 else
624 return 0;
625
626 if(shared->CountBits() >= 2) nsharebit=1;
627 else nsharebit=0;
628
8c1c76e9 629 return nsharebit;
630}
631
3a72645a 632//______________________________________________________
e3ae862b 633UInt_t AliHFEextraCuts::GetTPCncls(AliVTrack *track){
3a72645a 634 //
635 // Get Number of findable clusters in the TPC
636 //
3a72645a 637 Int_t nClusters = 0; // in case no Information available consider all clusters findable
11ff28c5 638 TClass *type = track->IsA();
639 if(TESTBIT(fTPCclusterDef, kFoundIter1)){
640 if(type == AliESDtrack::Class()){
641 AliDebug(2, ("Using def kFoundIter1"));
642 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
643 nClusters = esdtrack->GetTPCNclsIter1();
644 } else {
645 AliDebug(2, ("Number of clusters in iteration 1 not available for AOD tracks"));
bf892a6a 646 }
11ff28c5 647 } else if(TESTBIT(fTPCclusterDef, kCrossedRows)){
648 AliDebug(2, ("Using def kCrossedRows"));
649 nClusters = static_cast<UInt_t>(track->GetTPCClusterInfo(2,1));
650 } else{
651 AliDebug(2, ("Using def kFound"));
652 nClusters = track->GetTPCNcls();
3a72645a 653 }
3a72645a 654 return nClusters;
655}
656
9250ffbf 657//______________________________________________________
658Float_t AliHFEextraCuts::GetTPCsharedClustersRatio(AliVTrack *track){
659 //
660 // Get fraction of shared TPC clusters
661 //
662 Float_t fracClustersTPCShared = 0.0;
11ff28c5 663 Int_t nClustersTPC = track->GetTPCNcls();
664 Int_t nClustersTPCShared = 0;
665 TClass *type = track->IsA();
666 if(type == AliESDtrack::Class()){
667 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
668 nClustersTPCShared = esdtrack->GetTPCnclsS();
669 } else if(type == AliESDtrack::Class()){
670 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
671 const TBits &shared = aodtrack->GetTPCSharedMap();
672 nClustersTPCShared = shared.CountBits(0) - shared.CountBits(159);
673 }
674 if (nClustersTPC!=0) {
675 fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
9250ffbf 676 }
677 return fracClustersTPCShared;
678}
679
e3ae862b 680//______________________________________________________
681Double_t AliHFEextraCuts::GetTPCclusterRatio(AliVTrack *track){
682 //
683 // Get Ratio of found / findable clusters for different definitions
684 // Only implemented for ESD tracks
685 //
686 Double_t clusterRatio = 1.; // in case no Information available consider all clusters findable
11ff28c5 687 TClass *type = track->IsA();
688 if(TESTBIT(fTPCclusterRatioDef, kCROverFindable)){
689 AliDebug(2, "Using ratio def kCROverFindable");
690 clusterRatio = track->GetTPCNclsF() ? track->GetTPCClusterInfo(2,1)/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // crossed rows/findable
691 } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverCR)){
692 AliDebug(2, "Using ratio def kFoundOverCR");
693 clusterRatio = track->GetTPCClusterInfo(2,0); // found/crossed rows
694 } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindableIter1)){
695 if(type == AliESDtrack::Class()){
696 AliDebug(2, "Using ratio def kFoundOverFindableIter1");
697 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
698 clusterRatio = esdtrack->GetTPCNclsFIter1() ? static_cast<Double_t>(esdtrack->GetTPCNclsIter1())/static_cast<Double_t>(esdtrack->GetTPCNclsFIter1()) : 1.; // crossed
699 } else {
700 AliDebug(2, "Cluster ratio after iteration 1 not possible for AOD tracks");
701 clusterRatio = 1.;
e3ae862b 702 }
11ff28c5 703 } else {
704 AliDebug(2, "Using ratio def kFoundOverFindable");
705 clusterRatio = track->GetTPCNclsF() ? static_cast<Double_t>(track->GetTPCNcls())/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // found/findable
e3ae862b 706 }
707 return clusterRatio;
708}
709
3a72645a 710//______________________________________________________
711void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z){
11ff28c5 712 //
713 // Get impact parameter
714 //
715 TClass *type = track->IsA();
716 if(type == AliESDtrack::Class()){
717 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
718 esdtrack->GetImpactParameters(radial, z);
719 }
720 else if(type == AliAODTrack::Class()){
721 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
bf892a6a 722 if(aodtrack){
11ff28c5 723 Double_t xyz[3];
724 aodtrack->XYZAtDCA(xyz);
725 z = xyz[2];
726 radial = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]+xyz[1]);
bf892a6a 727 }
11ff28c5 728 }
729}
730//______________________________________________________
731Bool_t AliHFEextraCuts::IsKinkDaughter(AliVTrack *track){
732 //
733 // Is kink Daughter
734 //
735 TClass *type = track->IsA();
736 if(type == AliESDtrack::Class()){
737 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
738 if(!esdtrack) return kFALSE;
739 if(esdtrack->GetKinkIndex(0)<=0) return kFALSE;
740 else return kTRUE;
741
742 }
743 else if(type == AliAODTrack::Class()){
744 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
745 if(aodtrack){
746 AliAODVertex *aodvertex = aodtrack->GetProdVertex();
747 if(!aodvertex) return kFALSE;
748 if(aodvertex->GetType()==AliAODVertex::kKink) return kTRUE;
749 else return kFALSE;
750 }
751 else return kFALSE;
752 }
753
754 return kFALSE;
755}
756//______________________________________________________
757Bool_t AliHFEextraCuts::IsKinkMother(AliVTrack *track){
758 //
759 // Is kink Mother
760 //
761 TClass *type = track->IsA();
762 if(type == AliESDtrack::Class()){
763 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
764 if(!esdtrack) return kFALSE;
765 if(esdtrack->GetKinkIndex(0)!=0) return kTRUE;
766 else return kFALSE;
767 }
768 else if(type == AliAODTrack::Class()){
769 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
770 if(aodtrack){
771 //printf("find a track\n");
772 AliAODVertex *aodvertex = aodtrack->GetProdVertex();
773 if(!aodvertex) return kFALSE;
774 Int_t n = aodvertex->GetNDaughters();
775 //printf("Number of daughters %d\n",n);
776 for(Int_t k=0; k < n; k++) {
777 AliAODTrack *aodtrackdaughter = (AliAODTrack *) aodvertex->GetDaughter(k);
778 if(aodtrackdaughter){
779 AliAODVertex *aodvertexdaughter = aodtrackdaughter->GetProdVertex();
780 if(!aodvertexdaughter) continue;
781 //if(aodvertexdaughter->GetType()==AliAODVertex::kKink) printf("Daughter of type %d\n",(Int_t)aodvertexdaughter->GetType());
782 if(aodvertexdaughter->GetType()==AliAODVertex::kKink) return kTRUE;
3a72645a 783 }
11ff28c5 784 }
785 return kFALSE;
786 }
787 else return kFALSE;
788 }
789
790 return kFALSE;
791
3a72645a 792}
11ff28c5 793//______________________________________________________
794Int_t AliHFEextraCuts::GetITSNbOfcls(AliVTrack *track){
795 //
796 // Get ITS nb of clusters
797 //
798 TClass *type = track->IsA();
799 if(type == AliESDtrack::Class()){
800 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
801 return esdtrack->GetITSclusters(0);
3a72645a 802
11ff28c5 803 }
804 else if(type == AliAODTrack::Class()){
805 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
806 if(aodtrack){
807 return aodtrack->GetITSNcls();
808 }
809 }
810
811 return 0;
812}
3a72645a 813//______________________________________________________
814void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t &dcaxy, Double_t &dcansigmaxy){
815 //
816 // Get HFE impact parameter (with recalculated primary vertex)
817 //
818 dcaxy=0;
819 dcansigmaxy=0;
820 if(!fEvent){
821 AliDebug(1, "No Input event available\n");
822 return;
823 }
824 const Double_t kBeampiperadius=3.;
825 TString type = track->IsA()->GetName();
8c1c76e9 826 Double_t dcaD[2]={-999.,-999.};
827 Double_t covD[3]={-999.,-999.,-999.};
828 Float_t dcaF[2]={-999.,-999.};
829 Float_t covF[3]={-999.,-999.,-999.};
830
831 AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
aa1417ee 832 if(!esdevent) {
833 AliDebug(1, "No esd event available\n");
834 return;
835 }
8c1c76e9 836 const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
837 if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
838 // recalculate primary vertex for peri. and pp
839 AliVertexerTracks vertexer(fEvent->GetMagneticField());
840 vertexer.SetITSMode();
841 vertexer.SetMinClusters(4);
842 Int_t skipped[2];
843 skipped[0] = track->GetID();
844 vertexer.SetSkipTracks(1,skipped);
845
846 //----diamond constraint-----------------------------
847 vertexer.SetConstraintOn();
848 Float_t diamondcovxy[3];
849 esdevent->GetDiamondCovXY(diamondcovxy);
850 Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
851 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
852 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
853 vertexer.SetVtxStart(diamond);
854 delete diamond; diamond=NULL;
855 //----------------------------------------------------
856
857 vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
858
859 // Getting the DCA
860 // Propagation always done on a working copy to not disturb the track params of the original track
861 AliESDtrack *esdtrack = NULL;
862 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
863 // Case ESD track: take copy constructor
864 AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
865 if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
866 } else {
867 // Case AOD track: take different constructor
868 esdtrack = new AliESDtrack(track);
869 }
870 if(esdtrack && esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcaD, covD)){
871 // protection
872 dcaxy = dcaD[0];
873 if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
874 if(!covD[0]) dcansigmaxy = -49.;
875 }
876 delete esdtrack;
877 delete vtxESDSkip;
878 }
879 else{
880 AliESDtrack *esdtrack = NULL;
881 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
3a72645a 882 // Case ESD track: take copy constructor
bf892a6a 883 AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
884 if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
8c1c76e9 885 } else {
3a72645a 886 // Case AOD track: take different constructor
887 esdtrack = new AliESDtrack(track);
8c1c76e9 888 }
889 if(esdtrack) esdtrack->GetImpactParameters(dcaF, covF);
890 dcaxy = dcaF[0];
891 if(covF[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covF[0]);
892 if(!covF[0]) dcansigmaxy = -49.;
893 delete esdtrack;
3a72645a 894 }
3a72645a 895}
896
11ff28c5 897//______________________________________________________
898void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t dcaD[2], Double_t covD[3]){
899 //
900 // Get HFE impact parameter (with recalculated primary vertex)
901 //
902 if(!fEvent){
903 AliDebug(1, "No Input event available\n");
904 return;
905 }
906 const Double_t kBeampiperadius=3.;
907 TString type = track->IsA()->GetName();
908 Float_t dcaF[2]={-999.,-999.};
909 Float_t covF[3]={-999.,-999.,-999.};
910
911 AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
912 if(!esdevent) {
913 AliDebug(1, "No esd event available\n");
914 return;
915 }
916 const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
917 if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
918 // recalculate primary vertex for peri. and pp
919 AliVertexerTracks vertexer(fEvent->GetMagneticField());
920 vertexer.SetITSMode();
921 vertexer.SetMinClusters(4);
922 Int_t skipped[2];
923 skipped[0] = track->GetID();
924 vertexer.SetSkipTracks(1,skipped);
925
926 //----diamond constraint-----------------------------
927 vertexer.SetConstraintOn();
928 Float_t diamondcovxy[3];
929 esdevent->GetDiamondCovXY(diamondcovxy);
930 Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
931 Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
932 AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
933 vertexer.SetVtxStart(diamond);
934 delete diamond; diamond=NULL;
935 //----------------------------------------------------
936
937 vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
938
939 // Getting the DCA
940 // Propagation always done on a working copy to not disturb the track params of the original track
941 AliESDtrack *esdtrack = NULL;
942 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
943 // Case ESD track: take copy constructor
944 AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
945 if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
946 } else {
947 // Case AOD track: take different constructor
948 esdtrack = new AliESDtrack(track);
949 }
950 if(esdtrack) esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcaD, covD);
951 delete esdtrack;
952 delete vtxESDSkip;
953 }
954 else{
955 AliESDtrack *esdtrack = NULL;
956 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
957 // Case ESD track: take copy constructor
958 AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
959 if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
960 } else {
961 // Case AOD track: take different constructor
962 esdtrack = new AliESDtrack(track);
963 }
964 if(esdtrack) esdtrack->GetImpactParameters(dcaF, covF);
965 dcaD[0] = dcaF[0];
966 dcaD[1] = dcaF[1];
967 covD[0] = covF[0];
968 covD[1] = covF[1];
969 covD[2] = covF[2];
970 delete esdtrack;
971 }
972}
3a72645a 973
974//______________________________________________________
975void AliHFEextraCuts::GetHFEImpactParameterCuts(AliVTrack *track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
976 //
977 // Get HFE impact parameter cut(pt dependent)
978 //
979
bf892a6a 980 TString type = track->IsA()->GetName();
981 if(!type.CompareTo("AliESDtrack")){
982 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
983 if(!esdtrack) return;
984 Double_t pt = esdtrack->Pt();
c2690925 985 hfeimpactRcut = fIPcutParam[0]+fIPcutParam[1]*exp(fIPcutParam[2]*pt); // abs R cut
986 hfeimpactnsigmaRcut = fIPcutParam[3]; // sigma cut
3a72645a 987 }
988}
9250ffbf 989//______________________________________________________
990void AliHFEextraCuts::GetMaxImpactParameterCutR(AliVTrack *track, Double_t &maximpactRcut){
991 //
992 // Get max impact parameter cut r (pt dependent)
993 //
994
995 TString type = track->IsA()->GetName();
996 if(!type.CompareTo("AliESDtrack")){
997 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
998 if(!esdtrack) return;
999 Double_t pt = esdtrack->Pt();
1000 if(pt > 0.15) {
1001 maximpactRcut = 0.0182 + 0.035/TMath::Power(pt,1.01); // abs R cut
1002 }
1003 else maximpactRcut = 9999999999.0;
1004 }
1005}