[u/mrichter/AliRoot.git] / PWG / muon / AliESDMuonTrackCuts.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

#include "AliESDMuonTrackCuts.h"

#include <AliESDMuonTrack.h>
#include <AliESD.h>
#include <AliESDEvent.h>
#include <AliLog.h>

#include <TTree.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TDirectory.h>

//____________________________________________________________________
ClassImp(AliESDMuonTrackCuts)

// Cut names
const Char_t* AliESDMuonTrackCuts::fgkCutNames[kNCuts] = {
 "p",
 "p_{T}",
 "p_{x}",
 "p_{y}",
 "p_{z}",
 "y",
 "eta"
};

//____________________________________________________________________
AliESDMuonTrackCuts::AliESDMuonTrackCuts(const Char_t* name, const Char_t* title) : AliAnalysisCuts(name,title),
  fPMin(0),
  fPMax(0),
  fPtMin(0),
  fPtMax(0),
  fPxMin(0),
  fPxMax(0),
  fPyMin(0),
  fPyMax(0),
  fPzMin(0),
  fPzMax(0),
  fEtaMin(0),
  fEtaMax(0),
  fRapMin(0),
  fRapMax(0),
  fHistogramsOn(0),
  fhCutStatistics(0),         
  fhCutCorrelation(0)
{
  //
  // constructor
  //
  Init();

  //##############################################################################
  // setting default cuts
  SetPRange();
  SetPtRange();
  SetPxRange();
  SetPyRange();
  SetPzRange();
  SetEtaRange();
  SetRapRange();

  SetHistogramsOn();
}

//_____________________________________________________________________________
AliESDMuonTrackCuts::AliESDMuonTrackCuts(const AliESDMuonTrackCuts &c) : AliAnalysisCuts(c),
  fPMin(0),
  fPMax(0),
  fPtMin(0),
  fPtMax(0),
  fPxMin(0),
  fPxMax(0),
  fPyMin(0),
  fPyMax(0),
  fPzMin(0),
  fPzMax(0),
  fEtaMin(0),
  fEtaMax(0),
  fRapMin(0),
  fRapMax(0),
  fHistogramsOn(0),
  fhCutStatistics(0),         
  fhCutCorrelation(0)
{
  //
  // copy constructor
  //
  ((AliESDMuonTrackCuts &) c).Copy(*this);
}

AliESDMuonTrackCuts::~AliESDMuonTrackCuts()
{
  //
  // destructor
  //
  for (Int_t i=0; i<2; i++) {    
    if (fhPt[i])
      delete fhPt[i];
    if (fhEta[i])
      delete fhEta[i];
  }

  if (fhCutStatistics)
    delete fhCutStatistics;             
  if (fhCutCorrelation)
    delete fhCutCorrelation;            
}

void AliESDMuonTrackCuts::Init()
{
  //
  // sets everything to zero
  //
  fPMin = 0;
  fPMax = 0;
  fPtMin = 0;
  fPtMax = 0;
  fPxMin = 0;
  fPxMax = 0;
  fPyMin = 0;
  fPyMax = 0;
  fPzMin = 0;
  fPzMax = 0;
  fEtaMin = 0;
  fEtaMax = 0;
  fRapMin = 0;
  fRapMax = 0;

  fHistogramsOn = kFALSE;

  for (Int_t i=0; i<2; ++i)
  {
    fhPt[i] = 0;
    fhEta[i] = 0;
  }
  fhCutStatistics = 0;
  fhCutCorrelation = 0;
}

//_____________________________________________________________________________
AliESDMuonTrackCuts &AliESDMuonTrackCuts::operator=(const AliESDMuonTrackCuts &c)
{
  //
  // Assignment operator
  //
  if (this != &c) ((AliESDMuonTrackCuts &) c).Copy(*this);
  return *this;
}

//_____________________________________________________________________________
void AliESDMuonTrackCuts::Copy(TObject &c) const
{
  //
  // Copy function
  //

   AliESDMuonTrackCuts& target = (AliESDMuonTrackCuts &) c;
 
  target.Init();
// 
  target.fPMin = fPMin;
  target.fPMax = fPMax;
  target.fPtMin = fPtMin;
  target.fPtMax = fPtMax;
  target.fPxMin = fPxMin;
  target.fPxMax = fPxMax;
  target.fPyMin = fPyMin;
  target.fPyMax = fPyMax;
  target.fPzMin = fPzMin;
  target.fPzMax = fPzMax;
  target.fEtaMin = fEtaMin;
  target.fEtaMax = fEtaMax;
  target.fRapMin = fRapMin;
  target.fRapMax = fRapMax;

  target.fHistogramsOn = fHistogramsOn;

  for (Int_t i=0; i<2; ++i)
  {     
    if (fhPt[i]) target.fhPt[i] = (TH1F*) fhPt[i]->Clone();
    if (fhEta[i]) target.fhEta[i] = (TH1F*) fhEta[i]->Clone();
  }

  if (fhCutStatistics) target.fhCutStatistics = (TH1F*) fhCutStatistics->Clone();
  if (fhCutCorrelation) target.fhCutCorrelation = (TH2F*) fhCutCorrelation->Clone();

  TNamed::Copy(c);
}

//_____________________________________________________________________________
Long64_t AliESDMuonTrackCuts::Merge(TCollection* list) {
  // Merge a list of AliESDMuonTrackCuts objects with this (needed for PROOF)
  // Returns the number of merged objects (including this)

  if (!list)
    return 0;
  
  if (list->IsEmpty())
    return 1;

  if (!fHistogramsOn)
    return 0;

  TIterator* iter = list->MakeIterator();
  TObject* obj;

  // collection of measured and generated histograms
  Int_t count = 0;
  while ((obj = iter->Next())) {

    AliESDMuonTrackCuts* entry = dynamic_cast<AliESDMuonTrackCuts*>(obj);
    if (entry == 0)
      continue;

    if (!entry->fHistogramsOn)
      continue;
    
    for (Int_t i=0; i<2; i++) {
      fhPt[i]->Add(entry->fhPt[i]); 
      fhEta[i]->Add(entry->fhEta[i]); 
    }      

    fhCutStatistics->Add(entry->fhCutStatistics);        
    fhCutCorrelation ->Add(entry->fhCutCorrelation);      

    count++;
  }

  return count+1;
}

void AliESDMuonTrackCuts::EnableNeededBranches(TTree* tree)
{
  // enables the branches needed by AcceptTrack, for a list see comment of AcceptTrack

  tree->SetBranchStatus("fTracks.fFlags", 1);
  tree->SetBranchStatus("fTracks.fP*", 1);
  tree->SetBranchStatus("fTracks.fR*", 1);  //detector response probability
}

//____________________________________________________________________
Bool_t
AliESDMuonTrackCuts::AcceptTrack(AliESDMuonTrack* esdMuTrack) {
  // 
  // figure out if the tracks survives all the track cuts defined
  //
  // the different kinematic values are first
  // retrieved from the track. then it is found out what cuts the
  // track did not survive and finally the cuts are imposed.

  // getting the kinematic variables of the track
  // (assuming the mass is known)
  Double_t p[3];
  esdMuTrack->PxPyPz(p);
  Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
  Float_t pt       = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
  Float_t energy   = TMath::Sqrt(TMath::Power(esdMuTrack->M(),2) + TMath::Power(momentum,2));


  //y-eta related calculations
  Float_t eta = -100.;
  Float_t y   = -100.;
  if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
    eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
  if((energy != TMath::Abs(p[2]))&&(momentum != 0))
    y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));
     
  //########################################################################
  // cut the track?
  
  Bool_t cuts[kNCuts];
  for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
  
  // track kinematics cut
  if((momentum < fPMin) || (momentum > fPMax)) 
    cuts[0]=kTRUE;
  if((pt < fPtMin) || (pt > fPtMax)) 
    cuts[1] = kTRUE;
  if((p[0] < fPxMin) || (p[0] > fPxMax)) 
    cuts[2] = kTRUE;
  if((p[1] < fPyMin) || (p[1] > fPyMax)) 
    cuts[3] = kTRUE;
  if((p[2] < fPzMin) || (p[2] > fPzMax))
    cuts[4] = kTRUE;
  if((eta < fEtaMin) || (eta > fEtaMax)) 
    cuts[5] = kTRUE;
  if((y < fRapMin) || (y > fRapMax)) 
    cuts[6] = kTRUE;

  Bool_t cut=kFALSE;
  for (Int_t i=0; i<kNCuts; i++) 
    if (cuts[i]) cut = kTRUE;
  
  //########################################################################
  // filling histograms
  if (fHistogramsOn) {
    fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
    
    if (cut)
      fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
    
    for (Int_t i=0; i<kNCuts; i++) {
      if (cuts[i])
 	fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
      
      for (Int_t j=i; j<kNCuts; j++) {
 	if (cuts[i] && cuts[j]) {
 	  Float_t xC = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
 	  Float_t yC = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
 	  fhCutCorrelation->Fill(xC, yC);
 	}
      }
    }
    
    fhPt[0]->Fill(pt);
    fhEta[0]->Fill(eta);

  }

  //########################################################################
  // cut the track!
  if (cut) return kFALSE;

  //########################################################################
  // filling histograms after cut
  if (fHistogramsOn) {
// 
    fhPt[1]->Fill(pt);
    fhEta[1]->Fill(eta);
    
  }

  return kTRUE;
}

//____________________________________________________________________
TObjArray* AliESDMuonTrackCuts::GetAcceptedTracks(AliESD* esd)
{
  //
  // returns an array of all tracks that pass the cuts
  //

  TObjArray* acceptedTracks = new TObjArray();

  // loop over esd tracks
  for (Int_t iTrack = 0; iTrack < esd->GetNumberOfMuonTracks(); iTrack++) {
    AliESDMuonTrack* track = esd->GetMuonTrack(iTrack);

    if (AcceptTrack(track))
      acceptedTracks->Add(track);
  }

  return acceptedTracks;
}

//____________________________________________________________________
Int_t AliESDMuonTrackCuts::CountAcceptedTracks(AliESD* esd)
{
  //
  // returns an the number of tracks that pass the cuts
  //

  Int_t count = 0;

  // loop over esd tracks
  for (Int_t iTrack = 0; iTrack < esd->GetNumberOfMuonTracks(); iTrack++) {
    AliESDMuonTrack* track = esd->GetMuonTrack(iTrack);

    if (AcceptTrack(track))
      count++;
  }

  return count;
}

//____________________________________________________________________
TObjArray* AliESDMuonTrackCuts::GetAcceptedTracks(AliESDEvent* esd)
{
  //
  // returns an array of all tracks that pass the cuts
  //

  TObjArray* acceptedTracks = new TObjArray();

  // loop over esd tracks
  for (Int_t iTrack = 0; iTrack < esd->GetNumberOfMuonTracks(); iTrack++) {
    AliESDMuonTrack* track = esd->GetMuonTrack(iTrack);

    if (AcceptTrack(track))
      acceptedTracks->Add(track);
  }

  return acceptedTracks;
}

//____________________________________________________________________
Int_t AliESDMuonTrackCuts::CountAcceptedTracks(AliESDEvent* esd)
{
  //
  // returns an the number of tracks that pass the cuts
  //

  Int_t count = 0;

  // loop over esd tracks
  for (Int_t iTrack = 0; iTrack < esd->GetNumberOfMuonTracks(); iTrack++) {
    AliESDMuonTrack* track = esd->GetMuonTrack(iTrack);

    if (AcceptTrack(track))
      count++;
  }

  return count;
}

//____________________________________________________________________
 void AliESDMuonTrackCuts::DefineHistograms(Int_t color) {
   // 
   // diagnostics histograms are defined
   // 

   fHistogramsOn=kTRUE;
   
   Bool_t oldStatus = TH1::AddDirectoryStatus();
   TH1::AddDirectory(kFALSE);
   
   //###################################################################################
   // defining histograms

   fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);

   fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
   fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");

   fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
  
   for (Int_t i=0; i<kNCuts; i++) {
     fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
     fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
     fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
   } 

  fhCutStatistics  ->SetLineColor(color);
  fhCutCorrelation ->SetLineColor(color);
  fhCutStatistics  ->SetLineWidth(2);
  fhCutCorrelation ->SetLineWidth(2);

  Char_t str[256];
  for (Int_t i=0; i<2; i++) {
    if (i==0) snprintf(str,256," ");
    else snprintf(str,256,"_cut");

    fhPt[i]                  = new TH1F(Form("pt%s",str)     ,"p_{T} distribution;p_{T} (GeV/c)",500,0.0,100.0);
    fhEta[i]                 = new TH1F(Form("eta%s",str)     ,"#eta distribution;#eta",40,-2.0,2.0);
  }

  TH1::AddDirectory(oldStatus);
}

//____________________________________________________________________
Bool_t AliESDMuonTrackCuts::LoadHistograms(const Char_t* dir)
{
  //
  // loads the histograms from a file
  // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
  //

  if (!dir)
    dir = GetName();

  if (!gDirectory->cd(dir))
    return kFALSE;

  fhCutStatistics = dynamic_cast<TH1F*> (gDirectory->Get("cut_statistics"));
  fhCutCorrelation = dynamic_cast<TH2F*> (gDirectory->Get("cut_correlation"));

  Char_t str[5];
  for (Int_t i=0; i<2; i++) {
    if (i==0)
    {
      gDirectory->cd("before_cuts");
      str[0] = 0;
    }
    else
    {
      gDirectory->cd("after_cuts");
      snprintf(str,5,"_cut");
    }

    fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("pt%s",str)));
    fhEta[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("eta%s",str)));

    gDirectory->cd("../");
  }

  gDirectory->cd("..");

  return kTRUE;
}

//____________________________________________________________________
void AliESDMuonTrackCuts::SaveHistograms(const Char_t* dir) {
  //
  // saves the histograms in a directory (dir)
  //

  if (!fHistogramsOn) {
    AliDebug(0, "Histograms not on - cannot save histograms!!!");
    return;
  }

  if (!dir)
    dir = GetName();

  gDirectory->mkdir(dir);
  gDirectory->cd(dir);

  gDirectory->mkdir("before_cuts");
  gDirectory->mkdir("after_cuts");

  fhCutStatistics->Write();
  fhCutCorrelation->Write();

  for (Int_t i=0; i<2; i++) {
    if (i==0)
      gDirectory->cd("before_cuts");
    else
      gDirectory->cd("after_cuts");

    fhPt[i]                  ->Write();
    fhEta[i]                 ->Write();
    
    gDirectory->cd("../");
  }

  gDirectory->cd("../");
}

//____________________________________________________________________
void AliESDMuonTrackCuts::DrawHistograms()
{
  gStyle->SetPalette(1);
  gStyle->SetFrameFillColor(10);
  gStyle->SetCanvasColor(10);
  
  TCanvas* canvas1 = new TCanvas(Form("%s_1", GetName()), "Track Cut Results", 800, 500);
  canvas1->Divide(2, 1);

  canvas1->cd(1);
  fhCutStatistics->SetStats(kFALSE);
  fhCutStatistics->LabelsOption("v");
  gPad->SetBottomMargin(0.3);
  fhCutStatistics->Draw();

  canvas1->cd(2);
  fhCutCorrelation->SetStats(kFALSE);
  fhCutCorrelation->LabelsOption("v");
  gPad->SetBottomMargin(0.3);
  gPad->SetLeftMargin(0.3);
  fhCutCorrelation->Draw("COLZ");
  canvas1->Update();
  canvas1->SaveAs(Form("%s_%s.gif", GetName(), canvas1->GetName()));
  
}
Commit	Line	Data
f70a1b5d	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
27de2dfb	16	/* $Id$ */
27de2dfb	17
f70a1b5d	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++) {
fa48f7eb	467	if (i==0) snprintf(str,256," ");
fa48f7eb	468	else snprintf(str,256,"_cut");
f70a1b5d	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");
fa48f7eb	504	snprintf(str,5,"_cut");
f70a1b5d	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