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) sprintf(str," ");
 647     else sprintf(str,"_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       sprintf(str,"_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