]>
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> | |
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 | ||
51 | ClassImp(AliAnalysisTaskIPInfo) | |
52 | ||
d631e5e7 | 53 | const Char_t* AliAnalysisTaskIPInfo::fEstNames[kNEst] = {"ITSTPC","TPC","SPD"}; |
54 | ||
5b09c01f | 55 | //________________________________________________________________________ |
56 | AliAnalysisTaskIPInfo::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 | //________________________________________________________________________ | |
72 | AliAnalysisTaskIPInfo::~AliAnalysisTaskIPInfo() | |
73 | { | |
74 | // Destructor | |
75 | if (fOutput) { | |
76 | delete fOutput; | |
77 | fOutput = 0; | |
78 | } | |
79 | } | |
d631e5e7 | 80 | //________________________________________________________________________ |
81 | void 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 | //________________________________________________________________________ | |
91 | void 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 | //________________________________________________________________________ | |
116 | void 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 | //________________________________________________________________________ | |
137 | void 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 | //________________________________________________________________________ | |
171 | void 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 | ||
d631e5e7 | 226 | //________________________________________________________________________ |
227 | Int_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 | ||
5b09c01f | 292 | //________________________________________________________________________ |
293 | void 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 | } |