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