This commit was generated by cvs2svn to compensate for changes in r23278,
[u/mrichter/AliRoot.git] / CORRFW / AliCFTrackIsPrimaryCuts.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 // The class AliCFTrackIsPrimaryCut is designed to select reconstructed tracks
17 // with a small impact parameter and tracks which are (not) daughters of kink
18 // decays and to provide corresponding QA histograms.
19 // This class inherits from the Analysis' Framework abstract base class
20 // AliAnalysisCuts and is a part of the Correction Framework.
21 // This class acts on single, reconstructed tracks, it is applicable on
22 // ESD and AOD data.
23 // It mainly consists of a IsSelected function that returns a boolean.
24 // This function checks whether the considered track passes a set of cuts:
25 // - distance to main vertex in units of sigma (resolution)
26 // - require that the dca calculation doesn't fail
27 // - accept or not accept daughter tracks of kink decays
28 //
29 // The cut values for these cuts are set with the corresponding set functions.
30 // All cut classes provided by the correction framework are supposed to be
31 // added in the Analysis Framwork's class AliAnalysisFilter and applied by
32 // the filter via a loop.
33 //
34 // author: I. Kraus (Ingrid.Kraus@cern.ch)
35 // idea taken form
36 // AliESDtrackCuts writte by Jan Fiete Grosse-Oetringhaus and
37 // AliRsnDaughterCut class written by A. Pulvirenti.
38
39 #include <TCanvas.h>
40 #include <TDirectory.h>
41 #include <TBits.h>
42 #include <TH2.h>
43 #include <AliESDtrack.h>
44 #include <AliLog.h>
45 #include "AliCFTrackIsPrimaryCuts.h"
46
47 ClassImp(AliCFTrackIsPrimaryCuts)
48
49 //__________________________________________________________________________________
50 AliCFTrackIsPrimaryCuts::AliCFTrackIsPrimaryCuts() :
51   AliCFCutBase(),
52   fNSigmaToVertexMax(0),
53   fRequireSigmaToVertex(0),
54   fAcceptKinkDaughters(0),
55   fhCutStatistics(0),
56   fhCutCorrelation(0),
57   fBitmap(0x0),
58   fhNBinsNSigma(0),
59   fhNBinsRequireSigma(0),
60   fhNBinsAcceptKink(0),
61   fhNBinsDcaXY(0),
62   fhNBinsDcaZ(0),
63   fhNBinsDcaXYnorm(0),
64   fhNBinsDcaZnorm(0),
65   fhBinLimNSigma(0x0),
66   fhBinLimRequireSigma(0x0),
67   fhBinLimAcceptKink(0x0),
68   fhBinLimDcaXY(0x0),
69   fhBinLimDcaZ(0x0),
70   fhBinLimDcaXYnorm(0x0),
71   fhBinLimDcaZnorm(0x0)
72 {
73   //
74   // Default constructor
75   //
76   fBitmap=new TBits(0);
77   Initialise();
78 }
79 //__________________________________________________________________________________
80 AliCFTrackIsPrimaryCuts::AliCFTrackIsPrimaryCuts(Char_t* name, Char_t* title) :
81   AliCFCutBase(name,title),
82   fNSigmaToVertexMax(0),
83   fRequireSigmaToVertex(0),
84   fAcceptKinkDaughters(0),
85   fhCutStatistics(0),
86   fhCutCorrelation(0),
87   fBitmap(0x0),
88   fhNBinsNSigma(0),
89   fhNBinsRequireSigma(0),
90   fhNBinsAcceptKink(0),
91   fhNBinsDcaXY(0),
92   fhNBinsDcaZ(0),
93   fhNBinsDcaXYnorm(0),
94   fhNBinsDcaZnorm(0),
95   fhBinLimNSigma(0x0),
96   fhBinLimRequireSigma(0x0),
97   fhBinLimAcceptKink(0x0),
98   fhBinLimDcaXY(0x0),
99   fhBinLimDcaZ(0x0),
100   fhBinLimDcaXYnorm(0x0),
101   fhBinLimDcaZnorm(0x0)
102 {
103   //
104   // Constructor
105   //
106   fBitmap=new TBits(0);
107   Initialise();
108 }
109 //__________________________________________________________________________________
110 AliCFTrackIsPrimaryCuts::AliCFTrackIsPrimaryCuts(const AliCFTrackIsPrimaryCuts& c) :
111   AliCFCutBase(c),
112   fNSigmaToVertexMax(c.fNSigmaToVertexMax),
113   fRequireSigmaToVertex(c.fRequireSigmaToVertex),
114   fAcceptKinkDaughters(c.fAcceptKinkDaughters),
115   fhCutStatistics(c.fhCutStatistics),
116   fhCutCorrelation(c.fhCutCorrelation),
117   fBitmap(c.fBitmap),
118   fhNBinsNSigma(c.fhNBinsNSigma),
119   fhNBinsRequireSigma(c.fhNBinsRequireSigma),
120   fhNBinsAcceptKink(c.fhNBinsAcceptKink),
121   fhNBinsDcaXY(c.fhNBinsDcaXY),
122   fhNBinsDcaZ(c.fhNBinsDcaZ),
123   fhNBinsDcaXYnorm(c.fhNBinsDcaXYnorm),
124   fhNBinsDcaZnorm(c.fhNBinsDcaZnorm),
125   fhBinLimNSigma(c.fhBinLimNSigma),
126   fhBinLimRequireSigma(c.fhBinLimRequireSigma),
127   fhBinLimAcceptKink(c.fhBinLimAcceptKink),
128   fhBinLimDcaXY(c.fhBinLimDcaXY),
129   fhBinLimDcaZ(c.fhBinLimDcaZ),
130   fhBinLimDcaXYnorm(c.fhBinLimDcaXYnorm),
131   fhBinLimDcaZnorm(c.fhBinLimDcaZnorm)
132 {
133   //
134   // copy constructor
135   //
136   ((AliCFTrackIsPrimaryCuts &) c).Copy(*this);
137 }
138 //__________________________________________________________________________________
139 AliCFTrackIsPrimaryCuts& AliCFTrackIsPrimaryCuts::operator=(const AliCFTrackIsPrimaryCuts& c)
140 {
141   //
142   // Assignment operator
143   //
144   if (this != &c) {
145     AliCFCutBase::operator=(c) ;
146     fNSigmaToVertexMax = c.fNSigmaToVertexMax ;
147     fRequireSigmaToVertex = c.fRequireSigmaToVertex ;
148     fAcceptKinkDaughters = c.fAcceptKinkDaughters ;
149     fhCutStatistics = c.fhCutStatistics ;
150     fhCutCorrelation = c.fhCutCorrelation ;
151     fBitmap =  c.fBitmap;
152     fhNBinsNSigma = c.fhNBinsNSigma;
153     fhNBinsRequireSigma = c.fhNBinsRequireSigma;
154     fhNBinsAcceptKink = c.fhNBinsAcceptKink;
155     fhNBinsDcaXY = c.fhNBinsDcaXY;
156     fhNBinsDcaZ = c.fhNBinsDcaZ;
157     fhNBinsDcaXYnorm = c.fhNBinsDcaXYnorm;
158     fhNBinsDcaZnorm = c.fhNBinsDcaZnorm;
159     fhBinLimNSigma = c.fhBinLimNSigma;
160     fhBinLimRequireSigma = c.fhBinLimRequireSigma;
161     fhBinLimAcceptKink = c.fhBinLimAcceptKink;
162     fhBinLimDcaXY = c.fhBinLimDcaXY;
163     fhBinLimDcaZ = c.fhBinLimDcaZ;
164     fhBinLimDcaXYnorm = c.fhBinLimDcaXYnorm;
165     fhBinLimDcaZnorm = c.fhBinLimDcaZnorm;
166     
167     for (Int_t i=0; i<c.kNHist; i++){
168       for (Int_t j=0; j<c.kNStepQA; j++){
169         if(c.fhQA[i][j]) fhQA[i][j] = (TH1F*)c.fhQA[i][j]->Clone();
170         if(c.fhDcaXYvsDcaZ[j]) fhDcaXYvsDcaZ[j] = (TH2F*)c.fhDcaXYvsDcaZ[j]->Clone();
171         if(c.fhDcaXYvsDcaZnorm[j]) fhDcaXYvsDcaZnorm[j] = (TH2F*)c.fhDcaXYvsDcaZnorm[j]->Clone();
172       }
173     }
174
175     ((AliCFTrackIsPrimaryCuts &) c).Copy(*this);
176  }
177   return *this;
178 }
179 //__________________________________________________________________________________
180 AliCFTrackIsPrimaryCuts::~AliCFTrackIsPrimaryCuts()
181 {
182   //
183   // destructor
184   //
185   if (fhCutStatistics)                  delete fhCutStatistics;
186   if (fhCutCorrelation)                 delete fhCutCorrelation;
187
188   for (Int_t j=0; j<kNStepQA; j++){
189     for (Int_t i=0; i<kNHist; i++){
190       if(fhQA[i][j]) delete fhQA[i][j];
191     }
192     if(fhDcaXYvsDcaZ[j]) delete fhDcaXYvsDcaZ[j];
193     if(fhDcaXYvsDcaZnorm[j]) delete fhDcaXYvsDcaZnorm[j];
194   }
195
196   if (fBitmap)  delete fBitmap;
197   
198   if(fhBinLimNSigma) delete fhBinLimNSigma;
199   if(fhBinLimRequireSigma) delete fhBinLimRequireSigma;
200   if(fhBinLimAcceptKink) delete fhBinLimAcceptKink;
201   if(fhBinLimDcaXY) delete fhBinLimDcaXY;
202   if(fhBinLimDcaZ) delete fhBinLimDcaZ;
203   if(fhBinLimDcaXYnorm) delete fhBinLimDcaXYnorm;
204   if(fhBinLimDcaZnorm) delete fhBinLimDcaZnorm;
205 }
206 //__________________________________________________________________________________
207 void AliCFTrackIsPrimaryCuts::Initialise()
208 {
209   //
210   // sets everything to zero
211   //
212   fNSigmaToVertexMax = 0;
213   fRequireSigmaToVertex = 0;
214   fAcceptKinkDaughters = 0;
215
216   SetMaxNSigmaToVertex();
217   SetRequireSigmaToVertex();
218
219   for (Int_t j=0; j<kNStepQA; j++)  {
220     for (Int_t i=0; i<kNHist; i++){
221       fhQA[i][j] = 0x0;
222     }
223     fhDcaXYvsDcaZ[j] = 0x0;
224     fhDcaXYvsDcaZnorm[j] = 0x0;
225   }
226   fhCutStatistics = 0;
227   fhCutCorrelation = 0;
228   
229   //set default bining for QA histograms
230   SetHistogramBins(kCutNSigmaToVertex,500,0,50);
231   SetHistogramBins(kCutRequireSigmaToVertex,2,-0.5,1.5);
232   SetHistogramBins(kCutAcceptKinkDaughters,2,-0.5,1.5);
233   SetHistogramBins(kDcaXY,500,-10,10);
234   SetHistogramBins(kDcaZ,500,-10,10);
235   SetHistogramBins(kDcaXYnorm,500,-10,10);
236   SetHistogramBins(kDcaZnorm,500,-10,10);
237 }
238 //__________________________________________________________________________________
239 void AliCFTrackIsPrimaryCuts::Copy(TObject &c) const
240 {
241   //
242   // Copy function
243   //
244   AliCFTrackIsPrimaryCuts& target = (AliCFTrackIsPrimaryCuts &) c;
245
246   target.Initialise();
247
248   if (fhCutStatistics)  target.fhCutStatistics = (TH1F*) fhCutStatistics->Clone();
249   if (fhCutCorrelation) target.fhCutCorrelation = (TH2F*) fhCutCorrelation->Clone();
250
251   for (Int_t j=0; j<kNStepQA; j++){
252     for (Int_t i=0; i<kNHist; i++){
253       if(fhQA[i][j]) target.fhQA[i][j] = (TH1F*)fhQA[i][j]->Clone();
254       
255     }
256     if(fhDcaXYvsDcaZ[j]) target.fhDcaXYvsDcaZ[j] = (TH2F*)fhDcaXYvsDcaZ[j]->Clone();
257     if(fhDcaXYvsDcaZnorm[j]) target.fhDcaXYvsDcaZnorm[j] = (TH2F*)fhDcaXYvsDcaZnorm[j]->Clone();
258   }
259   
260   TNamed::Copy(c);
261 }
262 //____________________________________________________________________
263 Float_t AliCFTrackIsPrimaryCuts::GetSigmaToVertex(AliESDtrack* esdTrack) const
264 {
265   //
266   // Calculates the number of sigma to the vertex.
267   //
268   Float_t b[2];
269   Float_t bRes[2];
270   Float_t bCov[3];
271   esdTrack->GetImpactParameters(b,bCov);
272   if (bCov[0]<=0 || bCov[2]<=0) {
273     AliDebug(1, "Estimated b resolution lower or equal zero!");
274     bCov[0]=0; bCov[2]=0;
275   }
276   bRes[0] = TMath::Sqrt(bCov[0]);
277   bRes[1] = TMath::Sqrt(bCov[2]);
278
279   // -----------------------------------
280   // How to get to a n-sigma cut?
281   //
282   // The accumulated statistics from 0 to d is
283   //
284   // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
285   // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
286   //
287   // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
288   // Can this be expressed in a different way?
289
290   if (bRes[0] == 0 || bRes[1] ==0)
291     return -1;
292
293   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
294
295   // stupid rounding problem screws up everything:
296   // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
297   if (TMath::Exp(-d * d / 2) < 1e-10)
298     return 1000;
299
300   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
301   return d;
302 }
303 //__________________________________________________________________________________
304 void AliCFTrackIsPrimaryCuts::GetBitMap(TObject* obj, TBits *bitmap)  {
305   //
306   // retrieve the pointer to the bitmap
307   //
308   bitmap = SelectionBitMap(obj);
309 }
310 //__________________________________________________________________________________
311 TBits* AliCFTrackIsPrimaryCuts::SelectionBitMap(TObject* obj)
312 {
313   //
314   // test if the track passes the single cuts
315   // and store the information in a bitmap
316   //
317
318   // bitmap stores the decision of each single cut
319   for(Int_t i=0; i<kNCuts; i++)fBitmap->SetBitNumber(i,kFALSE);
320
321   // cast TObject into ESDtrack
322   AliESDtrack* esdTrack = dynamic_cast<AliESDtrack *>(obj);
323   if ( !esdTrack ) return fBitmap ;
324
325
326   for(Int_t i=0; i<kNCuts; i++)fBitmap->SetBitNumber(i,kTRUE);
327
328   // get the track to vertex parameter
329   Float_t nSigmaToVertex = GetSigmaToVertex(esdTrack);
330
331   // fill the bitmap
332   Int_t iCutBit = 0;
333   if (nSigmaToVertex > fNSigmaToVertexMax)
334         fBitmap->SetBitNumber(iCutBit,kFALSE);
335   iCutBit++;
336   if (nSigmaToVertex<0 && fRequireSigmaToVertex)
337         fBitmap->SetBitNumber(iCutBit,kFALSE);
338   iCutBit++;
339   if (!fAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
340         fBitmap->SetBitNumber(iCutBit,kFALSE);
341
342   return fBitmap;
343 }
344 //__________________________________________________________________________________
345 Bool_t AliCFTrackIsPrimaryCuts::IsSelected(TObject* obj) {
346   //
347   // loops over decisions of single cuts and returns if the track is accepted
348   //
349   TBits* bitmap = SelectionBitMap(obj);
350
351   Bool_t isSelected = kTRUE;
352
353   for (UInt_t icut=0; icut<bitmap->GetNbits();icut++)
354         if(!bitmap->TestBitNumber(icut)) isSelected = kFALSE;
355
356   return isSelected;
357 }
358 //__________________________________________________________________________________
359 void AliCFTrackIsPrimaryCuts::Init() {
360   //
361   // initialises all histograms and the TList which holds the histograms
362   //
363   if(fIsQAOn)
364     DefineHistograms();
365 }
366 //__________________________________________________________________________________
367 void AliCFTrackIsPrimaryCuts::SetHistogramBins(Int_t index, Int_t nbins, Double_t *bins)
368 {
369   //
370   // variable bin size
371   //
372   if(!fIsQAOn) return;
373
374   switch(index){
375   case kCutNSigmaToVertex:
376     fhNBinsNSigma=nbins;
377     fhBinLimNSigma=new Double_t[nbins+1];
378     for(Int_t i=0;i<nbins+1;i++)fhBinLimNSigma[i]=bins[i];
379     break;
380
381   case kCutRequireSigmaToVertex:
382     fhNBinsRequireSigma=nbins;
383     fhBinLimRequireSigma=new Double_t[nbins+1];
384     for(Int_t i=0;i<nbins+1;i++)fhBinLimRequireSigma[i]=bins[i];
385     break;
386
387   case kCutAcceptKinkDaughters:
388     fhNBinsAcceptKink=nbins;
389     fhBinLimAcceptKink=new Double_t[nbins+1];
390     for(Int_t i=0;i<nbins+1;i++)fhBinLimAcceptKink[i]=bins[i];
391     break;
392
393   case kDcaXY:
394     fhNBinsDcaXY=nbins;
395     fhBinLimDcaXY=new Double_t[nbins+1];
396     for(Int_t i=0;i<nbins+1;i++)fhBinLimDcaXY[i]=bins[i];
397     break;
398     
399   case kDcaZ:
400     fhNBinsDcaZ=nbins;
401     fhBinLimDcaZ=new Double_t[nbins+1];
402     for(Int_t i=0;i<nbins+1;i++)fhBinLimDcaZ[i]=bins[i];
403     break;
404     
405   case kDcaXYnorm:
406     fhNBinsDcaXYnorm=nbins;
407     fhBinLimDcaXYnorm=new Double_t[nbins+1];
408     for(Int_t i=0;i<nbins+1;i++)fhBinLimDcaXYnorm[i]=bins[i];
409     break;
410     
411   case kDcaZnorm:
412     fhNBinsDcaZnorm=nbins;
413     fhBinLimDcaZnorm=new Double_t[nbins+1];
414     for(Int_t i=0;i<nbins+1;i++)fhBinLimDcaZnorm[i]=bins[i];
415     break;
416   }
417 }
418 //__________________________________________________________________________________
419 void AliCFTrackIsPrimaryCuts::SetHistogramBins(Int_t index, Int_t nbins, Double_t xmin, Double_t xmax)
420 {
421   //
422   // fixed bin size
423   //
424   if(!fIsQAOn) return;
425
426   switch(index){
427   case kCutNSigmaToVertex:
428     fhNBinsNSigma=nbins;
429     fhBinLimNSigma=new Double_t[nbins+1];
430     for(Int_t i=0;i<nbins+1;i++)fhBinLimNSigma[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
431     break;
432
433   case kCutRequireSigmaToVertex:
434     fhNBinsRequireSigma=nbins;
435     fhBinLimRequireSigma=new Double_t[nbins+1];
436     for(Int_t i=0;i<nbins+1;i++)fhBinLimRequireSigma[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
437     break;
438
439   case kCutAcceptKinkDaughters:
440     fhNBinsAcceptKink=nbins;
441     fhBinLimAcceptKink=new Double_t[nbins+1];
442     for(Int_t i=0;i<nbins+1;i++)fhBinLimAcceptKink[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
443     break;
444
445   case kDcaXY:
446     fhNBinsDcaXY=nbins;
447     fhBinLimDcaXY=new Double_t[nbins+1];
448     for(Int_t i=0;i<nbins+1;i++)fhBinLimDcaXY[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
449     break;
450     
451   case kDcaZ:
452     fhNBinsDcaZ=nbins;
453     fhBinLimDcaZ=new Double_t[nbins+1];
454     for(Int_t i=0;i<nbins+1;i++)fhBinLimDcaZ[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
455     break;
456     
457   case kDcaXYnorm:
458     fhNBinsDcaXYnorm=nbins;
459     fhBinLimDcaXYnorm=new Double_t[nbins+1];
460     for(Int_t i=0;i<nbins+1;i++)fhBinLimDcaXYnorm[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
461     break;
462     
463   case kDcaZnorm:
464     fhNBinsDcaZnorm=nbins;
465     fhBinLimDcaZnorm=new Double_t[nbins+1];
466     for(Int_t i=0;i<nbins+1;i++)fhBinLimDcaZnorm[i]=xmin+i*(xmax-xmin)/Double_t(nbins);
467     break;
468   }
469 }
470 //__________________________________________________________________________________
471  void AliCFTrackIsPrimaryCuts::DefineHistograms() {
472   //
473   // histograms for cut variables, cut statistics and cut correlations
474   //
475
476   Int_t color = 2;
477
478   // book cut statistics and cut correlation histograms
479   fhCutStatistics = new TH1F(Form("%s_cut_statistics",GetName()), Form("%s cut statistics",GetName()), kNCuts,0.5,kNCuts+0.5);
480   fhCutStatistics->SetLineWidth(2);
481   fhCutStatistics->GetXaxis()->SetBinLabel(1,"n dca");
482   fhCutStatistics->GetXaxis()->SetBinLabel(2,"require dca");
483   fhCutStatistics->GetXaxis()->SetBinLabel(3,"kink daughter");
484
485   fhCutCorrelation = new TH2F(Form("%s_cut_correlation",GetName()), Form("%s cut  correlation",GetName()), kNCuts,0.5,kNCuts+0.5,kNCuts,0.5,kNCuts+0.5);
486   fhCutCorrelation->SetLineWidth(2);
487   fhCutCorrelation->GetXaxis()->SetBinLabel(1,"n dca");
488   fhCutCorrelation->GetXaxis()->SetBinLabel(2,"require dca");
489   fhCutCorrelation->GetXaxis()->SetBinLabel(3,"kink daughter");
490
491   fhCutCorrelation->GetYaxis()->SetBinLabel(1,"n dca");
492   fhCutCorrelation->GetYaxis()->SetBinLabel(2,"require dca");
493   fhCutCorrelation->GetYaxis()->SetBinLabel(3,"kink daughter");
494
495   // book QA histograms
496   Char_t str[256];
497   for (Int_t i=0; i<kNStepQA; i++) {
498     if (i==0) sprintf(str," ");
499     else sprintf(str,"_cut");
500   
501     fhDcaXYvsDcaZ[i]            = new  TH2F(Form("%s_dcaXYvsDcaZ%s",GetName(),str),"",200,-10,10,200,-10,10);
502     fhDcaXYvsDcaZnorm[i]        = new  TH2F(Form("%s_dcaXYvsDcaZnorm%s",GetName(),str),"",200,-10,10,200,-10,10);
503
504     fhQA[kCutNSigmaToVertex][i] = new TH1F(Form("%s_nSigmaToVertex%s",GetName(),str),"",fhNBinsNSigma,fhBinLimNSigma);
505     fhQA[kCutRequireSigmaToVertex][i] = new TH1F(Form("%s_requireSigmaToVertex%s",GetName(),str),"",fhNBinsRequireSigma,fhBinLimRequireSigma);
506     fhQA[kCutAcceptKinkDaughters][i] = new TH1F(Form("%s_acceptKinkDaughters%s",GetName(),str),"",fhNBinsAcceptKink,fhBinLimAcceptKink);
507     fhQA[kDcaXY][i]             = new TH1F(Form("%s_dcaXY%s",GetName(),str),"",fhNBinsDcaXY,fhBinLimDcaXY);
508     fhQA[kDcaZ][i]              = new TH1F(Form("%s_dcaZ%s",GetName(),str),"",fhNBinsDcaZ,fhBinLimDcaZ);
509     fhQA[kDcaXYnorm][i]         = new TH1F(Form("%s_dcaXYnorm%s",GetName(),str),"",fhNBinsDcaXYnorm,fhBinLimDcaXYnorm);
510     fhQA[kDcaZnorm][i]          = new TH1F(Form("%s_dcaZnorm%s",GetName(),str),"",fhNBinsDcaZnorm,fhBinLimDcaZnorm);
511
512     fhDcaXYvsDcaZ[i]->SetXTitle("impact par. d_{z}");
513     fhDcaXYvsDcaZ[i]->SetYTitle("impact par. d_{xy}");
514     fhDcaXYvsDcaZnorm[i]->SetXTitle("norm. impact par. d_{z} / #sigma_{z}");
515     fhDcaXYvsDcaZnorm[i]->SetYTitle("norm. impact par. d_{xy} / #sigma_{xy}");
516
517     fhQA[kCutNSigmaToVertex][i]->SetXTitle("n #sigma to vertex");
518     fhQA[kCutRequireSigmaToVertex][i]->SetXTitle("require #sigma to vertex");
519     fhQA[kCutAcceptKinkDaughters][i]->SetXTitle("accept kink daughters");
520     fhQA[kDcaXY][i]->SetXTitle("impact par. d_{xy}");
521     fhQA[kDcaZ][i]->SetXTitle("impact par. d_{z}");
522     fhQA[kDcaXYnorm][i]->SetXTitle("norm. impact par. d_{xy} / #sigma_{xy}");
523     fhQA[kDcaZnorm][i]->SetXTitle("norm. impact par. d_{z} / #sigma_{z}");
524   }
525
526   for(Int_t i=0; i<kNHist; i++) fhQA[i][1]->SetLineColor(color);
527
528 }
529 //__________________________________________________________________________________
530 void AliCFTrackIsPrimaryCuts::FillHistograms(TObject* obj, Bool_t f)
531 {
532   //
533   // fill the QA histograms
534   //
535   if(!fIsQAOn) return;
536
537   // cast TObject into ESDtrack
538   AliESDtrack* esdTrack = dynamic_cast<AliESDtrack *>(obj);
539   if ( !esdTrack ) return;
540
541   // index = 0: fill histograms before cuts
542   // index = 1: fill histograms after cuts
543   Int_t index = -1;
544   index = ((f) ? 1 : 0);
545
546         Float_t b[2];
547         Float_t bRes[2];
548         Float_t bCov[3];
549         esdTrack->GetImpactParameters(b,bCov);
550         if (bCov[0]<=0 || bCov[2]<=0) {
551            AliDebug(1, "Estimated b resolution lower or equal zero!");
552            bCov[0]=0; bCov[2]=0;
553         }
554         bRes[0] = TMath::Sqrt(bCov[0]);
555         bRes[1] = TMath::Sqrt(bCov[2]);
556
557         fhQA[kDcaZ][index]->Fill(b[1]);
558         fhQA[kDcaXY][index]->Fill(b[0]);
559         fhDcaXYvsDcaZ[index]->Fill(b[1],b[0]);
560
561         if (bRes[0]!=0 && bRes[1]!=0) {
562            fhQA[kDcaZnorm][index]->Fill(b[1]/bRes[1]);
563            fhQA[kDcaXYnorm][index]->Fill(b[0]/bRes[0]);
564            fhDcaXYvsDcaZnorm[index]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
565         }
566
567   // getting the track to vertex parameters
568   Float_t nSigmaToVertex = GetSigmaToVertex(esdTrack);
569
570   fhQA[kCutNSigmaToVertex][index]->Fill(nSigmaToVertex);
571   if (nSigmaToVertex<0 && fRequireSigmaToVertex) fhQA[kCutRequireSigmaToVertex][index]->Fill(0.);
572   if (!(nSigmaToVertex<0 && fRequireSigmaToVertex)) fhQA[kCutRequireSigmaToVertex][index]->Fill(1.);
573
574
575   if (!fAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0) fhQA[kCutAcceptKinkDaughters][index]->Fill(0.);
576   if (!(!fAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)) fhQA[kCutAcceptKinkDaughters][index]->Fill(0.);
577
578   // fill cut statistics and cut correlation histograms with information from the bitmap
579   if (f) return;
580
581   // Get the bitmap of the single cuts
582   if ( !obj ) return;
583   TBits* bitmap = SelectionBitMap(obj);
584
585   // Number of single cuts in this class
586   UInt_t ncuts = bitmap->GetNbits();
587   for(UInt_t bit=0; bit<ncuts;bit++) {
588     if (!bitmap->TestBitNumber(bit)) {
589         fhCutStatistics->Fill(bit+1);
590         for (UInt_t bit2=bit; bit2<ncuts;bit2++) {
591           if (!bitmap->TestBitNumber(bit2)) 
592             fhCutCorrelation->Fill(bit+1,bit2+1);
593         }
594     }
595   }
596 }
597 //__________________________________________________________________________________
598 void AliCFTrackIsPrimaryCuts::SaveHistograms(const Char_t* dir) {
599   //
600   // saves the histograms in a directory (dir)
601   //
602   if(!fIsQAOn) return;
603
604   if (!dir)
605     dir = GetName();
606
607   gDirectory->mkdir(dir);
608   gDirectory->cd(dir);
609
610   gDirectory->mkdir("before_cuts");
611   gDirectory->mkdir("after_cuts");
612
613   fhCutStatistics->Write();
614   fhCutCorrelation->Write();
615
616   for (Int_t j=0; j<kNStepQA; j++) {
617     if (j==0)
618       gDirectory->cd("before_cuts");
619     else
620       gDirectory->cd("after_cuts");
621
622     fhDcaXYvsDcaZ[j]    ->Write();
623     fhDcaXYvsDcaZnorm[j]->Write();
624
625     for(Int_t i=0; i<kNHist; i++) fhQA[i][j]->Write();
626
627     gDirectory->cd("../");
628   }
629   gDirectory->cd("../");
630 }
631 //__________________________________________________________________________________
632 void AliCFTrackIsPrimaryCuts::DrawHistograms()
633 {
634   //
635   // draws some histograms
636   //
637   if(!fIsQAOn) return;
638
639   // pad margins
640   Float_t right = 0.03;
641   Float_t left = 0.175;
642   Float_t top = 0.03;
643   Float_t bottom = 0.175;
644
645   TCanvas* canvas1 = new TCanvas("Track_QA_Primary_1", "Track QA Primary 1", 800, 500);
646   canvas1->Divide(2, 1);
647
648   canvas1->cd(1);
649   fhCutStatistics->SetStats(kFALSE);
650   fhCutStatistics->LabelsOption("v");
651   gPad->SetLeftMargin(left);
652   gPad->SetBottomMargin(0.25);
653   gPad->SetRightMargin(right);
654   gPad->SetTopMargin(0.1);
655   fhCutStatistics->Draw();
656
657   canvas1->cd(2);
658   fhCutCorrelation->SetStats(kFALSE);
659   fhCutCorrelation->LabelsOption("v");
660   gPad->SetLeftMargin(0.30);
661   gPad->SetRightMargin(bottom);
662   gPad->SetTopMargin(0.1);
663   gPad->SetBottomMargin(0.25);
664   fhCutCorrelation->Draw("COLZ");
665
666   canvas1->SaveAs(Form("%s.eps", canvas1->GetName()));
667   canvas1->SaveAs(Form("%s.ps", canvas1->GetName()));
668
669   // -----
670
671   TCanvas* canvas2 = new TCanvas("Track_QA_Primary_2", "Track QA Primary 2", 800, 800);
672   canvas2->Divide(2, 2);
673
674   canvas2->cd(1);
675   gPad->SetRightMargin(right);
676   gPad->SetLeftMargin(left);
677   gPad->SetTopMargin(top);
678   gPad->SetBottomMargin(bottom);
679   fhQA[kDcaXY][0]->SetStats(kFALSE);
680   fhQA[kDcaXY][0]->Draw();
681   fhQA[kDcaXY][1]->Draw("same");
682
683   canvas2->cd(2);
684   gPad->SetRightMargin(right);
685   gPad->SetLeftMargin(left);
686   gPad->SetTopMargin(top);
687   gPad->SetBottomMargin(bottom);
688   fhQA[kDcaZ][0]->SetStats(kFALSE);
689   fhQA[kDcaZ][0]->Draw();
690   fhQA[kDcaZ][1]->Draw("same");
691
692   canvas2->cd(3);
693 //   fhDXYvsDZ[0]->SetStats(kFALSE);
694   gPad->SetLogz();
695   gPad->SetLeftMargin(bottom);
696   gPad->SetTopMargin(0.1);
697   gPad->SetBottomMargin(bottom);
698   gPad->SetRightMargin(0.2);
699   fhDcaXYvsDcaZ[0]->Draw("COLZ");
700
701   canvas2->cd(4);
702 //   fhDXYvsDZ[1]->SetStats(kFALSE);
703   gPad->SetLogz();
704   gPad->SetLeftMargin(bottom);
705   gPad->SetTopMargin(0.1);
706   gPad->SetBottomMargin(bottom);
707   gPad->SetRightMargin(0.2);
708   fhDcaXYvsDcaZ[1]->Draw("COLZ");
709
710   canvas2->SaveAs(Form("%s.eps", canvas2->GetName()));
711   canvas2->SaveAs(Form("%s.ps", canvas2->GetName()));
712
713   // -----
714
715   TCanvas* canvas3 = new TCanvas("Track_QA_Primary_3", "Track QA Primary 3", 800, 800);
716   canvas3->Divide(2, 2);
717
718   canvas3->cd(1);
719   gPad->SetRightMargin(right);
720   gPad->SetLeftMargin(left);
721   gPad->SetTopMargin(top);
722   gPad->SetBottomMargin(bottom);
723   fhQA[kDcaXYnorm][0]->SetStats(kFALSE);
724   fhQA[kDcaXYnorm][0]->Draw();
725   fhQA[kDcaXYnorm][1]->Draw("same");
726
727   canvas3->cd(2);
728   gPad->SetRightMargin(right);
729   gPad->SetLeftMargin(left);
730   gPad->SetTopMargin(top);
731   gPad->SetBottomMargin(bottom);
732   fhQA[kDcaZnorm][0]->SetStats(kFALSE);
733   fhQA[kDcaZnorm][0]->Draw();
734   fhQA[kDcaZnorm][1]->Draw("same");
735
736   canvas3->cd(3);
737 //   fhDXYvsDZ[0]->SetStats(kFALSE);
738   gPad->SetLogz();
739   gPad->SetLeftMargin(bottom);
740   gPad->SetTopMargin(0.1);
741   gPad->SetBottomMargin(bottom);
742   gPad->SetRightMargin(0.2);
743   fhDcaXYvsDcaZnorm[0]->Draw("COLZ");
744
745   canvas3->cd(4);
746 //   fhDXYvsDZ[1]->SetStats(kFALSE);
747   gPad->SetLogz();
748   gPad->SetLeftMargin(bottom);
749   gPad->SetTopMargin(0.1);
750   gPad->SetBottomMargin(bottom);
751   gPad->SetRightMargin(0.2);
752   fhDcaXYvsDcaZnorm[1]->Draw("COLZ");
753
754   canvas3->SaveAs(Form("%s.eps", canvas3->GetName()));
755   canvas3->SaveAs(Form("%s.ps", canvas3->GetName()));
756
757   // -----
758
759   TCanvas* canvas4 = new TCanvas("Track_QA_Primary_4", "Track QA Primary 4", 1200, 500);
760   canvas4->Divide(3, 1);
761
762   canvas4->cd(1);
763   gPad->SetRightMargin(right);
764   gPad->SetLeftMargin(left);
765   gPad->SetTopMargin(top);
766   gPad->SetBottomMargin(bottom);
767   fhQA[kCutNSigmaToVertex][0]->SetStats(kFALSE);
768   fhQA[kCutNSigmaToVertex][0]->Draw();
769   fhQA[kCutNSigmaToVertex][1]->Draw("same");
770
771   canvas4->cd(2);
772   gPad->SetRightMargin(right);
773   gPad->SetLeftMargin(left);
774   gPad->SetTopMargin(top);
775   gPad->SetBottomMargin(bottom);
776   fhQA[kCutRequireSigmaToVertex][0]->SetStats(kFALSE);
777   fhQA[kCutRequireSigmaToVertex][0]->Draw();
778   fhQA[kCutRequireSigmaToVertex][1]->Draw("same");
779
780   canvas4->cd(3);
781   gPad->SetRightMargin(right);
782   gPad->SetLeftMargin(left);
783   gPad->SetTopMargin(top);
784   gPad->SetBottomMargin(bottom);
785   fhQA[kCutAcceptKinkDaughters][0]->SetStats(kFALSE);
786   fhQA[kCutAcceptKinkDaughters][0]->Draw();
787   fhQA[kCutAcceptKinkDaughters][1]->Draw("same");
788
789   canvas4->SaveAs(Form("%s.eps", canvas4->GetName()));
790   canvas4->SaveAs(Form("%s.ps", canvas4->GetName()));
791 }
792 //__________________________________________________________________________________
793 void AliCFTrackIsPrimaryCuts::AddQAHistograms(TList *qaList) const {
794   //
795   // saves the histograms in a TList
796   //
797   if(!fIsQAOn) return;
798
799   qaList->Add(fhCutStatistics);
800   qaList->Add(fhCutCorrelation);
801
802   for (Int_t j=0; j<kNStepQA; j++) {
803     qaList->Add(fhDcaXYvsDcaZ[j]);
804     qaList->Add(fhDcaXYvsDcaZnorm[j]);    
805     for(Int_t i=0; i<kNHist; i++)
806       qaList->Add(fhQA[i][j]);
807   }
808 }