AliPhysics  59e0e03 (59e0e03)
AliAnalysisTaskIPInfo.cxx
Go to the documentation of this file.
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 <TNtuple.h>
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"
43 #include "AliESDfriend.h"
44 #include "AliVertexerTracks.h"
45 #include "AliESDInputHandler.h"
46 #include "AliAnalysisTaskIPInfo.h"
47 #include "AliIntSpotEstimator.h"
48 #include "AliMultiplicity.h"
49 
50 
51 ClassImp(AliAnalysisTaskIPInfo)
52 
53 const Char_t* AliAnalysisTaskIPInfo::fEstNames[kNEst] = {"ITSTPC","TPC","SPD"};
54 
55 //________________________________________________________________________
57 : AliAnalysisTask(name, "IP analysis"),
58  fESD(0),fESDfriend(0),fOutput(0),fTracks(50)
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  //
67  for (int i=0;i<kNEst;i++) fIPEst[i] = 0;
68  //
69 }
70 
71 //________________________________________________________________________
73 {
74  // Destructor
75  if (fOutput) {
76  delete fOutput;
77  fOutput = 0;
78  }
79 }
80 //________________________________________________________________________
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 }
89 
90 //________________________________________________________________________
92  Double_t outcut,Int_t ntrIP,Int_t nPhiBins,Int_t nestb,
93  Double_t estmin,Double_t estmax,
94  Int_t ntrBins,Int_t ntMn,Int_t ntMx,
95  Int_t nPBins,Double_t pmn,Double_t pmx,Bool_t fillNt)
96 {
97  // set options for estimators
98  if (estID<0 || estID>= kNEst) return;
99  fNTrMinIP[estID] = ntrIP;
100  fRecoVtx[estID] = recoVtx;
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;
112  fFillNt[estID] = fillNt;
113 }
114 
115 //________________________________________________________________________
117 {
118  // Connect ESD or AOD here
119  // Called once
120  //
121  AliInfo("HERE");
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 {
127  tree->SetBranchAddress("ESDfriend.",&fESDfriend);
128  AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
129  if (!esdH) Printf("ERROR: Could not get ESDInputHandler");
130  else fESD = (AliESDEvent*)esdH->GetEvent();
131  }
132  //
133  return;
134 }
135 
136 //________________________________________________________________________
138 {
139  // Create estimators
140  // Several histograms are more conveniently managed in a TList
141  fOutput = new TList;
142  fOutput->SetOwner();
143  //
144  Double_t err[3]={0.0400,0.0400,7.5};
145  //
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  }
166  //
167  return;
168 }
169 
170 //________________________________________________________________________
172 {
173  static TClonesArray tracks("AliExternalTrackParam",50);
174  //
175  // Main loop
176  // Called for each event
177  if (!fESD) {
178  Printf("ERROR: fESD not available");
179  return;
180  }
181  fESD->SetESDfriend(fESDfriend);
182  //
183  const AliESDVertex *vtx;
184  UShort_t *trackID;
185  Int_t ntracks;
186  //
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());
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]));
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]);
217  fTracks.Clear();
218  }
219  }
220  //
221  PostData(0, fOutput);
222  //
223  return;
224 }
225 
226 //________________________________________________________________________
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->GetX();
248  double yv = spdVtx->GetY();
249  double zv = spdVtx->GetZ();
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 //________________________________________________________________________
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 }
double Double_t
Definition: External.C:58
Int_t fNTrMinIP[kNEst]
request to refit the vertex for given estimator
Double_t fEstMin[kNEst]
outliers cut level
AliVertexerTracks * GetVertexer() const
Double_t fOutCut[kNEst]
max vtx multuplicity
char Char_t
Definition: External.C:18
Int_t fNTrMax[kNEst]
min vtx multuplicity
AliIntSpotEstimator * fIPEst[kNEst]
request to fill ntuple
virtual void Terminate(Option_t *)
void SetOptions(Int_t estID, Bool_t recoVtx=kFALSE, Double_t outcut=1e-4, Int_t ntrIP=2, Int_t nPhiBins=12, Int_t nestb=1000, Double_t estmin=-4e-2, Double_t estmax=6e-2, Int_t ntrBins=10, Int_t ntMn=2, Int_t ntMx=32, Int_t nPBins=14, Double_t pmn=0.2, Double_t pmx=3., Bool_t fillNt=kFALSE)
Int_t CreateSPDTracklets(TClonesArray &tracks)
Double_t fEstMax[kNEst]
lower estimator boundary
Int_t fNPhiBins[kNEst]
min tracks for IP estimator
int Int_t
Definition: External.C:63
Bool_t ProcessEvent(const AliESDEvent *esd, const AliESDVertex *vtx=0)
TObjArray fTracks
list send on output slot 0
Int_t fNPBins[kNEst]
n of vtx.mult. bins
Int_t fNTrMin[kNEst]
n of track P bins
Bool_t fFillNt[kNEst]
initial estimate of IP Center
virtual void Exec(Option_t *option)
void SetIPCenIni(Double_t *xyz)
TList * fOutput
ESD friend object.
Double_t fPMin[kNEst]
upper estimator boundary
Double_t fIPCenIni[kNEst][3]
upper P cut
AliESDEvent * fESD
estimators
Double_t fPMax[kNEst]
lower P cut
void SetIPCenIni(Int_t esdID, Double_t x=0, Double_t y=0, Double_t z=0)
Int_t fNEstb[kNEst]
n bins in phi for IP
AliESDfriend * fESDfriend
ESD object.
virtual void ConnectInputData(Option_t *)
unsigned short UShort_t
Definition: External.C:28
static const Char_t * fEstNames[kNEst]
temporary storage for extracted tracks
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
AliAnalysisTaskIPInfo(const char *name="IPInfo")
Int_t fNTrBins[kNEst]
n of estimator bins