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