Adding more bins in QA (Alis)
[u/mrichter/AliRoot.git] / PWGPP / AliAnalysisTaskIPInfo.cxx
CommitLineData
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>
d631e5e7 33#include <TNtuple.h>
5b09c01f 34#include <TCanvas.h>
35
36#include "AliAnalysisTask.h"
37#include "AliAnalysisManager.h"
38
39#include "AliESDtrack.h"
40#include "AliExternalTrackParam.h"
41#include "AliESDVertex.h"
42#include "AliESDEvent.h"
d631e5e7 43#include "AliESDfriend.h"
5b09c01f 44#include "AliVertexerTracks.h"
45#include "AliESDInputHandler.h"
46#include "AliAnalysisTaskIPInfo.h"
47#include "AliIntSpotEstimator.h"
d631e5e7 48#include "AliMultiplicity.h"
5b09c01f 49
50
51ClassImp(AliAnalysisTaskIPInfo)
52
d631e5e7 53const Char_t* AliAnalysisTaskIPInfo::fEstNames[kNEst] = {"ITSTPC","TPC","SPD"};
54
5b09c01f 55//________________________________________________________________________
56AliAnalysisTaskIPInfo::AliAnalysisTaskIPInfo(const char *name)
57: AliAnalysisTask(name, "IP analysis"),
d631e5e7 58 fESD(0),fESDfriend(0),fOutput(0),fTracks(50)
5b09c01f 59{
60
61 // Define input and output slots here
62 // Input slot #0 works with a TChain
63 DefineInput(0, TChain::Class());
64 // Output slot #0 writes into a TList container
65 DefineOutput(0, TList::Class()); //My private output
66 //
d631e5e7 67 for (int i=0;i<kNEst;i++) fIPEst[i] = 0;
68 //
5b09c01f 69}
70
71//________________________________________________________________________
72AliAnalysisTaskIPInfo::~AliAnalysisTaskIPInfo()
73{
74 // Destructor
75 if (fOutput) {
76 delete fOutput;
77 fOutput = 0;
78 }
79}
d631e5e7 80//________________________________________________________________________
81void AliAnalysisTaskIPInfo::SetIPCenIni(Int_t estID, Double_t x,Double_t y,Double_t z)
82{
83 // set initial estimate of the IP center
84 if (estID<0 || estID>= kNEst) return;
85 fIPCenIni[estID][0] = x;
86 fIPCenIni[estID][1] = y;
87 fIPCenIni[estID][2] = z;
88}
5b09c01f 89
90//________________________________________________________________________
91void AliAnalysisTaskIPInfo::SetOptions(Int_t estID, Bool_t recoVtx,
d631e5e7 92 Double_t outcut,Int_t ntrIP,Int_t nPhiBins,Int_t nestb,
5b09c01f 93 Double_t estmin,Double_t estmax,
94 Int_t ntrBins,Int_t ntMn,Int_t ntMx,
d631e5e7 95 Int_t nPBins,Double_t pmn,Double_t pmx,Bool_t fillNt)
5b09c01f 96{
97 // set options for estimators
98 if (estID<0 || estID>= kNEst) return;
d631e5e7 99 fNTrMinIP[estID] = ntrIP;
100 fRecoVtx[estID] = recoVtx;
5b09c01f 101 fNPhiBins[estID] = nPhiBins;
102 fNEstb[estID] = nestb;
103 fNTrBins[estID] = ntrBins;
104 fNPBins[estID] = nPBins;
105 fNTrMin[estID] = ntMn;
106 fNTrMax[estID] = ntMx;
107 fOutCut[estID] = outcut;
108 fEstMin[estID] = estmin;
109 fEstMax[estID] = estmax;
110 fPMin[estID] = pmn;
111 fPMax[estID] = pmx;
d631e5e7 112 fFillNt[estID] = fillNt;
5b09c01f 113}
114
115//________________________________________________________________________
116void AliAnalysisTaskIPInfo::ConnectInputData(Option_t *)
117{
118 // Connect ESD or AOD here
119 // Called once
120 //
d631e5e7 121 AliInfo("HERE");
5b09c01f 122 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
123 if (!tree) {
124 Printf("ERROR: Could not read chain from input slot 0");
125 }
126 else {
d631e5e7 127 tree->SetBranchAddress("ESDfriend.",&fESDfriend);
5b09c01f 128 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
129 if (!esdH) Printf("ERROR: Could not get ESDInputHandler");
130 else fESD = esdH->GetEvent();
131 }
132 //
133 return;
134}
135
136//________________________________________________________________________
137void AliAnalysisTaskIPInfo::CreateOutputObjects()
138{
139 // Create estimators
140 // Several histograms are more conveniently managed in a TList
141 fOutput = new TList;
142 fOutput->SetOwner();
143 //
d631e5e7 144 Double_t err[3]={0.0400,0.0400,7.5};
5b09c01f 145 //
d631e5e7 146 for (int i=0;i<kNEst;i++) if (fNEstb[i]>1) {
147 TString nm = GetName();
148 TString nmes = fEstNames[i];
149 nm += nmes;
150 fIPEst[i] = new AliIntSpotEstimator(nm.Data(),fOutCut[i],fNTrMinIP[i],fNPhiBins[i],
151 fNEstb[i],fEstMin[i],fEstMax[i],
152 fNTrBins[i],fNTrMin[i],fNTrMax[i],
153 fNPBins[i],fPMin[i],fPMax[i],fFillNt[i]);
154 AliESDVertex *initVertex = new AliESDVertex(fIPCenIni[i],err);
155 fIPEst[i]->GetVertexer()->SetVtxStart(initVertex);
156 delete initVertex;
157 fIPEst[i]->GetVertexer()->SetConstraintOff();
158 fIPEst[i]->GetVertexer()->SetMinClusters(2);
159 fIPEst[i]->SetIPCenIni(fIPCenIni[i]);
160 if (nmes == "TPC") fIPEst[i]->GetVertexer()->SetTPCMode();
161 else fIPEst[i]->GetVertexer()->SetITSMode();
162 //
163 fOutput->Add(fIPEst[i]);
164 if (fIPEst[i]->GetNtuple()) fOutput->Add(fIPEst[i]->GetNtuple());
165 }
5b09c01f 166 //
167 return;
168}
169
170//________________________________________________________________________
171void AliAnalysisTaskIPInfo::Exec(Option_t *)
172{
d631e5e7 173 static TClonesArray tracks("AliExternalTrackParam",50);
174 //
5b09c01f 175 // Main loop
176 // Called for each event
5b09c01f 177 if (!fESD) {
178 Printf("ERROR: fESD not available");
179 return;
180 }
d631e5e7 181 fESD->SetESDfriend(fESDfriend);
5b09c01f 182 //
5b09c01f 183 const AliESDVertex *vtx;
184 UShort_t *trackID;
185 Int_t ntracks;
186 //
d631e5e7 187 for (int ie=0;ie<kNEst;ie++) {
188 //
189 if (!fIPEst[ie]) continue;
190 if (ie==kTPC) {
191 fIPEst[kTPC]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
192 vtx = fRecoVtx[kTPC] ? fIPEst[kTPC]->GetVertexer()->FindPrimaryVertex(fESD) : fESD->GetPrimaryVertexTPC();
193 if (vtx) {
194 ntracks = vtx->GetNIndices();
195 trackID = (UShort_t*)vtx->GetIndices();
196 for (int i=ntracks;i--;) fTracks.Add((TObject*)fESD->GetTrack(trackID[i])->GetTPCInnerParam());
197 fIPEst[kTPC]->ProcessEvent(&fTracks);
198 fTracks.Clear();
199 }
200 }
201 else if (ie==kITSTPC) {
202 fIPEst[kITSTPC]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
203 vtx = fRecoVtx[kITSTPC] ? fIPEst[kITSTPC]->GetVertexer()->FindPrimaryVertex(fESD) : fESD->GetPrimaryVertex();
204 if (vtx) {
205 ntracks = vtx->GetNIndices();
206 trackID = (UShort_t*)vtx->GetIndices();
207 for (int i=ntracks;i--;) fTracks.Add((TObject*)fESD->GetTrack(trackID[i]));
208 fIPEst[kITSTPC]->ProcessEvent(&fTracks);
209 fTracks.Clear();
210 }
211 }
212 else if (ie==kSPD) {
213 fIPEst[kSPD]->GetVertexer()->SetFieldkG( fESD->GetMagneticField() );
214 ntracks = CreateSPDTracklets(tracks);
215 for (int i=ntracks;i--;) fTracks.Add((TObject*)tracks[i]);
216 fIPEst[kSPD]->ProcessEvent(&fTracks);
217 fTracks.Clear();
218 }
5b09c01f 219 }
220 //
5b09c01f 221 PostData(0, fOutput);
222 //
223 return;
224}
225
226//________________________________________________________________________
d631e5e7 227Int_t AliAnalysisTaskIPInfo::CreateSPDTracklets(TClonesArray& tracks)
228{
229 // create traclets from multiplicity class
230 double cv[21] = {
231 25e-4,
232 0 , 25e-4,
233 0 , 0, 40e-2,
234 0 , 0, 0, 1e-2,
235 0 , 0, 0, 0, 1e-2,
236 0 , 0, 0, 0, 0, 1e-2
237 };
238 //
239 double xyzt[3],pxyz[3];
240 int nSPDtracks = 0;
241 tracks.Delete();
242 const AliMultiplicity *mult = fESD->GetMultiplicity();
243 const AliESDVertex *spdVtx = fESD->GetPrimaryVertexSPD();
244 int nTracklets = 0;
245 if (mult && spdVtx && (nTracklets=mult->GetNumberOfTracklets())>2 ) {
246 const Double_t kRLay1=3.9, kRLay2=7.6;
247 double xv = spdVtx->GetXv();
248 double yv = spdVtx->GetYv();
249 double zv = spdVtx->GetZv();
250 for (int i=0;i<nTracklets;i++) { // get cluster coordinates from tracklet
251 double phi1 = mult->GetPhi(i);
252 double tht1 = mult->GetTheta(i);
253 double phi2 = phi1 - mult->GetDeltaPhi(i);
254 double tht2 = tht1 - mult->GetDeltaTheta(i);
255 double cs = TMath::Cos(phi1);
256 double sn = TMath::Sin(tht1);
257 double det = xv*sn+yv*sn;
258 det = det*det + kRLay1*kRLay1;
259 double t = -(xv*cs+yv*sn) + TMath::Sqrt(det);
260 double x1 = cs*t;
261 double y1 = sn*t;
262 double z1 = zv + TMath::Sqrt(x1*x1+y1*y1)/TMath::Tan(tht1);
263 x1 += xv;
264 y1 += yv;
265 //
266 cs = TMath::Cos(phi2);
267 sn = TMath::Sin(tht2);
268 det = xv*sn+yv*sn;
269 det = det*det + kRLay2*kRLay2;
270 t = -(xv*cs+yv*sn) + TMath::Sqrt(det);
271 double dx = cs*t;
272 double dy = sn*t;
273 double dz = zv + TMath::Sqrt(dx*dx+dy*dy)/TMath::Tan(tht2);
274 dx += xv-x1;
275 dy += yv-y1;
276 dz += -z1;
277 double dr = TMath::Sqrt(dx*dx+dy*dy+dz*dz);
278 pxyz[0] = dx/dr; // direction cosines
279 pxyz[1] = dy/dr;
280 pxyz[2] = dz/dr;
281 t = (xv-x1)*pxyz[0] + (yv-y1)*pxyz[1] + (zv-z1)*pxyz[2];
282 xyzt[0] = x1 + t*pxyz[0]; // PCA to vertex
283 xyzt[1] = y1 + t*pxyz[1];
284 xyzt[2] = z1 + t*pxyz[2];
285 //
286 new(tracks[nSPDtracks++]) AliExternalTrackParam(xyzt,pxyz,cv,0);
287 }
288 }
289 return nSPDtracks;
290}
291
292//________________________________________________________________________
5b09c01f 293void AliAnalysisTaskIPInfo::Terminate(Option_t *)
294{
295 // Draw result to the screen
296 // Called once at the end of the query
297 fOutput = dynamic_cast<TList*> (GetOutputData(0));
298 if (!fOutput) {
299 Printf("ERROR: fOutput not available");
300 return;
301 }
302 //
303 return;
304}