Since it contains fixes of coding rule violations, all classes are involved. Further...
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEextraCuts.cxx
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>
33 #include <TMath.h>
34
35 #include "AliESDtrack.h"
36 #include "AliLog.h"
37 #include "AliMCParticle.h"
38
39 #include "AliHFEextraCuts.h"
40
41 ClassImp(AliHFEextraCuts)
42
43 //______________________________________________________
44 AliHFEextraCuts::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),
51   fCheck(kTRUE),
52   fQAlist(0x0),
53   fDebugLevel(0)
54 {
55   //
56   // Default Constructor
57   //
58   memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
59 }
60
61 //______________________________________________________
62 AliHFEextraCuts::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),
69   fCheck(c.fCheck),
70   fQAlist(0x0),
71   fDebugLevel(c.fDebugLevel)
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 //______________________________________________________
86 AliHFEextraCuts &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;
97     fCheck = c.fCheck;
98     fDebugLevel = c.fDebugLevel;
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 //______________________________________________________
111 AliHFEextraCuts::~AliHFEextraCuts(){
112   //
113   // Destructor
114   //
115   if(fQAlist){
116     fQAlist->Delete();
117     delete fQAlist;
118   }
119 }
120
121 //______________________________________________________
122 Bool_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 //______________________________________________________
133 Bool_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
143   Float_t impactR, impactZ, ratioTPC;
144   track->GetImpactParameters(impactR, impactZ);
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;
148   trdTracklets = track->GetTRDntrackletsPID();
149   UChar_t itsPixel = track->GetITSClusterMap();
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);
154   if(TESTBIT(fRequirements, kMinImpactParamR)){
155     // cut on min. Impact Parameter in Radial direction
156     if(TMath::Abs(impactR) >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
157   }
158   if(TESTBIT(fRequirements, kMinImpactParamZ)){
159     // cut on min. Impact Parameter in Z direction
160     if(TMath::Abs(impactZ) >= fImpactParamCut[1]) SETBIT(survivedCut, kMinImpactParamZ);
161   }
162   if(TESTBIT(fRequirements, kMaxImpactParamR)){
163     // cut on max. Impact Parameter in Radial direction
164     if(TMath::Abs(impactR) <= fImpactParamCut[2]) SETBIT(survivedCut, kMaxImpactParamR);
165   }
166   if(TESTBIT(fRequirements, kMaxImpactParamZ)){
167     // cut on max. Impact Parameter in Z direction
168     if(TMath::Abs(impactZ) <= fImpactParamCut[3]) SETBIT(survivedCut, kMaxImpactParamZ);
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
180     if(fDebugLevel > 0){
181       printf("ITS cluster Map: ");
182       PrintBitMap(itsPixel);
183     }
184     switch(fPixelITS){
185       case kFirst: 
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         }
196                     break;
197       case kSecond: 
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         }
208                     break;
209       case kBoth: 
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;
235       case kAny: 
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         }
246                     break;
247       default: break;
248     }
249     if(fDebugLevel > 0)printf("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO");
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 //______________________________________________________
263 Bool_t AliHFEextraCuts::CheckMCCuts(AliMCParticle */*track*/) const {
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 //______________________________________________________
274 void 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));
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);
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.);
286   (dynamic_cast<TH1F *>(container->At(3)))->Fill(track->GetTRDntrackletsPID());
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 //______________________________________________________
315 void 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 //______________________________________________________
331 void 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 }
385
386 //______________________________________________________
387 void 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 //______________________________________________________
394 Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
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 }