Since it contains fixes of coding rule violations, all classes are involved. Further...
[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**************************************************************************/
15/************************************************************************
16 * *
17 * Extra cuts implemented by the ALICE Heavy Flavour Electron Group *
18 * Cuts stored here: *
19 * - ITS pixels *
20 * - TPC cluster ratio *
21 * - TRD tracklets *
22 * *
23 * Authors: *
24 * Markus Fasel <M.Fasel@gsi.de> *
25 * *
26 ************************************************************************/
27#include <TAxis.h>
28#include <TClass.h>
29#include <TH1F.h>
30#include <TH2F.h>
31#include <TList.h>
32#include <TString.h>
75d81601 33#include <TMath.h>
809a4336 34
35#include "AliESDtrack.h"
36#include "AliLog.h"
37#include "AliMCParticle.h"
38
39#include "AliHFEextraCuts.h"
40
41ClassImp(AliHFEextraCuts)
42
43//______________________________________________________
44AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
45 AliCFCutBase(name, title),
46 fCutCorrelation(0),
47 fRequirements(0),
48 fClusterRatioTPC(0.),
49 fMinTrackletsTRD(0),
50 fPixelITS(0),
75d81601 51 fCheck(kTRUE),
dbe3abbe 52 fQAlist(0x0),
53 fDebugLevel(0)
809a4336 54{
55 //
56 // Default Constructor
57 //
58 memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
59}
60
61//______________________________________________________
62AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
63 AliCFCutBase(c),
64 fCutCorrelation(c.fCutCorrelation),
65 fRequirements(c.fRequirements),
66 fClusterRatioTPC(c.fClusterRatioTPC),
67 fMinTrackletsTRD(c.fMinTrackletsTRD),
68 fPixelITS(c.fPixelITS),
75d81601 69 fCheck(c.fCheck),
dbe3abbe 70 fQAlist(0x0),
71 fDebugLevel(c.fDebugLevel)
809a4336 72{
73 //
74 // Copy constructor
75 // Performs a deep copy
76 //
77 memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
78 if(IsQAOn()){
79 fIsQAOn = kTRUE;
80 fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
81 fQAlist->SetOwner();
82 }
83}
84
85//______________________________________________________
86AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
87 //
88 // Assignment operator
89 //
90 if(this != &c){
91 AliCFCutBase::operator=(c);
92 fCutCorrelation = c.fCutCorrelation;
93 fRequirements = c.fRequirements;
94 fClusterRatioTPC = c.fClusterRatioTPC;
95 fMinTrackletsTRD = c.fMinTrackletsTRD;
96 fPixelITS = c.fPixelITS;
75d81601 97 fCheck = c.fCheck;
dbe3abbe 98 fDebugLevel = c.fDebugLevel;
809a4336 99
100 memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
101 if(IsQAOn()){
102 fIsQAOn = kTRUE;
103 fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
104 fQAlist->SetOwner();
105 }else fQAlist = 0x0;
106 }
107 return *this;
108}
109
110//______________________________________________________
111AliHFEextraCuts::~AliHFEextraCuts(){
112 //
113 // Destructor
114 //
115 if(fQAlist){
116 fQAlist->Delete();
117 delete fQAlist;
118 }
119}
120
121//______________________________________________________
122Bool_t AliHFEextraCuts::IsSelected(TObject *o){
123 //
124 // Steering function for the track selection
125 //
126 if(TString(o->IsA()->GetName()).CompareTo("AliESDtrack") == 0){
127 return CheckESDCuts(dynamic_cast<AliESDtrack *>(o));
128 }
129 return CheckMCCuts(dynamic_cast<AliMCParticle *>(o));
130}
131
132//______________________________________________________
133Bool_t AliHFEextraCuts::CheckESDCuts(AliESDtrack *track){
134 //
135 // Checks cuts on reconstructed tracks
136 // returns true if track is selected
137 // QA histograms are filled before track selection and for
138 // selected tracks after track selection
139 //
140 ULong64_t survivedCut = 0; // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
141 if(IsQAOn()) FillQAhistosESD(track, kBeforeCuts);
142 // Apply cuts
75d81601 143 Float_t impactR, impactZ, ratioTPC;
144 track->GetImpactParameters(impactR, impactZ);
809a4336 145 // printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
146 ratioTPC = track->GetTPCNclsF() > 0. ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.;
147 UChar_t trdTracklets;
809a4336 148 trdTracklets = track->GetTRDntrackletsPID();
809a4336 149 UChar_t itsPixel = track->GetITSClusterMap();
dbe3abbe 150 Int_t det, status1, status2;
151 Float_t xloc, zloc;
152 track->GetITSModuleIndexInfo(0, det, status1, xloc, zloc);
153 track->GetITSModuleIndexInfo(1, det, status2, xloc, zloc);
809a4336 154 if(TESTBIT(fRequirements, kMinImpactParamR)){
155 // cut on min. Impact Parameter in Radial direction
75d81601 156 if(TMath::Abs(impactR) >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
809a4336 157 }
158 if(TESTBIT(fRequirements, kMinImpactParamZ)){
159 // cut on min. Impact Parameter in Z direction
75d81601 160 if(TMath::Abs(impactZ) >= fImpactParamCut[1]) SETBIT(survivedCut, kMinImpactParamZ);
809a4336 161 }
162 if(TESTBIT(fRequirements, kMaxImpactParamR)){
163 // cut on max. Impact Parameter in Radial direction
75d81601 164 if(TMath::Abs(impactR) <= fImpactParamCut[2]) SETBIT(survivedCut, kMaxImpactParamR);
809a4336 165 }
166 if(TESTBIT(fRequirements, kMaxImpactParamZ)){
167 // cut on max. Impact Parameter in Z direction
75d81601 168 if(TMath::Abs(impactZ) <= fImpactParamCut[3]) SETBIT(survivedCut, kMaxImpactParamZ);
809a4336 169 }
170 if(TESTBIT(fRequirements, kClusterRatioTPC)){
171 // cut on min ratio of found TPC clusters vs findable TPC clusters
172 if(ratioTPC >= fClusterRatioTPC) SETBIT(survivedCut, kClusterRatioTPC);
173 }
174 if(TESTBIT(fRequirements, kMinTrackletsTRD)){
175 // cut on minimum number of TRD tracklets
176 if(trdTracklets >= fMinTrackletsTRD) SETBIT(survivedCut, kMinTrackletsTRD);
177 }
178 if(TESTBIT(fRequirements, kPixelITS)){
179 // cut on ITS pixel layers
dbe3abbe 180 if(fDebugLevel > 0){
181 printf("ITS cluster Map: ");
182 PrintBitMap(itsPixel);
183 }
809a4336 184 switch(fPixelITS){
dbe3abbe 185 case kFirst:
75d81601 186 if(!TESTBIT(itsPixel, 0)) {
187 if(fCheck){
188 if(!CheckITSstatus(status1)) {
189 SETBIT(survivedCut, kPixelITS);
190 }
191 }
192 }
193 else {
194 SETBIT(survivedCut, kPixelITS);
195 }
dbe3abbe 196 break;
197 case kSecond:
75d81601 198 if(!TESTBIT(itsPixel, 1)) {
199 if(fCheck) {
200 if(!CheckITSstatus(status2)) {
201 SETBIT(survivedCut, kPixelITS);
202 }
203 }
204 }
205 else {
206 SETBIT(survivedCut, kPixelITS);
207 }
dbe3abbe 208 break;
209 case kBoth:
75d81601 210 if(!(TESTBIT(track->GetITSClusterMap(),0))) {
211 if(fCheck) {
212 if(!CheckITSstatus(status1)) {
213 if(!(TESTBIT(track->GetITSClusterMap(),1))) {
214 if(!CheckITSstatus(status2)) {
215 SETBIT(survivedCut, kPixelITS);
216 }
217 }
218 else SETBIT(survivedCut, kPixelITS);
219 }
220 }
221 }
222 else {
223
224 if(!(TESTBIT(track->GetITSClusterMap(),1))) {
225 if(fCheck) {
226 if(!CheckITSstatus(status2)) {
227 SETBIT(survivedCut, kPixelITS);
228 }
229 }
230 }
231 else SETBIT(survivedCut, kPixelITS);
232
233 }
234 break;
dbe3abbe 235 case kAny:
75d81601 236 if((!TESTBIT(itsPixel, 0)) && (!TESTBIT(itsPixel, 1))){
237 if(fCheck){
238 if(!CheckITSstatus(status1) || (!CheckITSstatus(status2))) {
239 SETBIT(survivedCut, kPixelITS);
240 }
241 }
242 }
243 else {
244 SETBIT(survivedCut, kPixelITS);
245 }
dbe3abbe 246 break;
809a4336 247 default: break;
248 }
dbe3abbe 249 if(fDebugLevel > 0)printf("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO");
809a4336 250 }
251 if(fRequirements == survivedCut){
252 //
253 // Track selected
254 //
255 if(IsQAOn()) FillQAhistosESD(track, kAfterCuts);
256 return kTRUE;
257 }
258 if(IsQAOn()) FillCutCorrelation(survivedCut);
259 return kFALSE;
260}
261
262//______________________________________________________
75d81601 263Bool_t AliHFEextraCuts::CheckMCCuts(AliMCParticle */*track*/) const {
809a4336 264 //
265 // Checks cuts on Monte Carlo tracks
266 // returns true if track is selected
267 // QA histograms are filled before track selection and for
268 // selected tracks after track selection
269 //
270 return kTRUE; // not yet implemented
271}
272
273//______________________________________________________
274void AliHFEextraCuts::FillQAhistosESD(AliESDtrack *track, UInt_t when){
275 //
276 // Fill the QA histograms for ESD tracks
277 // Function can be called before cuts or after cut application (second argument)
278 //
279 TList *container = dynamic_cast<TList *>(fQAlist->At(when));
75d81601 280 Float_t impactR, impactZ;
281 track->GetImpactParameters(impactR, impactZ);
282 (dynamic_cast<TH1F *>(container->At(0)))->Fill(impactR);
283 (dynamic_cast<TH1F *>(container->At(1)))->Fill(impactZ);
809a4336 284 // printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
285 (dynamic_cast<TH1F *>(container->At(2)))->Fill(track->GetTPCNclsF() > 0. ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.);
809a4336 286 (dynamic_cast<TH1F *>(container->At(3)))->Fill(track->GetTRDntrackletsPID());
809a4336 287 UChar_t itsPixel = track->GetITSClusterMap();
288 TH1 *pixelHist = dynamic_cast<TH1F *>(container->At(4));
289 Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
290 if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
291 pixelHist->Fill(firstEntry + 3);
292 else{
293 if(itsPixel & BIT(0)){
294 pixelHist->Fill(firstEntry);
295 if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
296 else pixelHist->Fill(firstEntry + 4);
297 }
298 if(itsPixel & BIT(1)){
299 pixelHist->Fill(firstEntry + 1);
300 if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
301 }
302 }
303}
304
305// //______________________________________________________
306// void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
307// //
308// // Fill the QA histograms for MC tracks
309// // Function can be called before cuts or after cut application (second argument)
310// // Not yet implemented
311// //
312// }
313
314//______________________________________________________
315void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
316 //
317 // Fill cut correlation histograms for tracks that didn't pass cuts
318 //
319 TH2 *correlation = dynamic_cast<TH2F *>(fQAlist->At(2));
320 for(Int_t icut = 0; icut < kNcuts; icut++){
321 if(!TESTBIT(fRequirements, icut)) continue;
322 for(Int_t jcut = icut; jcut < kNcuts; jcut++){
323 if(!TESTBIT(fRequirements, jcut)) continue;
324 if(TESTBIT(survivedCut, icut) && TESTBIT(survivedCut, jcut))
325 correlation->Fill(icut, jcut);
326 }
327 }
328}
329
330//______________________________________________________
331void AliHFEextraCuts::AddQAHistograms(TList *qaList){
332 //
333 // Add QA histograms
334 // For each cut a histogram before and after track cut is created
335 // Histos before respectively after cut are stored in different lists
336 // Additionally a histogram with the cut correlation is created and stored
337 // in the top directory
338 //
339 TList *histos[2];
340 TH1 *histo1D = 0x0;
341 TH2 *histo2D = 0x0;
342 histos[0] = new TList();
343 histos[0]->SetName("BeforeCut");
344 histos[0]->SetOwner();
345 histos[1] = new TList();
346 histos[1]->SetName("AfterCut");
347 histos[1]->SetOwner();
348 TString cutstr[2] = {"before", "after"};
349 for(Int_t icond = 0; icond < 2; icond++){
350 histos[icond]->AddAt((histo1D = new TH1F(Form("impactParamR%s", cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0);
351 histo1D->GetXaxis()->SetTitle("Impact Parameter");
352 histo1D->GetYaxis()->SetTitle("Number of Tracks");
353 histos[icond]->AddAt((histo1D = new TH1F(Form("impactParamZ%s", cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1);
354 histo1D->GetXaxis()->SetTitle("Impact Parameter");
355 histo1D->GetYaxis()->SetTitle("Number of Tracks");
356 histos[icond]->AddAt((histo1D = new TH1F(Form("tpcClr%s", cutstr[icond].Data()), "Cluster Ratio TPC", 10, 0, 1)), 2);
357 histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
358 histo1D->GetYaxis()->SetTitle("Number of Tracks");
359 histos[icond]->AddAt((histo1D = new TH1F(Form("trdTracklets%s", cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3);
360 histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
361 histo1D->GetYaxis()->SetTitle("Number of Tracks");
362 histos[icond]->AddAt((histo1D = new TH1F(Form("itsPixel%s", cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 4);
363 histo1D->GetXaxis()->SetTitle("ITS Pixel");
364 histo1D->GetYaxis()->SetTitle("Number of Tracks");
365 Int_t first = histo1D->GetXaxis()->GetFirst();
366 TString binNames[6] = { "First", "Second", "Both", "None", "Exclusive First", "Exclusive Second"};
367 for(Int_t ilabel = 0; ilabel < 6; ilabel++)
368 histo1D->GetXaxis()->SetBinLabel(first + ilabel, binNames[ilabel].Data());
369 }
370 fQAlist = new TList();
371 fQAlist->SetOwner();
372 fQAlist->SetName("HFelectronExtraCuts");
373 fQAlist->AddAt(histos[0], 0);
374 fQAlist->AddAt(histos[1], 1);
375 // Add cut correlation
376 fQAlist->AddAt((histo2D = new TH2F("cutcorrellation", "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2);
377 TString labels[kNcuts] = {"MinImpactParamR", "MaxImpactParamR", "MinImpactParamZ", "MaxImpactParamZ", "ClusterRatioTPC", "MinTrackletsTRD", "ITSpixel"};
378 Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
379 for(Int_t icut = 0; icut < kNcuts; icut++){
380 histo2D->GetXaxis()->SetBinLabel(firstx + icut, labels[icut].Data());
381 histo2D->GetYaxis()->SetBinLabel(firsty + icut, labels[icut].Data());
382 }
383 qaList->AddLast(fQAlist);
384}
dbe3abbe 385
386//______________________________________________________
387void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
388 for(Int_t ibit = 32; ibit--; )
389 printf("%d", bitmap & BIT(ibit) ? 1 : 0);
390 printf("\n");
391}
392
393//______________________________________________________
75d81601 394Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
dbe3abbe 395 //
396 // Check whether ITS area is dead
397 //
398 Bool_t status;
399 switch(itsStatus){
400 case 2: status = kFALSE; break;
401 case 3: status = kFALSE; break;
402 case 7: status = kFALSE; break;
403 default: status = kTRUE;
404 }
405 return status;
406}