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