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