[u/mrichter/AliRoot.git] / PWG3 / 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.                  *
 **************************************************************************/

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