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