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