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