]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/hfe/AliHFEextraCuts.cxx
Commit modifications done to take care of the problems
[u/mrichter/AliRoot.git] / PWG3 / 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"
35#include "AliESDEvent.h"
36#include "AliESDVertex.h"
809a4336 37#include "AliESDtrack.h"
70da6c5a 38#include "AliLog.h"
809a4336 39#include "AliMCParticle.h"
3a72645a 40#include "AliVEvent.h"
41#include "AliVTrack.h"
42#include "AliVParticle.h"
43#include "AliVertexerTracks.h"
44#include "AliVVertex.h"
809a4336 45
46#include "AliHFEextraCuts.h"
47
48ClassImp(AliHFEextraCuts)
49
50//______________________________________________________
51AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
52 AliCFCutBase(name, title),
3a72645a 53 fEvent(NULL),
809a4336 54 fCutCorrelation(0),
55 fRequirements(0),
3a72645a 56 fTPCiter1(kFALSE),
57 fMinNClustersTPC(0),
809a4336 58 fClusterRatioTPC(0.),
59 fMinTrackletsTRD(0),
60 fPixelITS(0),
70da6c5a 61 fCheck(kFALSE),
62 fQAlist(0x0) ,
dbe3abbe 63 fDebugLevel(0)
809a4336 64{
65 //
66 // Default Constructor
67 //
68 memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
69}
70
71//______________________________________________________
72AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
73 AliCFCutBase(c),
3a72645a 74 fEvent(c.fEvent),
809a4336 75 fCutCorrelation(c.fCutCorrelation),
76 fRequirements(c.fRequirements),
3a72645a 77 fTPCiter1(c.fTPCiter1),
78 fMinNClustersTPC(c.fMinNClustersTPC),
809a4336 79 fClusterRatioTPC(c.fClusterRatioTPC),
80 fMinTrackletsTRD(c.fMinTrackletsTRD),
81 fPixelITS(c.fPixelITS),
75d81601 82 fCheck(c.fCheck),
dbe3abbe 83 fQAlist(0x0),
70da6c5a 84 fDebugLevel(0)
809a4336 85{
86 //
87 // Copy constructor
88 // Performs a deep copy
89 //
90 memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
91 if(IsQAOn()){
92 fIsQAOn = kTRUE;
93 fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
bf892a6a 94 if(fQAlist) fQAlist->SetOwner();
809a4336 95 }
96}
97
98//______________________________________________________
99AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
100 //
101 // Assignment operator
102 //
103 if(this != &c){
104 AliCFCutBase::operator=(c);
3a72645a 105 fEvent = c.fEvent;
809a4336 106 fCutCorrelation = c.fCutCorrelation;
107 fRequirements = c.fRequirements;
3a72645a 108 fTPCiter1 = c.fTPCiter1;
809a4336 109 fClusterRatioTPC = c.fClusterRatioTPC;
3a72645a 110 fMinNClustersTPC = c.fMinNClustersTPC;
809a4336 111 fMinTrackletsTRD = c.fMinTrackletsTRD;
112 fPixelITS = c.fPixelITS;
75d81601 113 fCheck = c.fCheck;
dbe3abbe 114 fDebugLevel = c.fDebugLevel;
809a4336 115 memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
116 if(IsQAOn()){
117 fIsQAOn = kTRUE;
118 fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
bf892a6a 119 if(fQAlist) fQAlist->SetOwner();
809a4336 120 }else fQAlist = 0x0;
121 }
122 return *this;
123}
124
125//______________________________________________________
126AliHFEextraCuts::~AliHFEextraCuts(){
127 //
128 // Destructor
129 //
3a72645a 130 if(fQAlist) delete fQAlist;
131}
132
133//______________________________________________________
134void AliHFEextraCuts::SetRecEventInfo(const TObject *event){
135 //
136 // Set Virtual event an make a copy
137 //
138 if (!event) {
139 AliError("Pointer to AliVEvent !");
140 return;
141 }
142 TString className(event->ClassName());
143 if (! (className.CompareTo("AliESDEvent")==0 || className.CompareTo("AliAODEvent")==0)) {
144 AliError("argument must point to an AliESDEvent or AliAODEvent !");
145 return ;
809a4336 146 }
3a72645a 147 fEvent = (AliVEvent*) event;
148
809a4336 149}
150
151//______________________________________________________
9eeae5d5 152Bool_t AliHFEextraCuts::IsSelected(TObject *o){
809a4336 153 //
154 // Steering function for the track selection
155 //
3a72645a 156 TString type = o->IsA()->GetName();
157 AliDebug(2, Form("Object type %s", type.Data()));
158 if(!type.CompareTo("AliESDtrack") || !type.CompareTo("AliAODTrack")){
159 return CheckRecCuts(dynamic_cast<AliVTrack *>(o));
809a4336 160 }
3a72645a 161 return CheckMCCuts(dynamic_cast<AliVParticle *>(o));
809a4336 162}
163
164//______________________________________________________
3a72645a 165Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
809a4336 166 //
167 // Checks cuts on reconstructed tracks
168 // returns true if track is selected
169 // QA histograms are filled before track selection and for
170 // selected tracks after track selection
171 //
70da6c5a 172 AliDebug(1, "Called");
809a4336 173 ULong64_t survivedCut = 0; // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
3a72645a 174 if(IsQAOn()) FillQAhistosRec(track, kBeforeCuts);
809a4336 175 // Apply cuts
75d81601 176 Float_t impactR, impactZ, ratioTPC;
3a72645a 177 Double_t hfeimpactR, hfeimpactnsigmaR;
178 Double_t hfeimpactRcut, hfeimpactnsigmaRcut;
179 GetImpactParameters(track, impactR, impactZ);
180 if(TESTBIT(fRequirements, kMinHFEImpactParamR) || TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
181 // Protection for PbPb
182 GetHFEImpactParameterCuts(track, hfeimpactRcut, hfeimpactnsigmaRcut);
183 GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
184 }
185 UInt_t nTPCf = GetTPCfindableClusters(track, fTPCiter1), nclsTPC = GetTPCncls(track, fTPCiter1);
809a4336 186 // printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
3a72645a 187 ratioTPC = nTPCf > 0. ? static_cast<Float_t>(nclsTPC)/static_cast<Float_t>(nTPCf) : 1.;
809a4336 188 UChar_t trdTracklets;
3a72645a 189 trdTracklets = GetTRDnTrackletsPID(track);
809a4336 190 UChar_t itsPixel = track->GetITSClusterMap();
3a72645a 191 Int_t status1 = GetITSstatus(track, 0);
192 Int_t status2 = GetITSstatus(track, 1);
70da6c5a 193 Bool_t statusL0 = CheckITSstatus(status1);
194 Bool_t statusL1 = CheckITSstatus(status2);
809a4336 195 if(TESTBIT(fRequirements, kMinImpactParamR)){
196 // cut on min. Impact Parameter in Radial direction
75d81601 197 if(TMath::Abs(impactR) >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
809a4336 198 }
199 if(TESTBIT(fRequirements, kMinImpactParamZ)){
200 // cut on min. Impact Parameter in Z direction
75d81601 201 if(TMath::Abs(impactZ) >= fImpactParamCut[1]) SETBIT(survivedCut, kMinImpactParamZ);
809a4336 202 }
203 if(TESTBIT(fRequirements, kMaxImpactParamR)){
204 // cut on max. Impact Parameter in Radial direction
75d81601 205 if(TMath::Abs(impactR) <= fImpactParamCut[2]) SETBIT(survivedCut, kMaxImpactParamR);
809a4336 206 }
207 if(TESTBIT(fRequirements, kMaxImpactParamZ)){
208 // cut on max. Impact Parameter in Z direction
75d81601 209 if(TMath::Abs(impactZ) <= fImpactParamCut[3]) SETBIT(survivedCut, kMaxImpactParamZ);
809a4336 210 }
3a72645a 211 if(TESTBIT(fRequirements, kMinHFEImpactParamR)){
212 // cut on min. HFE Impact Parameter in Radial direction
213 if(TMath::Abs(hfeimpactR) >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamR);
214 }
215 if(TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
216 // cut on max. HFE Impact Parameter n sigma in Radial direction
217 if(TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
218 }
809a4336 219 if(TESTBIT(fRequirements, kClusterRatioTPC)){
220 // cut on min ratio of found TPC clusters vs findable TPC clusters
221 if(ratioTPC >= fClusterRatioTPC) SETBIT(survivedCut, kClusterRatioTPC);
222 }
223 if(TESTBIT(fRequirements, kMinTrackletsTRD)){
224 // cut on minimum number of TRD tracklets
70da6c5a 225 AliDebug(1, Form("Min TRD cut: [%d|%d]\n", fMinTrackletsTRD, trdTracklets));
809a4336 226 if(trdTracklets >= fMinTrackletsTRD) SETBIT(survivedCut, kMinTrackletsTRD);
227 }
3a72645a 228 if(TESTBIT(fRequirements, kMinNClustersTPC)){
229 // cut on minimum number of TRD tracklets
230 AliDebug(1, Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
231 if(nclsTPC >= fMinNClustersTPC) SETBIT(survivedCut, kMinNClustersTPC);
232 }
809a4336 233 if(TESTBIT(fRequirements, kPixelITS)){
234 // cut on ITS pixel layers
70da6c5a 235 AliDebug(1, "ITS cluster Map: ");
236 //PrintBitMap(itsPixel);
809a4336 237 switch(fPixelITS){
70da6c5a 238 case kFirst:
239 AliDebug(2, "First");
240 if(TESTBIT(itsPixel, 0))
241 SETBIT(survivedCut, kPixelITS);
242 else
243 if(fCheck && !statusL0)
244 SETBIT(survivedCut, kPixelITS);
dbe3abbe 245 break;
246 case kSecond:
70da6c5a 247 AliDebug(2, "Second");
248 if(TESTBIT(itsPixel, 1))
249 SETBIT(survivedCut, kPixelITS);
250 else
251 if(fCheck && !statusL1)
252 SETBIT(survivedCut, kPixelITS);
dbe3abbe 253 break;
254 case kBoth:
70da6c5a 255 AliDebug(2, "Both");
256 if(TESTBIT(itsPixel,0) && TESTBIT(itsPixel,1))
257 SETBIT(survivedCut, kPixelITS);
258 else
259 if(fCheck && !(statusL0 && statusL1))
260 SETBIT(survivedCut, kPixelITS);
261 break;
dbe3abbe 262 case kAny:
70da6c5a 263 AliDebug(2, "Any");
264 if(TESTBIT(itsPixel, 0) || TESTBIT(itsPixel, 1))
265 SETBIT(survivedCut, kPixelITS);
266 else
267 if(fCheck && !(statusL0 || statusL1))
268 SETBIT(survivedCut, kPixelITS);
dbe3abbe 269 break;
70da6c5a 270 default:
271 AliDebug(2, "None");
272 break;
809a4336 273 }
70da6c5a 274 AliDebug(1, Form("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"));
809a4336 275 }
276 if(fRequirements == survivedCut){
277 //
278 // Track selected
279 //
3a72645a 280 AliDebug(2, "Track Survived cuts\n");
281 if(IsQAOn()) FillQAhistosRec(track, kAfterCuts);
809a4336 282 return kTRUE;
283 }
3a72645a 284 AliDebug(2, "Track cut");
809a4336 285 if(IsQAOn()) FillCutCorrelation(survivedCut);
286 return kFALSE;
287}
288
289//______________________________________________________
3a72645a 290Bool_t AliHFEextraCuts::CheckMCCuts(AliVParticle */*track*/) const {
809a4336 291 //
292 // Checks cuts on Monte Carlo tracks
293 // returns true if track is selected
294 // QA histograms are filled before track selection and for
295 // selected tracks after track selection
296 //
297 return kTRUE; // not yet implemented
298}
299
300//______________________________________________________
3a72645a 301void AliHFEextraCuts::FillQAhistosRec(AliVTrack *track, UInt_t when){
809a4336 302 //
303 // Fill the QA histograms for ESD tracks
304 // Function can be called before cuts or after cut application (second argument)
305 //
3a72645a 306 const Int_t kNhistos = 6;
75d81601 307 Float_t impactR, impactZ;
3a72645a 308 GetImpactParameters(track, impactR, impactZ);
309 Int_t nTPCf = GetTPCfindableClusters(track, fTPCiter1), nclsTPC = GetTPCncls(track, fTPCiter1);
bf892a6a 310 TH1 *htmp = NULL;
311 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(0 + when * kNhistos)))) htmp->Fill(impactR);
312 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(1 + when * kNhistos)))) htmp->Fill(impactZ);
809a4336 313 // printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
bf892a6a 314 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(2 + when * kNhistos)))) htmp->Fill(nTPCf > 0. ? static_cast<Float_t>(nclsTPC)/static_cast<Float_t>(nTPCf) : 1.);
315 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(3 + when * kNhistos)))) htmp->Fill(GetTRDnTrackletsPID(track));
316 if((htmp = dynamic_cast<TH1F *>(fQAlist->At(4 + when * kNhistos)))) htmp->Fill(nclsTPC);
809a4336 317 UChar_t itsPixel = track->GetITSClusterMap();
3a72645a 318 TH1 *pixelHist = dynamic_cast<TH1F *>(fQAlist->At(5 + when * kNhistos));
e3fc062d 319 //Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
bf892a6a 320 if(pixelHist){
321 Double_t firstEntry = 0.5;
322 if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
323 pixelHist->Fill(firstEntry + 3);
324 else{
325 if(itsPixel & BIT(0)){
326 pixelHist->Fill(firstEntry);
327 if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
328 else pixelHist->Fill(firstEntry + 4);
329 }
330 if(itsPixel & BIT(1)){
331 pixelHist->Fill(firstEntry + 1);
332 if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
333 }
809a4336 334 }
335 }
336}
337
338// //______________________________________________________
339// void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
340// //
341// // Fill the QA histograms for MC tracks
342// // Function can be called before cuts or after cut application (second argument)
343// // Not yet implemented
344// //
345// }
346
347//______________________________________________________
348void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
349 //
350 // Fill cut correlation histograms for tracks that didn't pass cuts
351 //
3a72645a 352 const Int_t kNhistos = 6;
353 TH2 *correlation = dynamic_cast<TH2F *>(fQAlist->At(2 * kNhistos));
bf892a6a 354 if(correlation){
355 for(Int_t icut = 0; icut < kNcuts; icut++){
356 if(!TESTBIT(fRequirements, icut)) continue;
357 for(Int_t jcut = icut; jcut < kNcuts; jcut++){
358 if(!TESTBIT(fRequirements, jcut)) continue;
359 if(TESTBIT(survivedCut, icut) && TESTBIT(survivedCut, jcut))
360 correlation->Fill(icut, jcut);
361 }
809a4336 362 }
363 }
364}
365
366//______________________________________________________
367void AliHFEextraCuts::AddQAHistograms(TList *qaList){
368 //
369 // Add QA histograms
370 // For each cut a histogram before and after track cut is created
371 // Histos before respectively after cut are stored in different lists
372 // Additionally a histogram with the cut correlation is created and stored
373 // in the top directory
374 //
3a72645a 375
376 const Int_t kNhistos = 6;
809a4336 377 TH1 *histo1D = 0x0;
378 TH2 *histo2D = 0x0;
809a4336 379 TString cutstr[2] = {"before", "after"};
3a72645a 380
381 if(!fQAlist) fQAlist = new TList; // for internal representation, not owner
809a4336 382 for(Int_t icond = 0; icond < 2; icond++){
3a72645a 383 qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamR%s",GetName(),cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0 + icond * kNhistos);
384 fQAlist->AddAt(histo1D, 0 + icond * kNhistos);
809a4336 385 histo1D->GetXaxis()->SetTitle("Impact Parameter");
386 histo1D->GetYaxis()->SetTitle("Number of Tracks");
3a72645a 387 qaList->AddAt((histo1D = new TH1F(Form("%s_impactParamZ%s",GetName(),cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1 + icond * kNhistos);
388 fQAlist->AddAt(histo1D, 1 + icond * kNhistos);
809a4336 389 histo1D->GetXaxis()->SetTitle("Impact Parameter");
390 histo1D->GetYaxis()->SetTitle("Number of Tracks");
3a72645a 391 qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 10, 0, 1)), 2 + icond * kNhistos);
392 fQAlist->AddAt(histo1D, 2 + icond * kNhistos);
809a4336 393 histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
394 histo1D->GetYaxis()->SetTitle("Number of Tracks");
3a72645a 395 qaList->AddAt((histo1D = new TH1F(Form("%s_trdTracklets%s",GetName(),cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3 + icond * kNhistos);
396 fQAlist->AddAt(histo1D, 3 + icond * kNhistos);
809a4336 397 histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
398 histo1D->GetYaxis()->SetTitle("Number of Tracks");
3a72645a 399 qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClusters%s",GetName(),cutstr[icond].Data()), "Number of TPC clusters", 161, 0, 160)), 4 + icond * kNhistos);
400 fQAlist->AddAt(histo1D, 4 + icond * kNhistos);
401 histo1D->GetXaxis()->SetTitle("Number of TPC clusters");
402 histo1D->GetYaxis()->SetTitle("Number of Tracks");
403 qaList->AddAt((histo1D = new TH1F(Form("%s_itsPixel%s",GetName(),cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 5 + icond * kNhistos);
404 fQAlist->AddAt(histo1D, 5 + icond * kNhistos);
809a4336 405 histo1D->GetXaxis()->SetTitle("ITS Pixel");
406 histo1D->GetYaxis()->SetTitle("Number of Tracks");
407 Int_t first = histo1D->GetXaxis()->GetFirst();
408 TString binNames[6] = { "First", "Second", "Both", "None", "Exclusive First", "Exclusive Second"};
409 for(Int_t ilabel = 0; ilabel < 6; ilabel++)
410 histo1D->GetXaxis()->SetBinLabel(first + ilabel, binNames[ilabel].Data());
411 }
809a4336 412 // Add cut correlation
3a72645a 413 qaList->AddAt((histo2D = new TH2F(Form("%s_cutcorrelation",GetName()), "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2 * kNhistos);
414 fQAlist->AddAt(histo2D, 2 * kNhistos);
415 TString labels[kNcuts] = {"MinImpactParamR", "MaxImpactParamR", "MinImpactParamZ", "MaxImpactParamZ", "ClusterRatioTPC", "MinTrackletsTRD", "ITSpixel", "kMinHFEImpactParamR", "kMinHFEImpactParamNsigmaR", "TPC Number of clusters"};
809a4336 416 Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
417 for(Int_t icut = 0; icut < kNcuts; icut++){
418 histo2D->GetXaxis()->SetBinLabel(firstx + icut, labels[icut].Data());
419 histo2D->GetYaxis()->SetBinLabel(firsty + icut, labels[icut].Data());
420 }
809a4336 421}
dbe3abbe 422
423//______________________________________________________
424void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
425 for(Int_t ibit = 32; ibit--; )
426 printf("%d", bitmap & BIT(ibit) ? 1 : 0);
427 printf("\n");
428}
429
430//______________________________________________________
75d81601 431Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
dbe3abbe 432 //
433 // Check whether ITS area is dead
434 //
435 Bool_t status;
436 switch(itsStatus){
437 case 2: status = kFALSE; break;
438 case 3: status = kFALSE; break;
439 case 7: status = kFALSE; break;
440 default: status = kTRUE;
441 }
442 return status;
443}
3a72645a 444
445//______________________________________________________
446Int_t AliHFEextraCuts::GetTRDnTrackletsPID(AliVTrack *track){
447 //
448 // Get Number of TRD tracklets
449 //
450 Int_t nTracklets = 0;
451 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
452 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
bf892a6a 453 if(esdtrack) nTracklets = esdtrack->GetTRDntrackletsPID();
3a72645a 454 } else if(!TString(track->IsA()->GetName()).CompareTo("AliAODTrack")){
455 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
bf892a6a 456 AliAODPid *pidobject = NULL;
457 if(aodtrack) pidobject = aodtrack->GetDetPid();
3a72645a 458 // this is normally NOT the way to do this, but due to limitation in the
459 // AOD track it is not possible in a different way
460 if(pidobject){
461 Float_t *trdmom = pidobject->GetTRDmomentum();
462 for(Int_t ily = 0; ily < 6; ily++){
463 if(trdmom[ily] > -1) nTracklets++;
464 }
465 } else nTracklets = 6; // No Cut possible
466 }
467 return nTracklets;
468}
469
470//______________________________________________________
471Int_t AliHFEextraCuts::GetITSstatus(AliVTrack *track, Int_t layer){
472 //
473 // Check ITS layer status
474 //
475 Int_t status = 0;
476 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
477 Int_t det;
478 Float_t xloc, zloc;
479 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
bf892a6a 480 if(esdtrack) esdtrack->GetITSModuleIndexInfo(layer, det, status, xloc, zloc);
3a72645a 481 }
482 return status;
483}
484
485//______________________________________________________
486Int_t AliHFEextraCuts::GetTPCfindableClusters(AliVTrack *track, Bool_t iter1){
487 //
488 // Get Number of findable clusters in the TPC
489 //
490 AliDebug(1, Form("Using TPC clusters from iteration 1: %s", iter1 ? "Yes" : "No"));
491 Int_t nClusters = 159; // in case no Information available consider all clusters findable
492 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
493 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
bf892a6a 494 if(esdtrack) nClusters = esdtrack->GetTPCNclsF();
3a72645a 495 }
496 return nClusters;
497}
498
499//______________________________________________________
500Int_t AliHFEextraCuts::GetTPCncls(AliVTrack *track, Bool_t iter1){
501 //
502 // Get Number of findable clusters in the TPC
503 //
504 AliDebug(1, Form("Using TPC clusters from iteration 1: %s", iter1 ? "Yes" : "No"));
505 Int_t nClusters = 0; // in case no Information available consider all clusters findable
506 TString type = track->IsA()->GetName();
507 if(!type.CompareTo("AliESDtrack")){
508 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
bf892a6a 509 if(esdtrack){ // coverity
510 if(iter1)
511 nClusters = esdtrack->GetTPCNclsIter1();
512 else
513 nClusters = esdtrack->GetTPCNcls();
514 }
3a72645a 515 }
516 else if(!type.CompareTo("AliAODTrack")){
517 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
bf892a6a 518 const TBits &tpcmap = aodtrack->GetTPCClusterMap();
519 for(UInt_t ibit = 0; ibit < tpcmap.GetNbits(); ibit++)
520 if(tpcmap.TestBitNumber(ibit)) nClusters++;
3a72645a 521 }
522 return nClusters;
523}
524
525//______________________________________________________
526void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z){
527 //
528 // Get impact parameter
529 //
530 TString type = track->IsA()->GetName();
531 if(!type.CompareTo("AliESDtrack")){
532 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
bf892a6a 533 if(esdtrack) esdtrack->GetImpactParameters(radial, z);
3a72645a 534 }
535 else if(!type.CompareTo("AliAODTrack")){
536 AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
bf892a6a 537 if(aodtrack){
538 Double_t xyz[3];
539 aodtrack->XYZAtDCA(xyz);
540 z = xyz[2];
541 radial = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]+xyz[1]);
542 }
3a72645a 543 }
544}
545
546//______________________________________________________
547void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t &dcaxy, Double_t &dcansigmaxy){
548 //
549 // Get HFE impact parameter (with recalculated primary vertex)
550 //
551 dcaxy=0;
552 dcansigmaxy=0;
553 if(!fEvent){
554 AliDebug(1, "No Input event available\n");
555 return;
556 }
557 const Double_t kBeampiperadius=3.;
558 TString type = track->IsA()->GetName();
559 Double_t dca[2]={-999.,-999.};
560 Double_t cov[3]={-999.,-999.,-999.};
561
562 // recalculate primary vertex
563 AliVertexerTracks vertexer(fEvent->GetMagneticField());
564 vertexer.SetITSMode();
565 vertexer.SetMinClusters(4);
566 Int_t skipped[2];
567 skipped[0] = track->GetID();
568 vertexer.SetSkipTracks(1,skipped);
569 AliVVertex *vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
570 vertexer.SetSkipTracks(1,skipped);
571 if(vtxESDSkip->GetNContributors()<2) return;
572
573 // Getting the DCA
574 // Propagation always done on a working copy to not disturb the track params of the original track
575 AliESDtrack *esdtrack = NULL;
576 if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
577 // Case ESD track: take copy constructor
bf892a6a 578 AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
579 if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
3a72645a 580 } else {
581 // Case AOD track: take different constructor
582 esdtrack = new AliESDtrack(track);
583 }
584 if(esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dca, cov)){
585 // protection
586 dcaxy = dca[0];
587 if(cov[0]) dcansigmaxy = dcaxy/TMath::Sqrt(cov[0]);
588 if(!cov[0]) dcansigmaxy = -99.;
589 }
590 delete esdtrack;
591 delete vtxESDSkip;
592}
593
594
595//______________________________________________________
596void AliHFEextraCuts::GetHFEImpactParameterCuts(AliVTrack *track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
597 //
598 // Get HFE impact parameter cut(pt dependent)
599 //
600
bf892a6a 601 TString type = track->IsA()->GetName();
602 if(!type.CompareTo("AliESDtrack")){
603 AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
604 if(!esdtrack) return;
605 Double_t pt = esdtrack->Pt();
606 //hfeimpactRcut=0.0064+0.078*exp(-0.56*pt); // used Carlo's old parameter
607 hfeimpactRcut=0.011+0.077*exp(-0.65*pt); // used Carlo's new parameter
608 hfeimpactnsigmaRcut=3; // 3 sigma trail cut
3a72645a 609 }
610}
611