PWG/muon/AliESDMuonTrackCuts.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$ */
  17
  18 #include "AliESDMuonTrackCuts.h"
  19
  20 #include <AliESDMuonTrack.h>
  21 #include <AliESD.h>
  22 #include <AliESDEvent.h>
  23 #include <AliLog.h>
  24
  25 #include <TTree.h>
  26 #include <TStyle.h>
  27 #include <TCanvas.h>
  28 #include <TDirectory.h>
  29
  30 //____________________________________________________________________
  31 ClassImp(AliESDMuonTrackCuts)
  32
  33 // Cut names
  34 const Char_t* AliESDMuonTrackCuts::fgkCutNames[kNCuts] = {
  35  "p",
  36  "p_{T}",
  37  "p_{x}",
  38  "p_{y}",
  39  "p_{z}",
  40  "y",
  41  "eta"
  42 };
  43
  44 //____________________________________________________________________
  45 AliESDMuonTrackCuts::AliESDMuonTrackCuts(const Char_t* name, const Char_t* title) : AliAnalysisCuts(name,title),
  46   fPMin(0),
  47   fPMax(0),
  48   fPtMin(0),
  49   fPtMax(0),
  50   fPxMin(0),
  51   fPxMax(0),
  52   fPyMin(0),
  53   fPyMax(0),
  54   fPzMin(0),
  55   fPzMax(0),
  56   fEtaMin(0),
  57   fEtaMax(0),
  58   fRapMin(0),
  59   fRapMax(0),
  60   fHistogramsOn(0),
  61   fhCutStatistics(0),
  62   fhCutCorrelation(0)
  63 {
  64   //
  65   // constructor
  66   //
  67   Init();
  68
  69   //##############################################################################
  70   // setting default cuts
  71   SetPRange();
  72   SetPtRange();
  73   SetPxRange();
  74   SetPyRange();
  75   SetPzRange();
  76   SetEtaRange();
  77   SetRapRange();
  78
  79   SetHistogramsOn();
  80 }
  81
  82 //_____________________________________________________________________________
  83 AliESDMuonTrackCuts::AliESDMuonTrackCuts(const AliESDMuonTrackCuts &c) : AliAnalysisCuts(c),
  84   fPMin(0),
  85   fPMax(0),
  86   fPtMin(0),
  87   fPtMax(0),
  88   fPxMin(0),
  89   fPxMax(0),
  90   fPyMin(0),
  91   fPyMax(0),
  92   fPzMin(0),
  93   fPzMax(0),
  94   fEtaMin(0),
  95   fEtaMax(0),
  96   fRapMin(0),
  97   fRapMax(0),
  98   fHistogramsOn(0),
  99   fhCutStatistics(0),
 100   fhCutCorrelation(0)
 101 {
 102   //
 103   // copy constructor
 104   //
 105   ((AliESDMuonTrackCuts &) c).Copy(*this);
 106 }
 107
 108 AliESDMuonTrackCuts::~AliESDMuonTrackCuts()
 109 {
 110   //
 111   // destructor
 112   //
 113   for (Int_t i=0; i<2; i++) {
 114     if (fhPt[i])
 115       delete fhPt[i];
 116     if (fhEta[i])
 117       delete fhEta[i];
 118   }
 119
 120   if (fhCutStatistics)
 121     delete fhCutStatistics;
 122   if (fhCutCorrelation)
 123     delete fhCutCorrelation;
 124 }
 125
 126 void AliESDMuonTrackCuts::Init()
 127 {
 128   //
 129   // sets everything to zero
 130   //
 131   fPMin = 0;
 132   fPMax = 0;
 133   fPtMin = 0;
 134   fPtMax = 0;
 135   fPxMin = 0;
 136   fPxMax = 0;
 137   fPyMin = 0;
 138   fPyMax = 0;
 139   fPzMin = 0;
 140   fPzMax = 0;
 141   fEtaMin = 0;
 142   fEtaMax = 0;
 143   fRapMin = 0;
 144   fRapMax = 0;
 145
 146   fHistogramsOn = kFALSE;
 147
 148   for (Int_t i=0; i<2; ++i)
 149   {
 150     fhPt[i] = 0;
 151     fhEta[i] = 0;
 152   }
 153   fhCutStatistics = 0;
 154   fhCutCorrelation = 0;
 155 }
 156
 157 //_____________________________________________________________________________
 158 AliESDMuonTrackCuts &AliESDMuonTrackCuts::operator=(const AliESDMuonTrackCuts &c)
 159 {
 160   //
 161   // Assignment operator
 162   //
 163   if (this != &c) ((AliESDMuonTrackCuts &) c).Copy(*this);
 164   return *this;
 165 }
 166
 167 //_____________________________________________________________________________
 168 void AliESDMuonTrackCuts::Copy(TObject &c) const
 169 {
 170   //
 171   // Copy function
 172   //
 173
 174    AliESDMuonTrackCuts& target = (AliESDMuonTrackCuts &) c;
 175
 176   target.Init();
 177 //
 178   target.fPMin = fPMin;
 179   target.fPMax = fPMax;
 180   target.fPtMin = fPtMin;
 181   target.fPtMax = fPtMax;
 182   target.fPxMin = fPxMin;
 183   target.fPxMax = fPxMax;
 184   target.fPyMin = fPyMin;
 185   target.fPyMax = fPyMax;
 186   target.fPzMin = fPzMin;
 187   target.fPzMax = fPzMax;
 188   target.fEtaMin = fEtaMin;
 189   target.fEtaMax = fEtaMax;
 190   target.fRapMin = fRapMin;
 191   target.fRapMax = fRapMax;
 192
 193   target.fHistogramsOn = fHistogramsOn;
 194
 195   for (Int_t i=0; i<2; ++i)
 196   {
 197     if (fhPt[i]) target.fhPt[i] = (TH1F*) fhPt[i]->Clone();
 198     if (fhEta[i]) target.fhEta[i] = (TH1F*) fhEta[i]->Clone();
 199   }
 200
 201   if (fhCutStatistics) target.fhCutStatistics = (TH1F*) fhCutStatistics->Clone();
 202   if (fhCutCorrelation) target.fhCutCorrelation = (TH2F*) fhCutCorrelation->Clone();
 203
 204   TNamed::Copy(c);
 205 }
 206
 207 //_____________________________________________________________________________
 208 Long64_t AliESDMuonTrackCuts::Merge(TCollection* list) {
 209   // Merge a list of AliESDMuonTrackCuts objects with this (needed for PROOF)
 210   // Returns the number of merged objects (including this)
 211
 212   if (!list)
 213     return 0;
 214
 215   if (list->IsEmpty())
 216     return 1;
 217
 218   if (!fHistogramsOn)
 219     return 0;
 220
 221   TIterator* iter = list->MakeIterator();
 222   TObject* obj;
 223
 224   // collection of measured and generated histograms
 225   Int_t count = 0;
 226   while ((obj = iter->Next())) {
 227
 228     AliESDMuonTrackCuts* entry = dynamic_cast<AliESDMuonTrackCuts*>(obj);
 229     if (entry == 0)
 230       continue;
 231
 232     if (!entry->fHistogramsOn)
 233       continue;
 234
 235     for (Int_t i=0; i<2; i++) {
 236       fhPt[i]->Add(entry->fhPt[i]);
 237       fhEta[i]->Add(entry->fhEta[i]);
 238     }
 239
 240     fhCutStatistics->Add(entry->fhCutStatistics);
 241     fhCutCorrelation ->Add(entry->fhCutCorrelation);
 242
 243     count++;
 244   }
 245
 246   return count+1;
 247 }
 248
 249 void AliESDMuonTrackCuts::EnableNeededBranches(TTree* tree)
 250 {
 251   // enables the branches needed by AcceptTrack, for a list see comment of AcceptTrack
 252
 253   tree->SetBranchStatus("fTracks.fFlags", 1);
 254   tree->SetBranchStatus("fTracks.fP*", 1);
 255   tree->SetBranchStatus("fTracks.fR*", 1);  //detector response probability
 256 }
 257
 258 //____________________________________________________________________
 259 Bool_t
 260 AliESDMuonTrackCuts::AcceptTrack(AliESDMuonTrack* esdMuTrack) {
 261   //
 262   // figure out if the tracks survives all the track cuts defined
 263   //
 264   // the different kinematic values are first
 265   // retrieved from the track. then it is found out what cuts the
 266   // track did not survive and finally the cuts are imposed.
 267
 268   // getting the kinematic variables of the track
 269   // (assuming the mass is known)
 270   Double_t p[3];
 271   esdMuTrack->PxPyPz(p);
 272   Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
 273   Float_t pt       = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
 274   Float_t energy   = TMath::Sqrt(TMath::Power(esdMuTrack->M(),2) + TMath::Power(momentum,2));
 275
 276
 277   //y-eta related calculations
 278   Float_t eta = -100.;
 279   Float_t y   = -100.;
 280   if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
 281     eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
 282   if((energy != TMath::Abs(p[2]))&&(momentum != 0))
 283     y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));
 284
 285   //########################################################################
 286   // cut the track?
 287
 288   Bool_t cuts[kNCuts];
 289   for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
 290
 291   // track kinematics cut
 292   if((momentum < fPMin) || (momentum > fPMax))
 293     cuts[0]=kTRUE;
 294   if((pt < fPtMin) || (pt > fPtMax))
 295     cuts[1] = kTRUE;
 296   if((p[0] < fPxMin) || (p[0] > fPxMax))
 297     cuts[2] = kTRUE;
 298   if((p[1] < fPyMin) || (p[1] > fPyMax))
 299     cuts[3] = kTRUE;
 300   if((p[2] < fPzMin) || (p[2] > fPzMax))
 301     cuts[4] = kTRUE;
 302   if((eta < fEtaMin) || (eta > fEtaMax))
 303     cuts[5] = kTRUE;
 304   if((y < fRapMin) || (y > fRapMax))
 305     cuts[6] = kTRUE;
 306
 307   Bool_t cut=kFALSE;
 308   for (Int_t i=0; i<kNCuts; i++)
 309     if (cuts[i]) cut = kTRUE;
 310
 311   //########################################################################
 312   // filling histograms
 313   if (fHistogramsOn) {
 314     fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
 315
 316     if (cut)
 317       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
 318
 319     for (Int_t i=0; i<kNCuts; i++) {
 320       if (cuts[i])
 321         fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
 322
 323       for (Int_t j=i; j<kNCuts; j++) {
 324         if (cuts[i] && cuts[j]) {
 325           Float_t xC = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
 326           Float_t yC = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
 327           fhCutCorrelation->Fill(xC, yC);
 328         }
 329       }
 330     }
 331
 332     fhPt[0]->Fill(pt);
 333     fhEta[0]->Fill(eta);
 334
 335   }
 336
 337   //########################################################################
 338   // cut the track!
 339   if (cut) return kFALSE;
 340
 341   //########################################################################
 342   // filling histograms after cut
 343   if (fHistogramsOn) {
 344 //
 345     fhPt[1]->Fill(pt);
 346     fhEta[1]->Fill(eta);
 347
 348   }
 349
 350   return kTRUE;
 351 }
 352
 353 //____________________________________________________________________
 354 TObjArray* AliESDMuonTrackCuts::GetAcceptedTracks(AliESD* esd)
 355 {
 356   //
 357   // returns an array of all tracks that pass the cuts
 358   //
 359
 360   TObjArray* acceptedTracks = new TObjArray();
 361
 362   // loop over esd tracks
 363   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfMuonTracks(); iTrack++) {
 364     AliESDMuonTrack* track = esd->GetMuonTrack(iTrack);
 365
 366     if (AcceptTrack(track))
 367       acceptedTracks->Add(track);
 368   }
 369
 370   return acceptedTracks;
 371 }
 372
 373 //____________________________________________________________________
 374 Int_t AliESDMuonTrackCuts::CountAcceptedTracks(AliESD* esd)
 375 {
 376   //
 377   // returns an the number of tracks that pass the cuts
 378   //
 379
 380   Int_t count = 0;
 381
 382   // loop over esd tracks
 383   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfMuonTracks(); iTrack++) {
 384     AliESDMuonTrack* track = esd->GetMuonTrack(iTrack);
 385
 386     if (AcceptTrack(track))
 387       count++;
 388   }
 389
 390   return count;
 391 }
 392
 393 //____________________________________________________________________
 394 TObjArray* AliESDMuonTrackCuts::GetAcceptedTracks(AliESDEvent* esd)
 395 {
 396   //
 397   // returns an array of all tracks that pass the cuts
 398   //
 399
 400   TObjArray* acceptedTracks = new TObjArray();
 401
 402   // loop over esd tracks
 403   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfMuonTracks(); iTrack++) {
 404     AliESDMuonTrack* track = esd->GetMuonTrack(iTrack);
 405
 406     if (AcceptTrack(track))
 407       acceptedTracks->Add(track);
 408   }
 409
 410   return acceptedTracks;
 411 }
 412
 413 //____________________________________________________________________
 414 Int_t AliESDMuonTrackCuts::CountAcceptedTracks(AliESDEvent* esd)
 415 {
 416   //
 417   // returns an the number of tracks that pass the cuts
 418   //
 419
 420   Int_t count = 0;
 421
 422   // loop over esd tracks
 423   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfMuonTracks(); iTrack++) {
 424     AliESDMuonTrack* track = esd->GetMuonTrack(iTrack);
 425
 426     if (AcceptTrack(track))
 427       count++;
 428   }
 429
 430   return count;
 431 }
 432
 433 //____________________________________________________________________
 434  void AliESDMuonTrackCuts::DefineHistograms(Int_t color) {
 435    //
 436    // diagnostics histograms are defined
 437    //
 438
 439    fHistogramsOn=kTRUE;
 440
 441    Bool_t oldStatus = TH1::AddDirectoryStatus();
 442    TH1::AddDirectory(kFALSE);
 443
 444    //###################################################################################
 445    // defining histograms
 446
 447    fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
 448
 449    fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
 450    fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
 451
 452    fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
 453
 454    for (Int_t i=0; i<kNCuts; i++) {
 455      fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
 456      fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
 457      fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
 458    }
 459
 460   fhCutStatistics  ->SetLineColor(color);
 461   fhCutCorrelation ->SetLineColor(color);
 462   fhCutStatistics  ->SetLineWidth(2);
 463   fhCutCorrelation ->SetLineWidth(2);
 464
 465   Char_t str[256];
 466   for (Int_t i=0; i<2; i++) {
 467     if (i==0) snprintf(str,256," ");
 468     else snprintf(str,256,"_cut");
 469
 470     fhPt[i]                  = new TH1F(Form("pt%s",str)     ,"p_{T} distribution;p_{T} (GeV/c)",500,0.0,100.0);
 471     fhEta[i]                 = new TH1F(Form("eta%s",str)     ,"#eta distribution;#eta",40,-2.0,2.0);
 472   }
 473
 474   TH1::AddDirectory(oldStatus);
 475 }
 476
 477 //____________________________________________________________________
 478 Bool_t AliESDMuonTrackCuts::LoadHistograms(const Char_t* dir)
 479 {
 480   //
 481   // loads the histograms from a file
 482   // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
 483   //
 484
 485   if (!dir)
 486     dir = GetName();
 487
 488   if (!gDirectory->cd(dir))
 489     return kFALSE;
 490
 491   fhCutStatistics = dynamic_cast<TH1F*> (gDirectory->Get("cut_statistics"));
 492   fhCutCorrelation = dynamic_cast<TH2F*> (gDirectory->Get("cut_correlation"));
 493
 494   Char_t str[5];
 495   for (Int_t i=0; i<2; i++) {
 496     if (i==0)
 497     {
 498       gDirectory->cd("before_cuts");
 499       str[0] = 0;
 500     }
 501     else
 502     {
 503       gDirectory->cd("after_cuts");
 504       snprintf(str,5,"_cut");
 505     }
 506
 507     fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("pt%s",str)));
 508     fhEta[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("eta%s",str)));
 509
 510     gDirectory->cd("../");
 511   }
 512
 513   gDirectory->cd("..");
 514
 515   return kTRUE;
 516 }
 517
 518 //____________________________________________________________________
 519 void AliESDMuonTrackCuts::SaveHistograms(const Char_t* dir) {
 520   //
 521   // saves the histograms in a directory (dir)
 522   //
 523
 524   if (!fHistogramsOn) {
 525     AliDebug(0, "Histograms not on - cannot save histograms!!!");
 526     return;
 527   }
 528
 529   if (!dir)
 530     dir = GetName();
 531
 532   gDirectory->mkdir(dir);
 533   gDirectory->cd(dir);
 534
 535   gDirectory->mkdir("before_cuts");
 536   gDirectory->mkdir("after_cuts");
 537
 538   fhCutStatistics->Write();
 539   fhCutCorrelation->Write();
 540
 541   for (Int_t i=0; i<2; i++) {
 542     if (i==0)
 543       gDirectory->cd("before_cuts");
 544     else
 545       gDirectory->cd("after_cuts");
 546
 547     fhPt[i]                  ->Write();
 548     fhEta[i]                 ->Write();
 549
 550     gDirectory->cd("../");
 551   }
 552
 553   gDirectory->cd("../");
 554 }
 555
 556 //____________________________________________________________________
 557 void AliESDMuonTrackCuts::DrawHistograms()
 558 {
 559   gStyle->SetPalette(1);
 560   gStyle->SetFrameFillColor(10);
 561   gStyle->SetCanvasColor(10);
 562
 563   TCanvas* canvas1 = new TCanvas(Form("%s_1", GetName()), "Track Cut Results", 800, 500);
 564   canvas1->Divide(2, 1);
 565
 566   canvas1->cd(1);
 567   fhCutStatistics->SetStats(kFALSE);
 568   fhCutStatistics->LabelsOption("v");
 569   gPad->SetBottomMargin(0.3);
 570   fhCutStatistics->Draw();
 571
 572   canvas1->cd(2);
 573   fhCutCorrelation->SetStats(kFALSE);
 574   fhCutCorrelation->LabelsOption("v");
 575   gPad->SetBottomMargin(0.3);
 576   gPad->SetLeftMargin(0.3);
 577   fhCutCorrelation->Draw("COLZ");
 578   canvas1->Update();
 579   canvas1->SaveAs(Form("%s_%s.gif", GetName(), canvas1->GetName()));
 580
 581 }
 582