added error messages when task does not post data in UserCreateOutputObjects for...
[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   if(listObj->GetSize()!=4) return kFALSE;
347   AliESDv0           *esdV0     = (AliESDv0*)listObj->At(0);
348   AliESDtrack        *trackPos  = (AliESDtrack*)listObj->At(1);
349   AliESDtrack        *trackNeg  = (AliESDtrack*)listObj->At(2);
350   const AliESDVertex *esdVertex = (AliESDVertex*)listObj->At(3);
351   return AcceptV0(esdV0,trackPos,trackNeg,esdVertex);
352 }
353
354 //____________________________________________________________________
355 Bool_t
356 AliESDv0Cuts::AcceptV0(AliESDv0* esdV0, AliESDtrack* trackPos, AliESDtrack* trackNeg, const AliESDVertex* esdVertex) {
357   // 
358   // figure out if the v0s survives all the v0 cuts defined
359   //
360   // the different quality parameter and kinematic values are first
361   // retrieved from the v0. then it is found out what cuts the
362   // v0 did not survive and finally the cuts are imposed.
363
364   // this function needs the following branches... but this is not enough
365   // fV0s.fDcaV0Daughters
366   // fV0s.fChi2V0
367   // fV0s.fPos*
368   // fV0s.fNmom*
369   // fV0s.fPmom*
370   // fV0s.fRr
371   // fV0s.fPointAngle
372   // fV0s.fOnFlyStatus
373
374   Float_t  dcaPosToVertex = 0, dcaNegToVertex = 0;
375   Float_t  tdcaPosToVertex[2]={999,999};
376   Float_t  tdcaNegToVertex[2]={999,999};
377
378   if (trackPos) trackPos->GetImpactParameters(tdcaPosToVertex[0],tdcaPosToVertex[1]);
379   if (trackNeg) trackNeg->GetImpactParameters(tdcaNegToVertex[0],tdcaNegToVertex[1]);
380
381   dcaPosToVertex = TMath::Sqrt(tdcaPosToVertex[0]*tdcaPosToVertex[0]+tdcaPosToVertex[1]*tdcaPosToVertex[1]);
382   dcaNegToVertex = TMath::Sqrt(tdcaNegToVertex[0]*tdcaNegToVertex[0]+tdcaNegToVertex[1]*tdcaNegToVertex[1]);
383
384   UInt_t  status = esdV0->GetOnFlyStatus();
385   Float_t chi2  = esdV0->GetChi2V0();
386
387   Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
388
389   Double_t vtxPosition[3]; esdVertex->GetXYZ(vtxPosition);
390   Double_t dcaV0ToVertex  = esdV0->GetD(vtxPosition[0],vtxPosition[1],vtxPosition[2]);
391
392   Double_t v0Position[3];
393   esdV0->GetXYZ(v0Position[0],v0Position[1],v0Position[2]);
394   Double_t radius = TMath::Sqrt(TMath::Power(v0Position[0],2) + TMath::Power(v0Position[1],2));
395   Double_t v0CosinePointingAngle = esdV0->GetV0CosineOfPointingAngle();
396
397   // getting the kinematic variables of the v0
398   Double_t p[3];
399   esdV0->GetPxPyPz(p[0],p[1],p[2]);
400   Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
401   Float_t pt       = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
402
403   //########################################################################
404   // cut the v0?
405   
406   Bool_t cuts[kNCuts];
407   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
408   
409   // v0 quality cuts
410   if (dcaPosToVertex < fCutMinDcaPosToVertex)
411     cuts[0]=kTRUE;
412   if (dcaNegToVertex < fCutMinDcaNegToVertex) 
413     cuts[1]=kTRUE;
414   if (chi2 > fCutMaxChi2) 
415     cuts[2]=kTRUE; 
416   if (dcaV0Daughters > fCutMaxDcaV0Daughters) 
417     cuts[3]=kTRUE;  
418   if (radius  < fCutMinRadius) 
419     cuts[4]=kTRUE;  
420   if (radius  > fCutMaxRadius) 
421     cuts[5]=kTRUE;  
422   if (v0CosinePointingAngle < fCutMinCosinePointingAngle)
423     cuts[6]=kTRUE;  
424   if (fCutRequireOnFlyStatus && !status)
425     cuts[7]=kTRUE;  
426   if (dcaV0ToVertex > fCutMaxDcaV0ToVertex)
427     cuts[8] = kTRUE;
428
429   // v0 kinematics cut
430   if((momentum < fPMin) || (momentum > fPMax)) 
431     cuts[9]=kTRUE;
432   if((pt < fPtMin) || (pt > fPtMax)) 
433     cuts[10] = kTRUE;
434   if((p[0] < fPxMin) || (p[0] > fPxMax)) 
435     cuts[11] = kTRUE;
436   if((p[1] < fPyMin) || (p[1] > fPyMax)) 
437     cuts[12] = kTRUE;
438   if((p[2] < fPzMin) || (p[2] > fPzMax))
439     cuts[13] = kTRUE;
440
441   Bool_t cut=kFALSE;
442   for (Int_t i=0; i<kNCuts; i++) 
443     if (cuts[i]) cut = kTRUE;
444   
445   //########################################################################
446   // filling histograms
447   if (fHistogramsOn) {
448     fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n v0s")));
449     
450     if (cut)
451       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut v0s")));
452     
453     for (Int_t i=0; i<kNCuts; i++) {
454       if (cuts[i])
455         fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
456       
457       for (Int_t j=i; j<kNCuts; j++) {
458         if (cuts[i] && cuts[j]) {
459           Float_t x = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
460           Float_t y = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
461           fhCutCorrelation->Fill(x,y);
462         }
463       }
464     }
465     
466     fhDcaPosToVertex[0]->Fill(dcaPosToVertex);
467     fhDcaNegToVertex[0]->Fill(dcaNegToVertex);
468     fhChi2[0]->Fill(chi2);
469     fhDcaV0Daughters[0]->Fill(dcaV0Daughters);
470     fhRadius[0]->Fill(radius);
471     fhCosinePointingAngle[0]->Fill(v0CosinePointingAngle);
472     fhOnFlyStatus[0]->Fill(status);
473     fhDcaV0ToVertex[0]->Fill(dcaV0ToVertex);
474     
475     fhPt[0]->Fill(pt);
476   }
477
478   //########################################################################
479   // cut the v0!
480   if (cut) return kFALSE;
481
482   //########################################################################
483   // filling histograms after cut
484   if (fHistogramsOn) {
485     fhDcaPosToVertex[1]->Fill(dcaPosToVertex);
486     fhDcaNegToVertex[1]->Fill(dcaNegToVertex);
487     fhChi2[1]->Fill(chi2);
488     fhDcaV0Daughters[1]->Fill(dcaV0Daughters);
489     fhRadius[1]->Fill(radius);
490     fhCosinePointingAngle[1]->Fill(v0CosinePointingAngle);
491     fhOnFlyStatus[1]->Fill(status);
492     fhDcaV0ToVertex[1]->Fill(dcaV0ToVertex);
493     
494     fhPt[1]->Fill(pt);
495   }
496
497   return kTRUE;
498 }
499
500 //____________________________________________________________________
501 TObjArray* AliESDv0Cuts::GetAcceptedV0s(AliESD* esd)
502 {
503   //
504   // returns an array of all v0s that pass the cuts
505   //
506
507   TObjArray* acceptedV0s = new TObjArray();
508   //  const AliESDVertex *spdVertex = esd->GetVertex();
509   const AliESDVertex *primaryVertex = esd->GetPrimaryVertex();
510   Int_t    lIndexTrackPos       = 0, lIndexTrackNeg       = 0;
511
512   // loop over esd v0s
513   for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) {
514     AliESDv0* v0 = esd->GetV0(iV0);
515
516     lIndexTrackPos = TMath::Abs(v0->GetPindex());
517     lIndexTrackNeg = TMath::Abs(v0->GetNindex());
518     AliESDtrack *trackPos = esd->GetTrack(lIndexTrackPos);
519     AliESDtrack *trackNeg = esd->GetTrack(lIndexTrackNeg);
520
521     if (AcceptV0(v0,trackPos,trackNeg,primaryVertex))
522       acceptedV0s->Add(v0);
523   }
524
525   return acceptedV0s;
526 }
527
528 //____________________________________________________________________
529 Int_t AliESDv0Cuts::CountAcceptedV0s(AliESD* esd)
530 {
531   //
532   // returns an the number of v0s that pass the cuts
533   //
534
535   Int_t count = 0;
536   //  const AliESDVertex *spdVertex = esd->GetVertex();
537   const AliESDVertex *primaryVertex = esd->GetPrimaryVertex();
538   Int_t    lIndexTrackPos       = 0, lIndexTrackNeg       = 0;
539
540   // loop over esd v0s
541   for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) {
542     AliESDv0* v0 = esd->GetV0(iV0);
543
544     lIndexTrackPos = TMath::Abs(v0->GetPindex());
545     lIndexTrackNeg = TMath::Abs(v0->GetNindex());
546     AliESDtrack *trackPos = esd->GetTrack(lIndexTrackPos);
547     AliESDtrack *trackNeg = esd->GetTrack(lIndexTrackNeg);
548
549     if (AcceptV0(v0,trackPos,trackNeg,primaryVertex))
550       count++;
551   }
552
553   return count;
554 }
555
556 //____________________________________________________________________
557 TObjArray* AliESDv0Cuts::GetAcceptedV0s(AliESDEvent* esd)
558 {
559   //
560   // returns an array of all v0s that pass the cuts
561   //
562
563   TObjArray* acceptedV0s = new TObjArray();
564   //  const AliESDVertex *spdVertex = esd->GetVertex();
565   const AliESDVertex *primaryVertex = esd->GetPrimaryVertex();
566   Int_t    lIndexTrackPos       = 0, lIndexTrackNeg       = 0;
567
568   // loop over esd v0s
569   for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) {
570     AliESDv0* v0 = esd->GetV0(iV0);
571
572     lIndexTrackPos = TMath::Abs(v0->GetPindex());
573     lIndexTrackNeg = TMath::Abs(v0->GetNindex());
574     AliESDtrack *trackPos = esd->GetTrack(lIndexTrackPos);
575     AliESDtrack *trackNeg = esd->GetTrack(lIndexTrackNeg);
576
577     if (AcceptV0(v0,trackPos,trackNeg,primaryVertex))
578       acceptedV0s->Add(v0);
579   }
580
581   return acceptedV0s;
582 }
583
584 //____________________________________________________________________
585 Int_t AliESDv0Cuts::CountAcceptedV0s(AliESDEvent* esd)
586 {
587   //
588   // returns an the number of v0s that pass the cuts
589   //
590
591   Int_t count = 0;
592   //  const AliESDVertex *spdVertex = esd->GetVertex();
593   const AliESDVertex *primaryVertex = esd->GetPrimaryVertex();
594   Int_t    lIndexTrackPos       = 0, lIndexTrackNeg       = 0;
595
596   // loop over esd v0s
597   for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) {
598     AliESDv0* v0 = esd->GetV0(iV0);
599
600     lIndexTrackPos = TMath::Abs(v0->GetPindex());
601     lIndexTrackNeg = TMath::Abs(v0->GetNindex());
602     AliESDtrack *trackPos = esd->GetTrack(lIndexTrackPos);
603     AliESDtrack *trackNeg = esd->GetTrack(lIndexTrackNeg);
604
605     if (AcceptV0(v0,trackPos,trackNeg,primaryVertex))
606       count++;
607   }
608
609   return count;
610 }
611
612 //____________________________________________________________________
613  void AliESDv0Cuts::DefineHistograms(Int_t color) {
614    // 
615    // diagnostics histograms are defined
616    // 
617
618    fHistogramsOn=kTRUE;
619    
620    Bool_t oldStatus = TH1::AddDirectoryStatus();
621    TH1::AddDirectory(kFALSE);
622    
623    //###################################################################################
624    // defining histograms
625
626    fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
627
628    fhCutStatistics->GetXaxis()->SetBinLabel(1,"n v0s");
629    fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut v0s");
630
631    fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
632   
633    for (Int_t i=0; i<kNCuts; i++) {
634      fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
635      fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
636      fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
637    } 
638
639   fhCutStatistics  ->SetLineColor(color);
640   fhCutCorrelation ->SetLineColor(color);
641   fhCutStatistics  ->SetLineWidth(2);
642   fhCutCorrelation ->SetLineWidth(2);
643
644   Char_t str[256];
645   for (Int_t i=0; i<2; i++) {
646     if (i==0) snprintf(str,256, " ");
647     else snprintf(str,256, "_cut");
648
649     fhDcaPosToVertex[i]      = new TH1F(Form("dcaPosToVertex%s",str),"",120,0,3);
650     fhDcaNegToVertex[i]      = new TH1F(Form("dcaNegToVertex%s",str),"",120,0,3);
651     fhChi2[i]                = new TH1F(Form("chi2%s",str),"",50,0,50);
652     fhDcaV0Daughters[i]      = new  TH1F(Form("dcaV0Daughters%s",str),"",200,0,5);
653     fhRadius[i]              = new  TH1F(Form("decayRadius%s",str),"",300,0,150);
654     fhCosinePointingAngle[i] = new  TH1F(Form("cosinePointingAngle%s",str),"",100,-1,1);
655     fhOnFlyStatus[i]         = new  TH1F(Form("onflystatus%s",str),"",5,0,5);
656     fhDcaV0ToVertex[i]       = new  TH1F(Form("dcaV0ToVertex%s",str),"",100,0,5);
657
658     fhPt[i]                  = new TH1F(Form("pt%s",str)     ,"p_{T} distribution;p_{T} (GeV/c)",500,0.0,100.0);
659     
660     fhDcaPosToVertex[i]->SetTitle("Dca of positive daughter to parent vertex");
661     fhDcaNegToVertex[i]->SetTitle("Dca of negative daughter to parent vertex");
662     fhChi2[i]->SetTitle("#Chi^{2} for v0");
663     fhDcaV0Daughters[i]->SetTitle("Dca between daughter tracks");
664     fhRadius[i]->SetTitle("Decay radius of the v0");
665     fhCosinePointingAngle[i]->SetTitle("Cosine of the Pointing Angle");
666     fhOnFlyStatus[i]->SetTitle("On-the-Fly Status");
667     fhDcaV0ToVertex[i]->SetTitle("Dca of v0 to parent vertex");
668
669     fhDcaPosToVertex[i]->SetLineColor(color);      fhDcaPosToVertex[i]->SetLineWidth(2);
670     fhDcaNegToVertex[i]->SetLineColor(color);      fhDcaNegToVertex[i]->SetLineWidth(2);
671     fhChi2[i]->SetLineColor(color);                fhChi2[i]->SetLineWidth(2);
672     fhDcaV0Daughters[i]->SetLineColor(color);      fhDcaV0Daughters[i]->SetLineWidth(2);
673     fhRadius[i]->SetLineColor(color);              fhRadius[i]->SetLineWidth(2);
674     fhCosinePointingAngle[i]->SetLineColor(color); fhCosinePointingAngle[i]->SetLineWidth(2);
675     fhOnFlyStatus[i]->SetLineColor(color);         fhOnFlyStatus[i]->SetLineWidth(2);
676     fhDcaV0ToVertex[i]->SetLineColor(color);       fhDcaV0ToVertex[i]->SetLineWidth(2); 
677   }
678
679   TH1::AddDirectory(oldStatus);
680 }
681
682 //____________________________________________________________________
683 Bool_t AliESDv0Cuts::LoadHistograms(const Char_t* dir)
684 {
685   //
686   // loads the histograms from a file
687   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
688   //
689
690   if (!dir)
691     dir = GetName();
692
693   if (!gDirectory->cd(dir))
694     return kFALSE;
695
696   fhCutStatistics = dynamic_cast<TH1F*> (gDirectory->Get("cut_statistics"));
697   fhCutCorrelation = dynamic_cast<TH2F*> (gDirectory->Get("cut_correlation"));
698
699   Char_t str[5];
700   for (Int_t i=0; i<2; i++) {
701     if (i==0)
702     {
703       gDirectory->cd("before_cuts");
704       str[0] = 0;
705     }
706     else
707     {
708       gDirectory->cd("after_cuts");
709       snprintf(str,5, "_cut");
710     }
711
712     fhDcaPosToVertex[i]      = dynamic_cast<TH1F*> (gDirectory->Get(Form("dcaPosToVertex%s",str)     ));
713     fhDcaNegToVertex[i]      = dynamic_cast<TH1F*> (gDirectory->Get(Form("dcaNegToVertex%s",str)     ));
714     fhChi2[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("chi2%s",str)));
715     fhDcaV0Daughters[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("dcaV0Daughters%s",str)));
716     fhRadius[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("decayRadius%s",str)));
717     fhCosinePointingAngle[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("cosinepointingangle%s",str)));
718     fhOnFlyStatus[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("onflystatus%s",str)));
719     fhDcaV0ToVertex[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("dcaV0ToVertex%s",str)));
720
721     fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("pt%s",str)));
722
723     gDirectory->cd("../");
724   }
725
726   gDirectory->cd("..");
727
728   return kTRUE;
729 }
730
731 //____________________________________________________________________
732 void AliESDv0Cuts::SaveHistograms(const Char_t* dir) {
733   //
734   // saves the histograms in a directory (dir)
735   //
736
737   if (!fHistogramsOn) {
738     AliDebug(0, "Histograms not on - cannot save histograms!!!");
739     return;
740   }
741
742   if (!dir)
743     dir = GetName();
744
745   gDirectory->mkdir(dir);
746   gDirectory->cd(dir);
747
748   gDirectory->mkdir("before_cuts");
749   gDirectory->mkdir("after_cuts");
750  
751   fhCutStatistics->Write();
752   fhCutCorrelation->Write();
753
754   for (Int_t i=0; i<2; i++) {
755     if (i==0)
756       gDirectory->cd("before_cuts");
757     else
758       gDirectory->cd("after_cuts");
759
760     fhDcaPosToVertex[i]      ->Write();
761     fhDcaNegToVertex[i]      ->Write();
762     fhChi2[i]                ->Write();
763     fhDcaV0Daughters[i]      ->Write();
764     fhRadius[i]              ->Write();
765     fhCosinePointingAngle[i] ->Write();
766     fhOnFlyStatus[i]         ->Write();
767     fhDcaV0ToVertex[i]       ->Write();
768
769     fhPt[i]                  ->Write();
770     
771     gDirectory->cd("../");
772   }
773
774   gDirectory->cd("../");
775 }
776
777 //____________________________________________________________________
778 void AliESDv0Cuts::DrawHistograms()
779 {
780   // draws some histograms
781
782   TCanvas* canvas1 = new TCanvas(Form("%s_1", GetName()), "V0 Quality Results1", 800, 800);
783   canvas1->Divide(2, 2);
784
785   canvas1->cd(1);
786   fhDcaPosToVertex[0]->SetStats(kFALSE);
787   fhDcaPosToVertex[0]->Draw();
788
789   canvas1->cd(2);
790   fhChi2[0]->SetStats(kFALSE);
791   fhChi2[0]->Draw();
792
793   canvas1->cd(3);
794   fhDcaV0ToVertex[0]->SetStats(kFALSE);
795   fhDcaV0ToVertex[0]->Draw();
796
797   canvas1->SaveAs(Form("%s_%s.gif", GetName(), canvas1->GetName()));
798
799   TCanvas* canvas2 = new TCanvas(Form("%s_2", GetName()), "V0 Quality Results2", 1200, 800);
800   canvas2->Divide(3, 2);
801
802   canvas2->cd(1);
803   fhDcaV0Daughters[0]->SetStats(kFALSE);
804   gPad->SetLogy();
805   fhDcaV0Daughters[0]->Draw();
806
807   canvas2->cd(2);
808   fhRadius[0]->SetStats(kFALSE);
809   gPad->SetLogy();
810   fhRadius[0]->Draw();
811
812
813   canvas2->cd(4);
814   fhCosinePointingAngle[0]->SetStats(kFALSE);
815   gPad->SetLogy();
816   fhCosinePointingAngle[0]->Draw();
817
818   canvas2->cd(5);
819   fhOnFlyStatus[0]->SetStats(kFALSE);
820   gPad->SetLogy();
821   fhOnFlyStatus[0]->Draw();
822
823   canvas2->SaveAs(Form("%s_%s.gif", GetName(), canvas2->GetName()));
824
825   TCanvas* canvas3 = new TCanvas(Form("%s_4", GetName()), "V0 Quality Results3", 800, 500);
826   canvas3->Divide(2, 1);
827
828   canvas3->cd(1);
829   fhCutStatistics->SetStats(kFALSE);
830   fhCutStatistics->LabelsOption("v");
831   gPad->SetBottomMargin(0.3);
832   fhCutStatistics->Draw();
833
834   canvas3->cd(2);
835   fhCutCorrelation->SetStats(kFALSE);
836   fhCutCorrelation->LabelsOption("v");
837   gPad->SetBottomMargin(0.3);
838   gPad->SetLeftMargin(0.3);
839   fhCutCorrelation->Draw("COLZ");
840
841   canvas3->SaveAs(Form("%s_%s.gif", GetName(), canvas3->GetName()));
842 }
843