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