Changing to centrality dependent corrections
[u/mrichter/AliRoot.git] / ANALYSIS / AliESDv0Cuts.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 /* $Id: AliESDv0Cuts.cxx 24534 2008-03-16 22:22:11Z fca $ */
17
18 #include "AliESDv0Cuts.h"
19
20 #include <AliESDVertex.h>
21 #include <AliESDtrack.h>
22 #include <AliESDv0.h>
23 #include <AliESD.h>
24 #include <AliESDEvent.h>
25 #include <AliLog.h>
26
27 #include <TTree.h>
28 #include <TCanvas.h>
29 #include <TDirectory.h>
30
31 //____________________________________________________________________
32 ClassImp(AliESDv0Cuts)
33
34 // Cut names
35 const Char_t* AliESDv0Cuts::fgkCutNames[kNCuts] = {
36  "dca positive to pvtx",
37  "dca negative to pvtx",
38  "#Chi^{2}",
39  "dca v0 daughters",
40  "min decay radius",
41  "max decay radius",
42  "cosine pointing angle",
43  "on-the-fly status",
44  "dca v0 to pvtx",
45  "p",
46  "p_{T}",
47  "p_{x}",
48  "p_{y}",
49  "p_{z}"
50 };
51
52 //____________________________________________________________________
53 AliESDv0Cuts::AliESDv0Cuts(const Char_t* name, const Char_t* title) : AliAnalysisCuts(name,title),
54   fCutMinDcaPosToVertex(0),
55   fCutMinDcaNegToVertex(0),
56   fCutMaxChi2(40),
57   fCutMaxDcaV0Daughters(0),
58   fCutMinRadius(0),
59   fCutMaxRadius(0),
60   fCutMinCosinePointingAngle(0),
61   fCutRequireOnFlyStatus(0),
62   fCutMaxDcaV0ToVertex(0),
63   fPMin(0),
64   fPMax(0),
65   fPtMin(0),
66   fPtMax(0),
67   fPxMin(0),
68   fPxMax(0),
69   fPyMin(0),
70   fPyMax(0),
71   fPzMin(0),
72   fPzMax(0),
73   fHistogramsOn(0),
74   fhCutStatistics(0),         
75   fhCutCorrelation(0)
76 {
77   //
78   // constructor
79   //
80
81   Init();
82
83   //##############################################################################
84   // setting default cuts
85   SetMinDcaPosToVertex();
86   SetMinDcaNegToVertex();
87   SetMaxChi2();
88   SetMaxDcaV0Daughters();
89   SetMinRadius();
90   SetMaxRadius();
91   SetMinCosinePointingAngle();
92   SetRequireOnFlyStatus();
93   SetMaxDcaV0ToVertex();
94   SetPRange();
95   SetPtRange();
96   SetPxRange();
97   SetPyRange();
98   SetPzRange();
99
100   SetHistogramsOn();
101 }
102
103 //_____________________________________________________________________________
104 AliESDv0Cuts::AliESDv0Cuts(const AliESDv0Cuts &c) : AliAnalysisCuts(c),
105   fCutMinDcaPosToVertex(0),
106   fCutMinDcaNegToVertex(0),
107   fCutMaxChi2(0),
108   fCutMaxDcaV0Daughters(0),
109   fCutMinRadius(0),
110   fCutMaxRadius(0),
111   fCutMinCosinePointingAngle(0),
112   fCutRequireOnFlyStatus(0),
113   fCutMaxDcaV0ToVertex(0),
114   fPMin(0),
115   fPMax(0),
116   fPtMin(0),
117   fPtMax(0),
118   fPxMin(0),
119   fPxMax(0),
120   fPyMin(0),
121   fPyMax(0),
122   fPzMin(0),
123   fPzMax(0),
124   fHistogramsOn(0),
125   fhCutStatistics(0),         
126   fhCutCorrelation(0)
127 {
128   //
129   // copy constructor
130   //
131
132   ((AliESDv0Cuts &) c).Copy(*this);
133 }
134
135 AliESDv0Cuts::~AliESDv0Cuts()
136 {
137   //
138   // destructor
139   //
140
141   for (Int_t i=0; i<2; i++) {
142     
143     if (fhDcaPosToVertex[i])
144       delete fhDcaPosToVertex[i];
145     if (fhDcaNegToVertex[i])
146       delete fhDcaNegToVertex[i];
147     if (fhChi2[i])
148       delete fhChi2[i];
149     if (fhDcaV0Daughters[i])
150       delete fhDcaV0Daughters[i];
151     if (fhRadius[i])
152       delete fhRadius[i];             
153     if (fhCosinePointingAngle[i])
154       delete fhCosinePointingAngle[i];
155     if (fhOnFlyStatus[i])
156     delete fhOnFlyStatus[i];
157     if (fhDcaV0ToVertex[i])
158       delete fhDcaV0ToVertex[i];
159     if (fhPt[i])
160       delete fhPt[i];
161   }
162
163   if (fhCutStatistics)
164     delete fhCutStatistics;             
165   if (fhCutCorrelation)
166     delete fhCutCorrelation;            
167 }
168
169 void AliESDv0Cuts::Init()
170 {
171   //
172   // sets everything to zero
173   //
174   fCutMinDcaPosToVertex      = 0;
175   fCutMinDcaNegToVertex      = 0;
176   fCutMaxChi2                = 0;
177   fCutMaxDcaV0Daughters      = 0;
178   fCutMinRadius              = 0;
179   fCutMaxRadius              = 0;
180   fCutMinCosinePointingAngle = 0;
181   fCutRequireOnFlyStatus     = 0;
182   fCutMaxDcaV0ToVertex       = 0;
183
184   fPMin  = 0;
185   fPMax  = 0;
186   fPtMin = 0;
187   fPtMax = 0;
188   fPxMin = 0;
189   fPxMax = 0;
190   fPyMin = 0;
191   fPyMax = 0;
192   fPzMin = 0;
193   fPzMax = 0;
194
195   fHistogramsOn = kFALSE;
196
197   for (Int_t i=0; i<2; ++i)
198   {
199     fhDcaPosToVertex[i]      = 0;
200     fhDcaNegToVertex[i]      = 0;
201     fhChi2[i]                = 0;
202     fhDcaV0Daughters[i]      = 0;
203     fhRadius[i]              = 0;
204     fhCosinePointingAngle[i] = 0;
205     fhOnFlyStatus[i]         = 0;
206     fhDcaV0ToVertex[i]       = 0;
207     
208     fhPt[i]                  = 0;
209   }
210   fhCutStatistics = 0;
211   fhCutCorrelation = 0;
212 }
213
214 //_____________________________________________________________________________
215 AliESDv0Cuts &AliESDv0Cuts::operator=(const AliESDv0Cuts &c)
216 {
217   //
218   // Assignment operator
219   //
220
221   if (this != &c) ((AliESDv0Cuts &) c).Copy(*this);
222   return *this;
223 }
224
225 //_____________________________________________________________________________
226 void AliESDv0Cuts::Copy(TObject &c) const
227 {
228   //
229   // Copy function
230   //
231
232   AliESDv0Cuts& target = (AliESDv0Cuts &) c;
233
234   target.Init();
235
236   target.fCutMinDcaPosToVertex      = fCutMinDcaPosToVertex;
237   target.fCutMinDcaNegToVertex      = fCutMinDcaNegToVertex;
238   target.fCutMaxChi2                = fCutMaxChi2;
239   target.fCutMaxDcaV0Daughters      = fCutMaxDcaV0Daughters;
240   target.fCutMinRadius              = fCutMinRadius;
241   target.fCutMaxRadius              = fCutMaxRadius;
242   target.fCutMinCosinePointingAngle = fCutMinCosinePointingAngle;
243   target.fCutRequireOnFlyStatus     = fCutRequireOnFlyStatus;
244   target.fCutMaxDcaV0ToVertex       = fCutMaxDcaV0ToVertex;
245
246   target.fPMin  = fPMin;
247   target.fPMax  = fPMax;
248   target.fPtMin = fPtMin;
249   target.fPtMax = fPtMax;
250   target.fPxMin = fPxMin;
251   target.fPxMax = fPxMax;
252   target.fPyMin = fPyMin;
253   target.fPyMax = fPyMax;
254   target.fPzMin = fPzMin;
255   target.fPzMax = fPzMax;
256
257   target.fHistogramsOn = fHistogramsOn;
258
259   for (Int_t i=0; i<2; ++i)
260   {
261     if (fhDcaPosToVertex[i]) target.fhDcaPosToVertex[i] = (TH1F*) fhDcaPosToVertex[i]->Clone();
262     if (fhDcaNegToVertex[i]) target.fhDcaNegToVertex[i] = (TH1F*) fhDcaNegToVertex[i]->Clone();
263     if (fhChi2[i]) target.fhChi2[i] = (TH1F*) fhChi2[i]->Clone();
264     if (fhDcaV0Daughters[i]) target.fhDcaV0Daughters[i] = (TH1F*) fhDcaV0Daughters[i]->Clone();
265     if (fhRadius[i]) target.fhRadius[i] = (TH1F*) fhRadius[i]->Clone();
266     if (fhCosinePointingAngle[i]) target.fhCosinePointingAngle[i] = (TH1F*) fhCosinePointingAngle[i]->Clone();
267     if (fhOnFlyStatus[i]) target.fhOnFlyStatus[i] = (TH1F*) fhOnFlyStatus[i]->Clone();
268     if (fhDcaV0ToVertex[i]) target.fhDcaV0ToVertex[i] = (TH1F*) fhDcaV0ToVertex[i]->Clone();
269     
270     if (fhPt[i]) target.fhPt[i] = (TH1F*) fhPt[i]->Clone();
271   }
272   if (fhCutStatistics) target.fhCutStatistics = (TH1F*) fhCutStatistics->Clone();
273   if (fhCutCorrelation) target.fhCutCorrelation = (TH2F*) fhCutCorrelation->Clone();
274
275   TNamed::Copy(c);
276 }
277
278 //_____________________________________________________________________________
279 Long64_t AliESDv0Cuts::Merge(TCollection* list) {
280   // Merge a list of AliESDv0Cuts objects with this (needed for PROOF)
281   // Returns the number of merged objects (including this)
282
283   if (!list)
284     return 0;
285   
286   if (list->IsEmpty())
287     return 1;
288
289   if (!fHistogramsOn)
290     return 0;
291
292   TIterator* iter = list->MakeIterator();
293   TObject* obj;
294
295
296   // collection of measured and generated histograms
297   Int_t count = 0;
298   while ((obj = iter->Next())) {
299
300     AliESDv0Cuts* entry = dynamic_cast<AliESDv0Cuts*>(obj);
301     if (entry == 0)
302       continue;
303
304     if (!entry->fHistogramsOn)
305       continue;
306     
307     for (Int_t i=0; i<2; i++) {
308       
309       fhDcaPosToVertex[i]     ->Add(entry->fhDcaPosToVertex[i]     );
310       fhDcaNegToVertex[i]     ->Add(entry->fhDcaNegToVertex[i]     );
311       fhChi2[i]               ->Add(entry->fhChi2[i]               );
312       fhDcaV0Daughters[i]     ->Add(entry->fhDcaV0Daughters[i]     );
313       fhRadius[i]             ->Add(entry->fhRadius[i]             );
314       fhCosinePointingAngle[i]->Add(entry->fhCosinePointingAngle[i]);
315       fhOnFlyStatus[i]        ->Add(entry->fhOnFlyStatus[i]        );
316       fhDcaV0ToVertex[i]      ->Add(entry->fhDcaV0ToVertex[i]      );
317
318       fhPt[i]                ->Add(entry->fhPt[i]);
319     }      
320     fhCutStatistics  ->Add(entry->fhCutStatistics);        
321     fhCutCorrelation ->Add(entry->fhCutCorrelation);      
322
323     count++;
324   }
325
326   return count+1;
327 }
328
329 void AliESDv0Cuts::EnableNeededBranches(TTree* tree)
330 {
331   // enables the branches needed by AcceptV0, for a list see comment of AcceptV0
332
333   tree->SetBranchStatus("fV0s.fDcaV0Daughters", 1);
334   tree->SetBranchStatus("fV0s.fChi2V0", 1);
335   tree->SetBranchStatus("fV0s.fPos*", 1);
336   tree->SetBranchStatus("fV0s.fNmom*", 1);
337   tree->SetBranchStatus("fV0s.fPmom*", 1);
338   tree->SetBranchStatus("fV0s.fRr", 1);
339   tree->SetBranchStatus("fV0s.fPointAngle*", 1);
340   tree->SetBranchStatus("fV0s.fOnFlyStatus", 1);
341 }
342
343 //____________________________________________________________________
344 Bool_t
345 AliESDv0Cuts::IsSelected(TList* listObj) {
346 // Selection cuts
347   if(listObj->GetSize()!=4) return kFALSE;
348   AliESDv0           *esdV0     = (AliESDv0*)listObj->At(0);
349   AliESDtrack        *trackPos  = (AliESDtrack*)listObj->At(1);
350   AliESDtrack        *trackNeg  = (AliESDtrack*)listObj->At(2);
351   const AliESDVertex *esdVertex = (AliESDVertex*)listObj->At(3);
352   return AcceptV0(esdV0,trackPos,trackNeg,esdVertex);
353 }
354
355 //____________________________________________________________________
356 Bool_t
357 AliESDv0Cuts::AcceptV0(AliESDv0* const esdV0, AliESDtrack* const trackPos, AliESDtrack* const trackNeg, const AliESDVertex* esdVertex) {
358   // 
359   // figure out if the v0s survives all the v0 cuts defined
360   //
361   // the different quality parameter and kinematic values are first
362   // retrieved from the v0. then it is found out what cuts the
363   // v0 did not survive and finally the cuts are imposed.
364
365   // this function needs the following branches... but this is not enough
366   // fV0s.fDcaV0Daughters
367   // fV0s.fChi2V0
368   // fV0s.fPos*
369   // fV0s.fNmom*
370   // fV0s.fPmom*
371   // fV0s.fRr
372   // fV0s.fPointAngle
373   // fV0s.fOnFlyStatus
374
375   Float_t  dcaPosToVertex = 0, dcaNegToVertex = 0;
376   Float_t  tdcaPosToVertex[2]={999,999};
377   Float_t  tdcaNegToVertex[2]={999,999};
378
379   if (trackPos) trackPos->GetImpactParameters(tdcaPosToVertex[0],tdcaPosToVertex[1]);
380   if (trackNeg) trackNeg->GetImpactParameters(tdcaNegToVertex[0],tdcaNegToVertex[1]);
381
382   dcaPosToVertex = TMath::Sqrt(tdcaPosToVertex[0]*tdcaPosToVertex[0]+tdcaPosToVertex[1]*tdcaPosToVertex[1]);
383   dcaNegToVertex = TMath::Sqrt(tdcaNegToVertex[0]*tdcaNegToVertex[0]+tdcaNegToVertex[1]*tdcaNegToVertex[1]);
384
385   UInt_t  status = esdV0->GetOnFlyStatus();
386   Float_t chi2  = esdV0->GetChi2V0();
387
388   Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
389
390   Double_t vtxPosition[3]; esdVertex->GetXYZ(vtxPosition);
391   Double_t dcaV0ToVertex  = esdV0->GetD(vtxPosition[0],vtxPosition[1],vtxPosition[2]);
392
393   Double_t v0Position[3];
394   esdV0->GetXYZ(v0Position[0],v0Position[1],v0Position[2]);
395   Double_t radius = TMath::Sqrt(TMath::Power(v0Position[0],2) + TMath::Power(v0Position[1],2));
396   Double_t v0CosinePointingAngle = esdV0->GetV0CosineOfPointingAngle();
397
398   // getting the kinematic variables of the v0
399   Double_t p[3];
400   esdV0->GetPxPyPz(p[0],p[1],p[2]);
401   Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
402   Float_t pt       = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
403
404   //########################################################################
405   // cut the v0?
406   
407   Bool_t cuts[kNCuts];
408   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
409   
410   // v0 quality cuts
411   if (dcaPosToVertex < fCutMinDcaPosToVertex)
412     cuts[0]=kTRUE;
413   if (dcaNegToVertex < fCutMinDcaNegToVertex) 
414     cuts[1]=kTRUE;
415   if (chi2 > fCutMaxChi2) 
416     cuts[2]=kTRUE; 
417   if (dcaV0Daughters > fCutMaxDcaV0Daughters) 
418     cuts[3]=kTRUE;  
419   if (radius  < fCutMinRadius) 
420     cuts[4]=kTRUE;  
421   if (radius  > fCutMaxRadius) 
422     cuts[5]=kTRUE;  
423   if (v0CosinePointingAngle < fCutMinCosinePointingAngle)
424     cuts[6]=kTRUE;  
425   if (fCutRequireOnFlyStatus && !status)
426     cuts[7]=kTRUE;  
427   if (dcaV0ToVertex > fCutMaxDcaV0ToVertex)
428     cuts[8] = kTRUE;
429
430   // v0 kinematics cut
431   if((momentum < fPMin) || (momentum > fPMax)) 
432     cuts[9]=kTRUE;
433   if((pt < fPtMin) || (pt > fPtMax)) 
434     cuts[10] = kTRUE;
435   if((p[0] < fPxMin) || (p[0] > fPxMax)) 
436     cuts[11] = kTRUE;
437   if((p[1] < fPyMin) || (p[1] > fPyMax)) 
438     cuts[12] = kTRUE;
439   if((p[2] < fPzMin) || (p[2] > fPzMax))
440     cuts[13] = kTRUE;
441
442   Bool_t cut=kFALSE;
443   for (Int_t i=0; i<kNCuts; i++) 
444     if (cuts[i]) cut = kTRUE;
445   
446   //########################################################################
447   // filling histograms
448   if (fHistogramsOn) {
449     fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n v0s")));
450     
451     if (cut)
452       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut v0s")));
453     
454     for (Int_t i=0; i<kNCuts; i++) {
455       if (cuts[i])
456         fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
457       
458       for (Int_t j=i; j<kNCuts; j++) {
459         if (cuts[i] && cuts[j]) {
460           Float_t x = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
461           Float_t y = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
462           fhCutCorrelation->Fill(x,y);
463         }
464       }
465     }
466     
467     fhDcaPosToVertex[0]->Fill(dcaPosToVertex);
468     fhDcaNegToVertex[0]->Fill(dcaNegToVertex);
469     fhChi2[0]->Fill(chi2);
470     fhDcaV0Daughters[0]->Fill(dcaV0Daughters);
471     fhRadius[0]->Fill(radius);
472     fhCosinePointingAngle[0]->Fill(v0CosinePointingAngle);
473     fhOnFlyStatus[0]->Fill(status);
474     fhDcaV0ToVertex[0]->Fill(dcaV0ToVertex);
475     
476     fhPt[0]->Fill(pt);
477   }
478
479   //########################################################################
480   // cut the v0!
481   if (cut) return kFALSE;
482
483   //########################################################################
484   // filling histograms after cut
485   if (fHistogramsOn) {
486     fhDcaPosToVertex[1]->Fill(dcaPosToVertex);
487     fhDcaNegToVertex[1]->Fill(dcaNegToVertex);
488     fhChi2[1]->Fill(chi2);
489     fhDcaV0Daughters[1]->Fill(dcaV0Daughters);
490     fhRadius[1]->Fill(radius);
491     fhCosinePointingAngle[1]->Fill(v0CosinePointingAngle);
492     fhOnFlyStatus[1]->Fill(status);
493     fhDcaV0ToVertex[1]->Fill(dcaV0ToVertex);
494     
495     fhPt[1]->Fill(pt);
496   }
497
498   return kTRUE;
499 }
500
501 //____________________________________________________________________
502 TObjArray* AliESDv0Cuts::GetAcceptedV0s(const AliESD* esd)
503 {
504   //
505   // returns an array of all v0s that pass the cuts
506   //
507
508   TObjArray* acceptedV0s = new TObjArray();
509   //  const AliESDVertex *spdVertex = esd->GetVertex();
510   const AliESDVertex *primaryVertex = esd->GetPrimaryVertex();
511   Int_t    lIndexTrackPos       = 0, lIndexTrackNeg       = 0;
512
513   // loop over esd v0s
514   for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) {
515     AliESDv0* v0 = esd->GetV0(iV0);
516
517     lIndexTrackPos = TMath::Abs(v0->GetPindex());
518     lIndexTrackNeg = TMath::Abs(v0->GetNindex());
519     AliESDtrack *trackPos = esd->GetTrack(lIndexTrackPos);
520     AliESDtrack *trackNeg = esd->GetTrack(lIndexTrackNeg);
521
522     if (AcceptV0(v0,trackPos,trackNeg,primaryVertex))
523       acceptedV0s->Add(v0);
524   }
525
526   return acceptedV0s;
527 }
528
529 //____________________________________________________________________
530 Int_t AliESDv0Cuts::CountAcceptedV0s(const AliESD* esd)
531 {
532   //
533   // returns an the number of v0s that pass the cuts
534   //
535
536   Int_t count = 0;
537   //  const AliESDVertex *spdVertex = esd->GetVertex();
538   const AliESDVertex *primaryVertex = esd->GetPrimaryVertex();
539   Int_t    lIndexTrackPos       = 0, lIndexTrackNeg       = 0;
540
541   // loop over esd v0s
542   for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) {
543     AliESDv0* v0 = esd->GetV0(iV0);
544
545     lIndexTrackPos = TMath::Abs(v0->GetPindex());
546     lIndexTrackNeg = TMath::Abs(v0->GetNindex());
547     AliESDtrack *trackPos = esd->GetTrack(lIndexTrackPos);
548     AliESDtrack *trackNeg = esd->GetTrack(lIndexTrackNeg);
549
550     if (AcceptV0(v0,trackPos,trackNeg,primaryVertex))
551       count++;
552   }
553
554   return count;
555 }
556
557 //____________________________________________________________________
558 TObjArray* AliESDv0Cuts::GetAcceptedV0s(const AliESDEvent* esd)
559 {
560   //
561   // returns an array of all v0s that pass the cuts
562   //
563
564   TObjArray* acceptedV0s = new TObjArray();
565   //  const AliESDVertex *spdVertex = esd->GetVertex();
566   const AliESDVertex *primaryVertex = esd->GetPrimaryVertex();
567   Int_t    lIndexTrackPos       = 0, lIndexTrackNeg       = 0;
568
569   // loop over esd v0s
570   for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) {
571     AliESDv0* v0 = esd->GetV0(iV0);
572
573     lIndexTrackPos = TMath::Abs(v0->GetPindex());
574     lIndexTrackNeg = TMath::Abs(v0->GetNindex());
575     AliESDtrack *trackPos = esd->GetTrack(lIndexTrackPos);
576     AliESDtrack *trackNeg = esd->GetTrack(lIndexTrackNeg);
577
578     if (AcceptV0(v0,trackPos,trackNeg,primaryVertex))
579       acceptedV0s->Add(v0);
580   }
581
582   return acceptedV0s;
583 }
584
585 //____________________________________________________________________
586 Int_t AliESDv0Cuts::CountAcceptedV0s(const AliESDEvent* esd)
587 {
588   //
589   // returns an the number of v0s that pass the cuts
590   //
591
592   Int_t count = 0;
593   //  const AliESDVertex *spdVertex = esd->GetVertex();
594   const AliESDVertex *primaryVertex = esd->GetPrimaryVertex();
595   Int_t    lIndexTrackPos       = 0, lIndexTrackNeg       = 0;
596
597   // loop over esd v0s
598   for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) {
599     AliESDv0* v0 = esd->GetV0(iV0);
600
601     lIndexTrackPos = TMath::Abs(v0->GetPindex());
602     lIndexTrackNeg = TMath::Abs(v0->GetNindex());
603     AliESDtrack *trackPos = esd->GetTrack(lIndexTrackPos);
604     AliESDtrack *trackNeg = esd->GetTrack(lIndexTrackNeg);
605
606     if (AcceptV0(v0,trackPos,trackNeg,primaryVertex))
607       count++;
608   }
609
610   return count;
611 }
612
613 //____________________________________________________________________
614  void AliESDv0Cuts::DefineHistograms(Int_t color) {
615    // 
616    // diagnostics histograms are defined
617    // 
618
619    fHistogramsOn=kTRUE;
620    
621    Bool_t oldStatus = TH1::AddDirectoryStatus();
622    TH1::AddDirectory(kFALSE);
623    
624    //###################################################################################
625    // defining histograms
626
627    fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
628
629    fhCutStatistics->GetXaxis()->SetBinLabel(1,"n v0s");
630    fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut v0s");
631
632    fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
633   
634    for (Int_t i=0; i<kNCuts; i++) {
635      fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
636      fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
637      fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
638    } 
639
640   fhCutStatistics  ->SetLineColor(color);
641   fhCutCorrelation ->SetLineColor(color);
642   fhCutStatistics  ->SetLineWidth(2);
643   fhCutCorrelation ->SetLineWidth(2);
644
645   Char_t str[256];
646   for (Int_t i=0; i<2; i++) {
647     if (i==0) snprintf(str,256, " ");
648     else snprintf(str,256, "_cut");
649
650     fhDcaPosToVertex[i]      = new TH1F(Form("dcaPosToVertex%s",str),"",120,0,3);
651     fhDcaNegToVertex[i]      = new TH1F(Form("dcaNegToVertex%s",str),"",120,0,3);
652     fhChi2[i]                = new TH1F(Form("chi2%s",str),"",50,0,50);
653     fhDcaV0Daughters[i]      = new  TH1F(Form("dcaV0Daughters%s",str),"",200,0,5);
654     fhRadius[i]              = new  TH1F(Form("decayRadius%s",str),"",300,0,150);
655     fhCosinePointingAngle[i] = new  TH1F(Form("cosinePointingAngle%s",str),"",100,-1,1);
656     fhOnFlyStatus[i]         = new  TH1F(Form("onflystatus%s",str),"",5,0,5);
657     fhDcaV0ToVertex[i]       = new  TH1F(Form("dcaV0ToVertex%s",str),"",100,0,5);
658
659     fhPt[i]                  = new TH1F(Form("pt%s",str)     ,"p_{T} distribution;p_{T} (GeV/c)",500,0.0,100.0);
660     
661     fhDcaPosToVertex[i]->SetTitle("Dca of positive daughter to parent vertex");
662     fhDcaNegToVertex[i]->SetTitle("Dca of negative daughter to parent vertex");
663     fhChi2[i]->SetTitle("#Chi^{2} for v0");
664     fhDcaV0Daughters[i]->SetTitle("Dca between daughter tracks");
665     fhRadius[i]->SetTitle("Decay radius of the v0");
666     fhCosinePointingAngle[i]->SetTitle("Cosine of the Pointing Angle");
667     fhOnFlyStatus[i]->SetTitle("On-the-Fly Status");
668     fhDcaV0ToVertex[i]->SetTitle("Dca of v0 to parent vertex");
669
670     fhDcaPosToVertex[i]->SetLineColor(color);      fhDcaPosToVertex[i]->SetLineWidth(2);
671     fhDcaNegToVertex[i]->SetLineColor(color);      fhDcaNegToVertex[i]->SetLineWidth(2);
672     fhChi2[i]->SetLineColor(color);                fhChi2[i]->SetLineWidth(2);
673     fhDcaV0Daughters[i]->SetLineColor(color);      fhDcaV0Daughters[i]->SetLineWidth(2);
674     fhRadius[i]->SetLineColor(color);              fhRadius[i]->SetLineWidth(2);
675     fhCosinePointingAngle[i]->SetLineColor(color); fhCosinePointingAngle[i]->SetLineWidth(2);
676     fhOnFlyStatus[i]->SetLineColor(color);         fhOnFlyStatus[i]->SetLineWidth(2);
677     fhDcaV0ToVertex[i]->SetLineColor(color);       fhDcaV0ToVertex[i]->SetLineWidth(2); 
678   }
679
680   TH1::AddDirectory(oldStatus);
681 }
682
683 //____________________________________________________________________
684 Bool_t AliESDv0Cuts::LoadHistograms(const Char_t* dir)
685 {
686   //
687   // loads the histograms from a file
688   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
689   //
690
691   if (!dir)
692     dir = GetName();
693
694   if (!gDirectory->cd(dir))
695     return kFALSE;
696
697   fhCutStatistics = dynamic_cast<TH1F*> (gDirectory->Get("cut_statistics"));
698   fhCutCorrelation = dynamic_cast<TH2F*> (gDirectory->Get("cut_correlation"));
699
700   Char_t str[5];
701   for (Int_t i=0; i<2; i++) {
702     if (i==0)
703     {
704       gDirectory->cd("before_cuts");
705       str[0] = 0;
706     }
707     else
708     {
709       gDirectory->cd("after_cuts");
710       snprintf(str,5, "_cut");
711     }
712
713     fhDcaPosToVertex[i]      = dynamic_cast<TH1F*> (gDirectory->Get(Form("dcaPosToVertex%s",str)     ));
714     fhDcaNegToVertex[i]      = dynamic_cast<TH1F*> (gDirectory->Get(Form("dcaNegToVertex%s",str)     ));
715     fhChi2[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("chi2%s",str)));
716     fhDcaV0Daughters[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("dcaV0Daughters%s",str)));
717     fhRadius[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("decayRadius%s",str)));
718     fhCosinePointingAngle[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("cosinepointingangle%s",str)));
719     fhOnFlyStatus[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("onflystatus%s",str)));
720     fhDcaV0ToVertex[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("dcaV0ToVertex%s",str)));
721
722     fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("pt%s",str)));
723
724     gDirectory->cd("../");
725   }
726
727   gDirectory->cd("..");
728
729   return kTRUE;
730 }
731
732 //____________________________________________________________________
733 void AliESDv0Cuts::SaveHistograms(const Char_t* dir) {
734   //
735   // saves the histograms in a directory (dir)
736   //
737
738   if (!fHistogramsOn) {
739     AliDebug(0, "Histograms not on - cannot save histograms!!!");
740     return;
741   }
742
743   if (!dir)
744     dir = GetName();
745
746   gDirectory->mkdir(dir);
747   gDirectory->cd(dir);
748
749   gDirectory->mkdir("before_cuts");
750   gDirectory->mkdir("after_cuts");
751  
752   fhCutStatistics->Write();
753   fhCutCorrelation->Write();
754
755   for (Int_t i=0; i<2; i++) {
756     if (i==0)
757       gDirectory->cd("before_cuts");
758     else
759       gDirectory->cd("after_cuts");
760
761     fhDcaPosToVertex[i]      ->Write();
762     fhDcaNegToVertex[i]      ->Write();
763     fhChi2[i]                ->Write();
764     fhDcaV0Daughters[i]      ->Write();
765     fhRadius[i]              ->Write();
766     fhCosinePointingAngle[i] ->Write();
767     fhOnFlyStatus[i]         ->Write();
768     fhDcaV0ToVertex[i]       ->Write();
769
770     fhPt[i]                  ->Write();
771     
772     gDirectory->cd("../");
773   }
774
775   gDirectory->cd("../");
776 }
777
778 //____________________________________________________________________
779 void AliESDv0Cuts::DrawHistograms()
780 {
781   // draws some histograms
782
783   TCanvas* canvas1 = new TCanvas(Form("%s_1", GetName()), "V0 Quality Results1", 800, 800);
784   canvas1->Divide(2, 2);
785
786   canvas1->cd(1);
787   fhDcaPosToVertex[0]->SetStats(kFALSE);
788   fhDcaPosToVertex[0]->Draw();
789
790   canvas1->cd(2);
791   fhChi2[0]->SetStats(kFALSE);
792   fhChi2[0]->Draw();
793
794   canvas1->cd(3);
795   fhDcaV0ToVertex[0]->SetStats(kFALSE);
796   fhDcaV0ToVertex[0]->Draw();
797
798   canvas1->SaveAs(Form("%s_%s.gif", GetName(), canvas1->GetName()));
799
800   TCanvas* canvas2 = new TCanvas(Form("%s_2", GetName()), "V0 Quality Results2", 1200, 800);
801   canvas2->Divide(3, 2);
802
803   canvas2->cd(1);
804   fhDcaV0Daughters[0]->SetStats(kFALSE);
805   gPad->SetLogy();
806   fhDcaV0Daughters[0]->Draw();
807
808   canvas2->cd(2);
809   fhRadius[0]->SetStats(kFALSE);
810   gPad->SetLogy();
811   fhRadius[0]->Draw();
812
813
814   canvas2->cd(4);
815   fhCosinePointingAngle[0]->SetStats(kFALSE);
816   gPad->SetLogy();
817   fhCosinePointingAngle[0]->Draw();
818
819   canvas2->cd(5);
820   fhOnFlyStatus[0]->SetStats(kFALSE);
821   gPad->SetLogy();
822   fhOnFlyStatus[0]->Draw();
823
824   canvas2->SaveAs(Form("%s_%s.gif", GetName(), canvas2->GetName()));
825
826   TCanvas* canvas3 = new TCanvas(Form("%s_4", GetName()), "V0 Quality Results3", 800, 500);
827   canvas3->Divide(2, 1);
828
829   canvas3->cd(1);
830   fhCutStatistics->SetStats(kFALSE);
831   fhCutStatistics->LabelsOption("v");
832   gPad->SetBottomMargin(0.3);
833   fhCutStatistics->Draw();
834
835   canvas3->cd(2);
836   fhCutCorrelation->SetStats(kFALSE);
837   fhCutCorrelation->LabelsOption("v");
838   gPad->SetBottomMargin(0.3);
839   gPad->SetLeftMargin(0.3);
840   fhCutCorrelation->Draw("COLZ");
841
842   canvas3->SaveAs(Form("%s_%s.gif", GetName(), canvas3->GetName()));
843 }
844