1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // Extra cuts implemented by the ALICE Heavy Flavour Electron Group
19 // - TPC cluster ratio
23 // Markus Fasel <M.Fasel@gsi.de>
32 #include "AliESDtrack.h"
34 #include "AliMCParticle.h"
36 #include "AliHFEextraCuts.h"
38 ClassImp(AliHFEextraCuts)
40 //______________________________________________________
41 AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
42 AliCFCutBase(name, title),
53 // Default Constructor
55 memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
58 //______________________________________________________
59 AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
61 fCutCorrelation(c.fCutCorrelation),
62 fRequirements(c.fRequirements),
63 fClusterRatioTPC(c.fClusterRatioTPC),
64 fMinTrackletsTRD(c.fMinTrackletsTRD),
65 fPixelITS(c.fPixelITS),
72 // Performs a deep copy
74 memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
77 fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
82 //______________________________________________________
83 AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
85 // Assignment operator
88 AliCFCutBase::operator=(c);
89 fCutCorrelation = c.fCutCorrelation;
90 fRequirements = c.fRequirements;
91 fClusterRatioTPC = c.fClusterRatioTPC;
92 fMinTrackletsTRD = c.fMinTrackletsTRD;
93 fPixelITS = c.fPixelITS;
95 fDebugLevel = c.fDebugLevel;
96 memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
99 fQAlist = dynamic_cast<TList *>(c.fQAlist->Clone());
106 //______________________________________________________
107 AliHFEextraCuts::~AliHFEextraCuts(){
117 //______________________________________________________
118 Bool_t AliHFEextraCuts::IsSelected(TObject *o){
120 // Steering function for the track selection
122 if(TString(o->IsA()->GetName()).CompareTo("AliESDtrack") == 0){
123 return CheckESDCuts(dynamic_cast<AliESDtrack *>(o));
125 return CheckMCCuts(dynamic_cast<AliMCParticle *>(o));
128 //______________________________________________________
129 Bool_t AliHFEextraCuts::CheckESDCuts(AliESDtrack *track){
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
136 AliDebug(1, "Called");
137 ULong64_t survivedCut = 0; // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
138 if(IsQAOn()) FillQAhistosESD(track, kBeforeCuts);
140 Float_t impactR, impactZ, ratioTPC;
141 track->GetImpactParameters(impactR, impactZ);
142 // printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
143 ratioTPC = track->GetTPCNclsF() > 0. ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.;
144 UChar_t trdTracklets;
145 trdTracklets = track->GetTRDntrackletsPID();
146 UChar_t itsPixel = track->GetITSClusterMap();
147 Int_t det, status1, status2;
149 track->GetITSModuleIndexInfo(0, det, status1, xloc, zloc);
150 track->GetITSModuleIndexInfo(1, det, status2, xloc, zloc);
151 Bool_t statusL0 = CheckITSstatus(status1);
152 Bool_t statusL1 = CheckITSstatus(status2);
153 if(TESTBIT(fRequirements, kMinImpactParamR)){
154 // cut on min. Impact Parameter in Radial direction
155 if(TMath::Abs(impactR) >= fImpactParamCut[0]) SETBIT(survivedCut, kMinImpactParamR);
157 if(TESTBIT(fRequirements, kMinImpactParamZ)){
158 // cut on min. Impact Parameter in Z direction
159 if(TMath::Abs(impactZ) >= fImpactParamCut[1]) SETBIT(survivedCut, kMinImpactParamZ);
161 if(TESTBIT(fRequirements, kMaxImpactParamR)){
162 // cut on max. Impact Parameter in Radial direction
163 if(TMath::Abs(impactR) <= fImpactParamCut[2]) SETBIT(survivedCut, kMaxImpactParamR);
165 if(TESTBIT(fRequirements, kMaxImpactParamZ)){
166 // cut on max. Impact Parameter in Z direction
167 if(TMath::Abs(impactZ) <= fImpactParamCut[3]) SETBIT(survivedCut, kMaxImpactParamZ);
169 if(TESTBIT(fRequirements, kClusterRatioTPC)){
170 // cut on min ratio of found TPC clusters vs findable TPC clusters
171 if(ratioTPC >= fClusterRatioTPC) SETBIT(survivedCut, kClusterRatioTPC);
173 if(TESTBIT(fRequirements, kMinTrackletsTRD)){
174 // cut on minimum number of TRD tracklets
175 AliDebug(1, Form("Min TRD cut: [%d|%d]\n", fMinTrackletsTRD, trdTracklets));
176 if(trdTracklets >= fMinTrackletsTRD) SETBIT(survivedCut, kMinTrackletsTRD);
178 if(TESTBIT(fRequirements, kPixelITS)){
179 // cut on ITS pixel layers
180 AliDebug(1, "ITS cluster Map: ");
181 //PrintBitMap(itsPixel);
184 AliDebug(2, "First");
185 if(TESTBIT(itsPixel, 0))
186 SETBIT(survivedCut, kPixelITS);
188 if(fCheck && !statusL0)
189 SETBIT(survivedCut, kPixelITS);
192 AliDebug(2, "Second");
193 if(TESTBIT(itsPixel, 1))
194 SETBIT(survivedCut, kPixelITS);
196 if(fCheck && !statusL1)
197 SETBIT(survivedCut, kPixelITS);
201 if(TESTBIT(itsPixel,0) && TESTBIT(itsPixel,1))
202 SETBIT(survivedCut, kPixelITS);
204 if(fCheck && !(statusL0 && statusL1))
205 SETBIT(survivedCut, kPixelITS);
209 if(TESTBIT(itsPixel, 0) || TESTBIT(itsPixel, 1))
210 SETBIT(survivedCut, kPixelITS);
212 if(fCheck && !(statusL0 || statusL1))
213 SETBIT(survivedCut, kPixelITS);
219 AliDebug(1, Form("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"));
221 if(fRequirements == survivedCut){
225 if(IsQAOn()) FillQAhistosESD(track, kAfterCuts);
228 if(IsQAOn()) FillCutCorrelation(survivedCut);
232 //______________________________________________________
233 Bool_t AliHFEextraCuts::CheckMCCuts(AliMCParticle */*track*/) const {
235 // Checks cuts on Monte Carlo tracks
236 // returns true if track is selected
237 // QA histograms are filled before track selection and for
238 // selected tracks after track selection
240 return kTRUE; // not yet implemented
243 //______________________________________________________
244 void AliHFEextraCuts::FillQAhistosESD(AliESDtrack *track, UInt_t when){
246 // Fill the QA histograms for ESD tracks
247 // Function can be called before cuts or after cut application (second argument)
249 TList *container = dynamic_cast<TList *>(fQAlist->At(when));
250 Float_t impactR, impactZ;
251 track->GetImpactParameters(impactR, impactZ);
252 (dynamic_cast<TH1F *>(container->At(0)))->Fill(impactR);
253 (dynamic_cast<TH1F *>(container->At(1)))->Fill(impactZ);
254 // printf("TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
255 (dynamic_cast<TH1F *>(container->At(2)))->Fill(track->GetTPCNclsF() > 0. ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.);
256 (dynamic_cast<TH1F *>(container->At(3)))->Fill(track->GetTRDntrackletsPID());
257 UChar_t itsPixel = track->GetITSClusterMap();
258 TH1 *pixelHist = dynamic_cast<TH1F *>(container->At(4));
259 //Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
260 Double_t firstEntry = 0.5;
261 if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
262 pixelHist->Fill(firstEntry + 3);
264 if(itsPixel & BIT(0)){
265 pixelHist->Fill(firstEntry);
266 if(itsPixel & BIT(1)) pixelHist->Fill(firstEntry + 2);
267 else pixelHist->Fill(firstEntry + 4);
269 if(itsPixel & BIT(1)){
270 pixelHist->Fill(firstEntry + 1);
271 if(!(itsPixel & BIT(0))) pixelHist->Fill(firstEntry + 5);
276 // //______________________________________________________
277 // void AliHFEextraCuts::FillQAhistosMC(AliMCParticle *track, UInt_t when){
279 // // Fill the QA histograms for MC tracks
280 // // Function can be called before cuts or after cut application (second argument)
281 // // Not yet implemented
285 //______________________________________________________
286 void AliHFEextraCuts::FillCutCorrelation(ULong64_t survivedCut){
288 // Fill cut correlation histograms for tracks that didn't pass cuts
290 TH2 *correlation = dynamic_cast<TH2F *>(fQAlist->At(2));
291 for(Int_t icut = 0; icut < kNcuts; icut++){
292 if(!TESTBIT(fRequirements, icut)) continue;
293 for(Int_t jcut = icut; jcut < kNcuts; jcut++){
294 if(!TESTBIT(fRequirements, jcut)) continue;
295 if(TESTBIT(survivedCut, icut) && TESTBIT(survivedCut, jcut))
296 correlation->Fill(icut, jcut);
301 //______________________________________________________
302 void AliHFEextraCuts::AddQAHistograms(TList *qaList){
305 // For each cut a histogram before and after track cut is created
306 // Histos before respectively after cut are stored in different lists
307 // Additionally a histogram with the cut correlation is created and stored
308 // in the top directory
313 histos[0] = new TList();
314 histos[0]->SetName(Form("%s_BeforeCut",GetName()));
315 histos[0]->SetOwner();
316 histos[1] = new TList();
317 histos[1]->SetName(Form("%s_AfterCut",GetName()));
318 histos[1]->SetOwner();
319 TString cutstr[2] = {"before", "after"};
320 for(Int_t icond = 0; icond < 2; icond++){
321 histos[icond]->AddAt((histo1D = new TH1F(Form("%s_impactParamR%s",GetName(),cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0);
322 histo1D->GetXaxis()->SetTitle("Impact Parameter");
323 histo1D->GetYaxis()->SetTitle("Number of Tracks");
324 histos[icond]->AddAt((histo1D = new TH1F(Form("%s_impactParamZ%s",GetName(),cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1);
325 histo1D->GetXaxis()->SetTitle("Impact Parameter");
326 histo1D->GetYaxis()->SetTitle("Number of Tracks");
327 histos[icond]->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 10, 0, 1)), 2);
328 histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
329 histo1D->GetYaxis()->SetTitle("Number of Tracks");
330 histos[icond]->AddAt((histo1D = new TH1F(Form("%s_trdTracklets%s",GetName(),cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3);
331 histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
332 histo1D->GetYaxis()->SetTitle("Number of Tracks");
333 histos[icond]->AddAt((histo1D = new TH1F(Form("%s_itsPixel%s",GetName(),cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 4);
334 histo1D->GetXaxis()->SetTitle("ITS Pixel");
335 histo1D->GetYaxis()->SetTitle("Number of Tracks");
336 Int_t first = histo1D->GetXaxis()->GetFirst();
337 TString binNames[6] = { "First", "Second", "Both", "None", "Exclusive First", "Exclusive Second"};
338 for(Int_t ilabel = 0; ilabel < 6; ilabel++)
339 histo1D->GetXaxis()->SetBinLabel(first + ilabel, binNames[ilabel].Data());
341 fQAlist = new TList();
343 fQAlist->SetName(Form("%s_HFelectronExtraCuts",GetName()));
344 fQAlist->AddAt(histos[0], 0);
345 fQAlist->AddAt(histos[1], 1);
346 // Add cut correlation
347 fQAlist->AddAt((histo2D = new TH2F(Form("%s_cutcorrelation",GetName()), "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2);
348 TString labels[kNcuts] = {"MinImpactParamR", "MaxImpactParamR", "MinImpactParamZ", "MaxImpactParamZ", "ClusterRatioTPC", "MinTrackletsTRD", "ITSpixel"};
349 Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
350 for(Int_t icut = 0; icut < kNcuts; icut++){
351 histo2D->GetXaxis()->SetBinLabel(firstx + icut, labels[icut].Data());
352 histo2D->GetYaxis()->SetBinLabel(firsty + icut, labels[icut].Data());
354 qaList->AddLast(fQAlist);
357 //______________________________________________________
358 void AliHFEextraCuts::PrintBitMap(Int_t bitmap){
359 for(Int_t ibit = 32; ibit--; )
360 printf("%d", bitmap & BIT(ibit) ? 1 : 0);
364 //______________________________________________________
365 Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
367 // Check whether ITS area is dead
371 case 2: status = kFALSE; break;
372 case 3: status = kFALSE; break;
373 case 7: status = kFALSE; break;
374 default: status = kTRUE;