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