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