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