[u/mrichter/AliRoot.git] / PWG1 / AliAnalysisTaskIPInfo.cxx

/**************************************************************************
 * Copyright(c) 1998-2009, 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.                  *
 **************************************************************************/

//*************************************************************************
// Class AliAnalysisTaskIPInfo
// AliAnalysisTask to extract from the ESD the IP position and sigma
// as well as to estimate the primary vertex and tracks DCA resolution.
// Uses external class AliIntSpotEstimator
//
// Author: ruben.shahoyan@cern.ch
//*************************************************************************

#include <TChain.h>
#include <TTree.h>
#include <TNtuple.h>
#include <TBranch.h>
#include <TClonesArray.h>
#include <TObjArray.h>
#include <TH1F.h>
#include <TH2F.h>  
#include <TCanvas.h>

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"

#include "AliESDtrack.h"
#include "AliExternalTrackParam.h"
#include "AliESDVertex.h"
#include "AliESDEvent.h"
#include "AliVertexerTracks.h"
#include "AliESDInputHandler.h"
#include "AliAnalysisTaskIPInfo.h"
#include "AliIntSpotEstimator.h"


ClassImp(AliAnalysisTaskIPInfo)

//________________________________________________________________________
AliAnalysisTaskIPInfo::AliAnalysisTaskIPInfo(const char *name) 
: AliAnalysisTask(name, "IP analysis"),
  fESD(0),fOutput(0),fTracks(50)
{

  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #0 writes into a TList container
  DefineOutput(0, TList::Class());  //My private output
  //
  for (int i=0;i<kNEst;i++) { // default settings
    fIPEst[i] = 0;
    SetOptions(i);
  }
}

//________________________________________________________________________
AliAnalysisTaskIPInfo::~AliAnalysisTaskIPInfo()
{
  // Destructor
  if (fOutput) {
    delete fOutput;
    fOutput = 0;
  }
}

//________________________________________________________________________
void AliAnalysisTaskIPInfo::SetOptions(Int_t estID, Bool_t recoVtx,
				       Double_t outcut,Int_t nPhiBins,Int_t nestb,
				       Double_t estmin,Double_t estmax,
				       Int_t ntrBins,Int_t ntMn,Int_t ntMx,
				       Int_t nPBins,Double_t pmn,Double_t pmx)
{
  // set options for estimators
  if (estID<0 || estID>= kNEst) return;
  fRecoVtx[estID] = recoVtx;
  fNPhiBins[estID] = nPhiBins;
  fNEstb[estID]    = nestb;
  fNTrBins[estID]  = ntrBins;
  fNPBins[estID]   = nPBins;
  fNTrMin[estID]   = ntMn;
  fNTrMax[estID]   = ntMx;
  fOutCut[estID]   = outcut;
  fEstMin[estID]   = estmin;
  fEstMax[estID]   = estmax;
  fPMin[estID]     = pmn;
  fPMax[estID]     = pmx;
}

//________________________________________________________________________
void AliAnalysisTaskIPInfo::ConnectInputData(Option_t *) 
{
  // Connect ESD or AOD here
  // Called once
  //
  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  if (!tree) {
    Printf("ERROR: Could not read chain from input slot 0");
  } 
  else {
    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());    
    if (!esdH) Printf("ERROR: Could not get ESDInputHandler");
    else fESD = esdH->GetEvent();
  }
  //
  return;
}

//________________________________________________________________________
void AliAnalysisTaskIPInfo::CreateOutputObjects()
{
  // Create estimators
  // Several histograms are more conveniently managed in a TList
  fOutput = new TList;
  fOutput->SetOwner();
  //
  TString nm = GetName();
  nm += "_TPC";
  fIPEst[kTPC] = new AliIntSpotEstimator(nm.Data(),fOutCut[kTPC],fNPhiBins[kTPC],
					 fNEstb[kTPC],fEstMin[kTPC],fEstMax[kTPC],
					 fNTrBins[kTPC],fNTrMin[kTPC],fNTrMax[kTPC],
					 fNPBins[kTPC],fPMin[kTPC],fPMax[kTPC]);
  fIPEst[kTPC]->GetVertexer()->SetTPCMode();
  fIPEst[kTPC]->GetVertexer()->SetConstraintOff();
  //
  nm = GetName();
  nm += "_ITSTPC";
  fIPEst[kITS] = new AliIntSpotEstimator(nm.Data(),fOutCut[kITS],fNPhiBins[kITS],
					 fNEstb[kITS],fEstMin[kITS],fEstMax[kITS],
					 fNTrBins[kITS],fNTrMin[kITS],fNTrMax[kITS],
					 fNPBins[kITS],fPMin[kITS],fPMax[kITS]);
  //
  fIPEst[kITS]->GetVertexer()->SetITSMode();
  fIPEst[kITS]->GetVertexer()->SetConstraintOff();
  //
  fOutput->Add(fIPEst[kTPC]);
  fOutput->Add(fIPEst[kITS]);
  //
  return;
}

//________________________________________________________________________
void AliAnalysisTaskIPInfo::Exec(Option_t *) 
{
  // Main loop
  // Called for each event
  
  if (!fESD) {
    Printf("ERROR: fESD not available");
    return;
  }
  //
  fIPEst[kTPC]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
  fIPEst[kITS]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
  const AliESDVertex *vtx;
  UShort_t *trackID;
  Int_t ntracks;
  //
  // Process vertex made of TPC tracks only
  vtx = fRecoVtx[kTPC] ? fIPEst[kTPC]->GetVertexer()->FindPrimaryVertex(fESD) : fESD->GetPrimaryVertexTPC();
  if (vtx) {
    ntracks = vtx->GetNIndices();
    trackID = (UShort_t*)vtx->GetIndices();
    for (int i=ntracks;i--;) fTracks.Add((TObject*)fESD->GetTrack(trackID[i])->GetTPCInnerParam());
    fIPEst[kTPC]->ProcessEvent(&fTracks);
    fTracks.Clear();
  }
  //
  // Process vertex made of TPC+ITS tracks
  vtx = fRecoVtx[kITS] ? fIPEst[kITS]->GetVertexer()->FindPrimaryVertex(fESD) : fESD->GetPrimaryVertex();
  if (vtx) {
    ntracks = vtx->GetNIndices();
    trackID = (UShort_t*)vtx->GetIndices();
    for (int i=ntracks;i--;) fTracks.Add((TObject*)fESD->GetTrack(trackID[i]));
    fIPEst[kITS]->ProcessEvent(&fTracks);
    fTracks.Clear();
  }
  //  
  PostData(0, fOutput);
  //
  return;
}      

//________________________________________________________________________
void AliAnalysisTaskIPInfo::Terminate(Option_t *) 
{
  // Draw result to the screen
  // Called once at the end of the query
  fOutput = dynamic_cast<TList*> (GetOutputData(0));
  if (!fOutput) {     
    Printf("ERROR: fOutput not available");
    return;
  }
  //
  return;
}
Commit	Line	Data
5b09c01f	1	/**************************************************************************
	2	* Copyright(c) 1998-2009, 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	//*************************************************************************
	17	// Class AliAnalysisTaskIPInfo
	18	// AliAnalysisTask to extract from the ESD the IP position and sigma
	19	// as well as to estimate the primary vertex and tracks DCA resolution.
	20	// Uses external class AliIntSpotEstimator
	21	//
	22	// Author: ruben.shahoyan@cern.ch
	23	//*************************************************************************
	24
	25	#include <TChain.h>
	26	#include <TTree.h>
	27	#include <TNtuple.h>
	28	#include <TBranch.h>
	29	#include <TClonesArray.h>
	30	#include <TObjArray.h>
	31	#include <TH1F.h>
	32	#include <TH2F.h>
	33	#include <TCanvas.h>
	34
	35	#include "AliAnalysisTask.h"
	36	#include "AliAnalysisManager.h"
	37
	38	#include "AliESDtrack.h"
	39	#include "AliExternalTrackParam.h"
	40	#include "AliESDVertex.h"
	41	#include "AliESDEvent.h"
	42	#include "AliVertexerTracks.h"
	43	#include "AliESDInputHandler.h"
	44	#include "AliAnalysisTaskIPInfo.h"
	45	#include "AliIntSpotEstimator.h"
	46
	47
	48	ClassImp(AliAnalysisTaskIPInfo)
	49
	50	//________________________________________________________________________
	51	AliAnalysisTaskIPInfo::AliAnalysisTaskIPInfo(const char *name)
	52	: AliAnalysisTask(name, "IP analysis"),
	53	fESD(0),fOutput(0),fTracks(50)
	54	{
	55
	56	// Define input and output slots here
	57	// Input slot #0 works with a TChain
	58	DefineInput(0, TChain::Class());
	59	// Output slot #0 writes into a TList container
	60	DefineOutput(0, TList::Class()); //My private output
	61	//
	62	for (int i=0;i<kNEst;i++) { // default settings
	63	fIPEst[i] = 0;
	64	SetOptions(i);
65	}
66	}
67
68	//________________________________________________________________________
69	AliAnalysisTaskIPInfo::~AliAnalysisTaskIPInfo()
70	{
71	// Destructor
72	if (fOutput) {
73	delete fOutput;
74	fOutput = 0;
75	}
76	}
77
78	//________________________________________________________________________
79	void AliAnalysisTaskIPInfo::SetOptions(Int_t estID, Bool_t recoVtx,
80	Double_t outcut,Int_t nPhiBins,Int_t nestb,
81	Double_t estmin,Double_t estmax,
82	Int_t ntrBins,Int_t ntMn,Int_t ntMx,
83	Int_t nPBins,Double_t pmn,Double_t pmx)
84	{
85	// set options for estimators
86	if (estID<0 \|\| estID>= kNEst) return;
87	fRecoVtx[estID] = recoVtx;
88	fNPhiBins[estID] = nPhiBins;
89	fNEstb[estID] = nestb;
90	fNTrBins[estID] = ntrBins;
91	fNPBins[estID] = nPBins;
92	fNTrMin[estID] = ntMn;
93	fNTrMax[estID] = ntMx;
94	fOutCut[estID] = outcut;
95	fEstMin[estID] = estmin;
96	fEstMax[estID] = estmax;
97	fPMin[estID] = pmn;
98	fPMax[estID] = pmx;
99	}
100
101	//________________________________________________________________________
102	void AliAnalysisTaskIPInfo::ConnectInputData(Option_t *)
103	{
104	// Connect ESD or AOD here
105	// Called once
106	//
107	TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
108	if (!tree) {
109	Printf("ERROR: Could not read chain from input slot 0");
110	}
111	else {
112	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
113	if (!esdH) Printf("ERROR: Could not get ESDInputHandler");
114	else fESD = esdH->GetEvent();
115	}
116	//
117	return;
118	}
119
120	//________________________________________________________________________
121	void AliAnalysisTaskIPInfo::CreateOutputObjects()
122	{
123	// Create estimators
124	// Several histograms are more conveniently managed in a TList
125	fOutput = new TList;
126	fOutput->SetOwner();
127	//
128	TString nm = GetName();
129	nm += "_TPC";
130	fIPEst[kTPC] = new AliIntSpotEstimator(nm.Data(),fOutCut[kTPC],fNPhiBins[kTPC],
131	fNEstb[kTPC],fEstMin[kTPC],fEstMax[kTPC],
132	fNTrBins[kTPC],fNTrMin[kTPC],fNTrMax[kTPC],
133	fNPBins[kTPC],fPMin[kTPC],fPMax[kTPC]);
134	fIPEst[kTPC]->GetVertexer()->SetTPCMode();
135	fIPEst[kTPC]->GetVertexer()->SetConstraintOff();
136	//
137	nm = GetName();
138	nm += "_ITSTPC";
139	fIPEst[kITS] = new AliIntSpotEstimator(nm.Data(),fOutCut[kITS],fNPhiBins[kITS],
140	fNEstb[kITS],fEstMin[kITS],fEstMax[kITS],
141	fNTrBins[kITS],fNTrMin[kITS],fNTrMax[kITS],
142	fNPBins[kITS],fPMin[kITS],fPMax[kITS]);
143	//
144	fIPEst[kITS]->GetVertexer()->SetITSMode();
145	fIPEst[kITS]->GetVertexer()->SetConstraintOff();
146	//
147	fOutput->Add(fIPEst[kTPC]);
148	fOutput->Add(fIPEst[kITS]);
149	//
150	return;
151	}
152
153	//________________________________________________________________________
154	void AliAnalysisTaskIPInfo::Exec(Option_t *)
155	{
156	// Main loop
157	// Called for each event
158
159	if (!fESD) {
160	Printf("ERROR: fESD not available");
161	return;
162	}
163	//
164	fIPEst[kTPC]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
165	fIPEst[kITS]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
166	const AliESDVertex *vtx;
167	UShort_t *trackID;
168	Int_t ntracks;
169	//
170	// Process vertex made of TPC tracks only
171	vtx = fRecoVtx[kTPC] ? fIPEst[kTPC]->GetVertexer()->FindPrimaryVertex(fESD) : fESD->GetPrimaryVertexTPC();
172	if (vtx) {
173	ntracks = vtx->GetNIndices();
174	trackID = (UShort_t*)vtx->GetIndices();
175	for (int i=ntracks;i--;) fTracks.Add((TObject*)fESD->GetTrack(trackID[i])->GetTPCInnerParam());
176	fIPEst[kTPC]->ProcessEvent(&fTracks);
177	fTracks.Clear();
178	}
179	//
180	// Process vertex made of TPC+ITS tracks
181	vtx = fRecoVtx[kITS] ? fIPEst[kITS]->GetVertexer()->FindPrimaryVertex(fESD) : fESD->GetPrimaryVertex();
182	if (vtx) {
183	ntracks = vtx->GetNIndices();
184	trackID = (UShort_t*)vtx->GetIndices();
185	for (int i=ntracks;i--;) fTracks.Add((TObject*)fESD->GetTrack(trackID[i]));
186	fIPEst[kITS]->ProcessEvent(&fTracks);
187	fTracks.Clear();
188	}
189	//
190	PostData(0, fOutput);
191	//
192	return;
193	}
194
195	//________________________________________________________________________
196	void AliAnalysisTaskIPInfo::Terminate(Option_t *)
197	{
198	// Draw result to the screen
199	// Called once at the end of the query
200	fOutput = dynamic_cast<TList*> (GetOutputData(0));
201	if (!fOutput) {
202	Printf("ERROR: fOutput not available");
203	return;
204	}
205	//
206	return;
207	}