1 ////////////////////////////////////////////////////////////////////////////////
2 // AliAnalysisTaskITSTPCalignment
3 // Runs the relative ITS TPC alignment procedure and TPC vdrift calib
4 // Origin: Mikolaj Krzewicki, mikolaj.krzewicki@cern.ch
5 ////////////////////////////////////////////////////////////////////////////////
13 #include <TObjArray.h>
14 #include <TGraphErrors.h>
16 #include <AliAnalysisTaskSE.h>
17 #include <AliMCEventHandler.h>
18 #include <AliAnalysisManager.h>
19 #include <AliESDEvent.h>
20 #include <AliESDfriend.h>
21 #include <AliMCEvent.h>
23 #include <AliExternalTrackParam.h>
24 #include <AliESDtrack.h>
25 #include <AliESDfriendTrack.h>
26 #include <AliESDInputHandler.h>
27 #include "AliRelAlignerKalman.h"
28 #include "AliRelAlignerKalmanArray.h"
29 #include "AliAnalysisTaskITSTPCalignment.h"
32 ClassImp(AliAnalysisTaskITSTPCalignment)
34 //________________________________________________________________________
35 AliAnalysisTaskITSTPCalignment::AliAnalysisTaskITSTPCalignment():
42 fFillDebugTree(kFALSE),
50 fRejectOutliers(kTRUE),
52 fRejectOutliersSigma2Median(kFALSE),
53 fOutRejSigma2Median(5.),
54 fOutRejSigmaOnMerge(10.),
55 fUseITSoutGlobalTrack(kFALSE),
56 fUseITSoutITSSAtrack(kTRUE)
59 // Define input and output slots here
60 // Input slot #0 works with a TChain
61 //DefineInput(0, TChain::Class());
62 //DefineOutput(0, TTree::Class());
63 //DefineOutput(1, TList::Class());
64 //DefineOutput(2, AliRelAlignerKalmanArray::Class());
67 //________________________________________________________________________
68 AliAnalysisTaskITSTPCalignment::AliAnalysisTaskITSTPCalignment(const char *name):
69 AliAnalysisTaskSE(name),
75 fFillDebugTree(kFALSE),
83 fRejectOutliers(kTRUE),
85 fRejectOutliersSigma2Median(kFALSE),
86 fOutRejSigma2Median(5.),
87 fOutRejSigmaOnMerge(10.),
88 fUseITSoutGlobalTrack(kFALSE),
89 fUseITSoutITSSAtrack(kTRUE)
93 // Define input and output slots here
94 // Input slot #0 is a TChain
95 // Output slot #0 is a TTree
96 DefineOutput(1, TList::Class());
97 DefineOutput(2, AliRelAlignerKalmanArray::Class());
98 DefineOutput(3, AliRelAlignerKalmanArray::Class());
101 //________________________________________________________________________
102 Bool_t AliAnalysisTaskITSTPCalignment::UserNotify()
104 //ON INPUT FILE/TREE CHANGE
109 if (fAligner->GetNUpdates()>0) fDebugTree->Fill();
115 //________________________________________________________________________
116 void AliAnalysisTaskITSTPCalignment::UserCreateOutputObjects()
118 // Create output objects
120 fArrayITSglobal = new AliRelAlignerKalmanArray();
121 fArrayITSglobal->SetName("outputArrayITSglobal");
122 fArrayITSglobal->SetupArray(fT0,fTend,fSlotWidth);
123 fArrayITSglobal->SetOutRejSigmaOnMerge(fOutRejSigmaOnMerge);
124 //set up the template
125 AliRelAlignerKalman* templ = fArrayITSglobal->GetAlignerTemplate();
126 templ->SetRejectOutliers(fRejectOutliers);
127 templ->SetOutRejSigma(fOutRejSigma);
128 templ->SetRejectOutliersSigma2Median(fRejectOutliersSigma2Median);
129 templ->SetOutRejSigma2Median(fOutRejSigma2Median);
130 fArrayITSsa = new AliRelAlignerKalmanArray();
131 fArrayITSsa->SetName("outputArrayITSsa");
132 fArrayITSsa->SetupArray(fT0,fTend,fSlotWidth);
133 fArrayITSsa->SetOutRejSigmaOnMerge(fOutRejSigmaOnMerge);
134 //set up the template
135 templ = fArrayITSsa->GetAlignerTemplate();
136 templ->SetRejectOutliers(fRejectOutliers);
137 templ->SetOutRejSigma(fOutRejSigma);
138 templ->SetRejectOutliersSigma2Median(fRejectOutliersSigma2Median);
139 templ->SetOutRejSigma2Median(fOutRejSigma2Median);
143 TH2F* pZYAResidualsHistBpos = new TH2F("fZYAResidualsHistBpos","z-r\\phi residuals side A (z>0), Bpos", 100, -4., 4., 100, -1.5, 1.5 );
144 pZYAResidualsHistBpos->SetXTitle("\\deltaz [cm]");
145 pZYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
146 TH2F* pZYCResidualsHistBpos = new TH2F("fZYCResidualsHistBpos","z-r\\phi residuals side C (z<0), Bpos", 100, -4., 4., 100, -1.5, 1.5 );
147 pZYCResidualsHistBpos->SetXTitle("\\deltaz [cm]");
148 pZYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
149 TH2F* pLPAResidualsHistBpos = new TH2F("fLPAResidualsHistBpos", "sin(\\phi) tan(\\lambda) residuals side A (z>0), Bpos", 100, -.05, 0.05, 100, -0.05, 0.05 );
150 pLPAResidualsHistBpos->SetXTitle("\\deltasin(\\phi)");
151 pLPAResidualsHistBpos->SetYTitle("\\deltatan(\\lambda)");
152 TH2F* pLPCResidualsHistBpos = new TH2F("fLPCResidualsHistBpos", "sin(\\phi) tan(\\lambda) residuals side C (z<0), Bpos", 100, -.05, 0.05, 100, -0.05, 0.05 );
153 pLPCResidualsHistBpos->SetXTitle("\\deltasin(\\phi)");
154 pLPCResidualsHistBpos->SetYTitle("\\deltatan(\\lambda)");
155 TH2F* pPhiYAResidualsHistBpos = new TH2F("fPhiYAResidualsHistBpos","\\phi-\\deltar\\phi side A (z>0), Bpos",36,0.0,6.3,100,-1.5,1.5);
156 pPhiYAResidualsHistBpos->SetXTitle("\\phi [rad]");
157 pPhiYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
158 TH2F* pPhiYCResidualsHistBpos = new TH2F("fPhiYCResidualsHistBpos","\\phi-\\deltar\\phi side C (z<0), Bpos",36,0.0,6.3,100,-1.5,1.5);
159 pPhiYCResidualsHistBpos->SetXTitle("\\phi [rad]");
160 pPhiYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
161 TH2F* pPhiZAResidualsHistBpos = new TH2F("fPhiZAResidualsHistBpos","\\phi-\\deltaz side A (z>0), Bpos",36,0.0,6.3,100,-3.0,3.0);
162 pPhiZAResidualsHistBpos->SetXTitle("\\phi [rad]");
163 pPhiZAResidualsHistBpos->SetYTitle("\\deltaz [cm]");
164 TH2F* pPhiZCResidualsHistBpos = new TH2F("fPhiZCResidualsHistBpos","\\phi-\\deltaz side C (z<0), Bpos",36,0.0,6.3,100,-3.0,3.0);
165 pPhiZCResidualsHistBpos->SetXTitle("\\phi [rad]");
166 pPhiZCResidualsHistBpos->SetYTitle("\\deltaz [cm]");
167 TH2F* pPtYAResidualsHistBpos = new TH2F("fPtYAResidualsHistBpos","Pt-\\deltar\\phi side A (z>0), Bpos",20,.3,8.,80,-1.5,1.5);
168 pPtYAResidualsHistBpos->SetXTitle("Pt [rad]");
169 pPtYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
170 TH2F* pPtYCResidualsHistBpos = new TH2F("fPtYCResidualsHistBpos","Pt-\\deltar\\phi side C (z<0), Bpos",20,.3,8.,80,-1.5,1.5);
171 pPtYCResidualsHistBpos->SetXTitle("Pt [rad]");
172 pPtYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
173 TH2F* pPtZAResidualsHistBpos = new TH2F("fPtZAResidualsHistBpos","Pt-\\deltaz side A (z>0), Bpos",20,.3,8.,80,-3.0,3.0);
174 pPtZAResidualsHistBpos->SetXTitle("Pt");
175 pPtZAResidualsHistBpos->SetYTitle("\\deltaz [cm]");
176 TH2F* pPtZCResidualsHistBpos = new TH2F("fPtZCResidualsHistBpos","Pt-\\deltaz side C (z<0), Bpos",20,.3,8.,80,-3.0,3.0);
177 pPtZCResidualsHistBpos->SetXTitle("Pt");
178 pPtZCResidualsHistBpos->SetYTitle("\\deltaz [cm]");
179 TH2F* pLowPtYAResidualsHistBpos = new TH2F("fLowPtYAResidualsHistBpos","Pt-\\deltar\\phi side A (z>0), Bpos",100,0.0,1.5,80,-1.5,1.5);
180 pLowPtYAResidualsHistBpos->SetXTitle("Pt [rad]");
181 pLowPtYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
182 TH2F* pLowPtYCResidualsHistBpos = new TH2F("fLowPtYCResidualsHistBpos","Pt-\\deltar\\phi side C (z<0), Bpos",100,0.0,1.5,80,-1.5,1.5);
183 pLowPtYCResidualsHistBpos->SetXTitle("Pt [rad]");
184 pLowPtYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
185 TH2F* pLowPtZAResidualsHistBpos = new TH2F("fLowPtZAResidualsHistBpos","Pt-\\deltaz side A (z>0), Bpos",100,0.0,1.5,80,-3.0,3.0);
186 pLowPtZAResidualsHistBpos->SetXTitle("Pt");
187 pLowPtZAResidualsHistBpos->SetYTitle("\\deltaz [cm]");
188 TH2F* pLowPtZCResidualsHistBpos = new TH2F("fLowPtZCResidualsHistBpos","Pt-\\deltaz side C (z<0), Bpos",100,0.0,1.5,80,-3.0,3.0);
189 pLowPtZCResidualsHistBpos->SetXTitle("Pt");
190 pLowPtZCResidualsHistBpos->SetYTitle("\\deltaz [cm]");
191 TList* listBpos = new TList();
192 fList->Add(listBpos); //0
193 listBpos->Add(pZYAResidualsHistBpos);
194 listBpos->Add(pZYCResidualsHistBpos);
195 listBpos->Add(pLPAResidualsHistBpos);
196 listBpos->Add(pLPCResidualsHistBpos);
197 listBpos->Add(pPhiYAResidualsHistBpos);
198 listBpos->Add(pPhiYCResidualsHistBpos);
199 listBpos->Add(pPhiZAResidualsHistBpos);
200 listBpos->Add(pPhiZCResidualsHistBpos);
201 listBpos->Add(pPtYAResidualsHistBpos);
202 listBpos->Add(pPtYCResidualsHistBpos);
203 listBpos->Add(pPtZAResidualsHistBpos);
204 listBpos->Add(pPtZCResidualsHistBpos);
205 listBpos->Add(pLowPtYAResidualsHistBpos);
206 listBpos->Add(pLowPtYCResidualsHistBpos);
207 listBpos->Add(pLowPtZAResidualsHistBpos);
208 listBpos->Add(pLowPtZCResidualsHistBpos);
210 TH2F* pZYAResidualsHistBneg = new TH2F("fZYAResidualsHistBneg","z-r\\phi residuals side A (z>0), Bneg", 100, -4., 4., 100, -1.5, 1.5 );
211 pZYAResidualsHistBneg->SetXTitle("\\deltaz [cm]");
212 pZYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
213 TH2F* pZYCResidualsHistBneg = new TH2F("fZYCResidualsHistBneg","z-r\\phi residuals side C (z<0), Bneg", 100, -4., 4., 100, -1.5, 1.5 );
214 pZYCResidualsHistBneg->SetXTitle("\\deltaz [cm]");
215 pZYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
216 TH2F* pLPAResidualsHistBneg = new TH2F("fLPAResidualsHistBneg", "sin(\\phi) tan(\\lambda) residuals side A (z>0), Bneg", 100, -.05, 0.05, 100, -0.05, 0.05 );
217 pLPAResidualsHistBneg->SetXTitle("\\deltasin(\\phi)");
218 pLPAResidualsHistBneg->SetYTitle("\\deltatan(\\lambda)");
219 TH2F* pLPCResidualsHistBneg = new TH2F("fLPCResidualsHistBneg", "sin(\\phi) tan(\\lambda) residuals side C (z<0), Bneg", 100, -.05, 0.05, 100, -0.05, 0.05 );
220 pLPCResidualsHistBneg->SetXTitle("\\deltasin(\\phi)");
221 pLPCResidualsHistBneg->SetYTitle("\\deltatan(\\lambda)");
222 TH2F* pPhiYAResidualsHistBneg = new TH2F("fPhiYAResidualsHistBneg","\\phi-\\deltar\\phi side A (z>0), Bneg",36,0.0,6.3,100,-1.5,1.5);
223 pPhiYAResidualsHistBneg->SetXTitle("\\phi [rad]");
224 pPhiYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
225 TH2F* pPhiYCResidualsHistBneg = new TH2F("fPhiYCResidualsHistBneg","\\phi-\\deltar\\phi side C (z<0), Bneg",36,0.0,6.3,100,-1.5,1.5);
226 pPhiYCResidualsHistBneg->SetXTitle("\\phi [rad]");
227 pPhiYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
228 TH2F* pPhiZAResidualsHistBneg = new TH2F("fPhiZAResidualsHistBneg","\\phi-\\deltaz side A (z>0), Bneg",36,0.0,6.3,100,-3.0,3.0);
229 pPhiZAResidualsHistBneg->SetXTitle("\\phi [rad]");
230 pPhiZAResidualsHistBneg->SetYTitle("\\deltaz [cm]");
231 TH2F* pPhiZCResidualsHistBneg = new TH2F("fPhiZCResidualsHistBneg","\\Pt-\\deltaz side C (z<0), Bneg",36,0.0,6.3,100,-3.0,3.0);
232 pPhiZCResidualsHistBneg->SetXTitle("\\phi [rad]");
233 pPhiZCResidualsHistBneg->SetYTitle("\\deltaz [cm]");
234 TH2F* pPtYAResidualsHistBneg = new TH2F("fPtYAResidualsHistBneg","Pt-\\deltar\\phi side A (z>0), Bneg",20,.3,8.,80,-1.5,1.5);
235 pPtYAResidualsHistBneg->SetXTitle("Pt [rad]");
236 pPtYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
237 TH2F* pPtYCResidualsHistBneg = new TH2F("fPtYCResidualsHistBneg","Pt-\\deltar\\phi side C (z<0), Bneg",20,.3,8.,80,-1.5,1.5);
238 pPtYCResidualsHistBneg->SetXTitle("Pt [rad]");
239 pPtYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
240 TH2F* pPtZAResidualsHistBneg = new TH2F("fPtZAResidualsHistBneg","Pt-\\deltaz side A (z>0), Bneg",20,.3,8.,80,-3.0,3.0);
241 pPtZAResidualsHistBneg->SetXTitle("Pt");
242 pPtZAResidualsHistBneg->SetYTitle("\\deltaz [cm]");
243 TH2F* pPtZCResidualsHistBneg = new TH2F("fPtZCResidualsHistBneg","Pt-\\deltaz side C (z<0), Bneg",20,.3,8.,80,-3.0,3.0);
244 pPtZCResidualsHistBneg->SetXTitle("Pt");
245 pPtZCResidualsHistBneg->SetYTitle("\\deltaz [cm]");
246 TH2F* pLowPtYAResidualsHistBneg = new TH2F("fLowPtYAResidualsHistBneg","Pt-\\deltar\\phi side A (z>0), Bneg",100,0.0,1.5,80,-1.5,1.5);
247 pLowPtYAResidualsHistBneg->SetXTitle("Pt [rad]");
248 pLowPtYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
249 TH2F* pLowPtYCResidualsHistBneg = new TH2F("fLowPtYCResidualsHistBneg","Pt-\\deltar\\phi side C (z<0), Bneg",100,0.0,1.5,80,-1.5,1.5);
250 pLowPtYCResidualsHistBneg->SetXTitle("Pt [rad]");
251 pLowPtYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
252 TH2F* pLowPtZAResidualsHistBneg = new TH2F("fLowPtZAResidualsHistBneg","Pt-\\deltaz side A (z>0), Bneg",100,0.0,1.5,80,-3.0,3.0);
253 pLowPtZAResidualsHistBneg->SetXTitle("Pt");
254 pLowPtZAResidualsHistBneg->SetYTitle("\\deltaz [cm]");
255 TH2F* pLowPtZCResidualsHistBneg = new TH2F("fLowPtZCResidualsHistBneg","Pt-\\deltaz side C (z<0), Bneg",100,0.0,1.5,80,-3.0,3.0);
256 pLowPtZCResidualsHistBneg->SetXTitle("Pt");
257 pLowPtZCResidualsHistBneg->SetYTitle("\\deltaz [cm]");
258 TList* listBneg = new TList();
259 fList->Add(listBneg); //1
260 listBneg->Add(pZYAResidualsHistBneg);
261 listBneg->Add(pZYCResidualsHistBneg);
262 listBneg->Add(pLPAResidualsHistBneg);
263 listBneg->Add(pLPCResidualsHistBneg);
264 listBneg->Add(pPhiYAResidualsHistBneg);
265 listBneg->Add(pPhiYCResidualsHistBneg);
266 listBneg->Add(pPhiZAResidualsHistBneg);
267 listBneg->Add(pPhiZCResidualsHistBneg);
268 listBneg->Add(pPtYAResidualsHistBneg);
269 listBneg->Add(pPtYCResidualsHistBneg);
270 listBneg->Add(pPtZAResidualsHistBneg);
271 listBneg->Add(pPtZCResidualsHistBneg);
272 listBneg->Add(pLowPtYAResidualsHistBneg);
273 listBneg->Add(pLowPtYCResidualsHistBneg);
274 listBneg->Add(pLowPtZAResidualsHistBneg);
275 listBneg->Add(pLowPtZCResidualsHistBneg);
277 TH2F* pZYAResidualsHistBnil = new TH2F("fZYAResidualsHistBnil","z-r\\phi residuals side A (z>0), Bnil", 100, -3., 3., 100, -1.5, 1.5 );
278 pZYAResidualsHistBnil->SetXTitle("\\deltaz [cm]");
279 pZYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
280 TH2F* pZYCResidualsHistBnil = new TH2F("fZYCResidualsHistBnil","z-r\\phi residuals side C (z<0), Bnil", 100, -3., 3., 100, -1.5, 1.5 );
281 pZYCResidualsHistBnil->SetXTitle("\\deltaz [cm]");
282 pZYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
283 TH2F* pLPAResidualsHistBnil = new TH2F("fLPAResidualsHistBnil", "sin(\\phi) tan(\\lambda) residuals side A (z>0), Bnil", 100, -.05, 0.05, 100, -0.05, 0.05 );
284 pLPAResidualsHistBnil->SetXTitle("\\deltasin(\\phi)");
285 pLPAResidualsHistBnil->SetYTitle("\\deltatan(\\lambda)");
286 TH2F* pLPCResidualsHistBnil = new TH2F("fLPCResidualsHistBnil", "sin(\\phi) tan(\\lambda) residuals side C (z<0), Bnil", 100, -.05, 0.05, 100, -0.05, 0.05 );
287 pLPCResidualsHistBnil->SetXTitle("\\deltasin(\\phi)");
288 pLPCResidualsHistBnil->SetYTitle("\\deltatan(\\lambda)");
289 TH2F* pPhiYAResidualsHistBnil = new TH2F("fPhiYAResidualsHistBnil","\\phi-\\deltar\\phi side A (z>0), Bnil",36,0.0,6.3,100,-1.5,1.5);
290 pPhiYAResidualsHistBnil->SetXTitle("\\phi [rad]");
291 pPhiYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
292 TH2F* pPhiYCResidualsHistBnil = new TH2F("fPhiYCResidualsHistBnil","\\phi-\\deltar\\phi side C (z<0), Bnil",36,0.0,6.3,100,-1.5,1.5);
293 pPhiYCResidualsHistBnil->SetXTitle("\\phi [rad]");
294 pPhiYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
295 TH2F* pPhiZAResidualsHistBnil = new TH2F("fPhiZAResidualsHistBnil","\\phi-\\deltaz side A (z>0), Bnil",36,0.0,6.3,100,-3.0,3.0);
296 pPhiZAResidualsHistBnil->SetXTitle("\\phi [rad]");
297 pPhiZAResidualsHistBnil->SetYTitle("\\deltaz [cm]");
298 TH2F* pPhiZCResidualsHistBnil = new TH2F("fPhiZCResidualsHistBnil","\\Pt-\\deltaz side C (z<0), Bnil",36,0.0,6.3,100,-3.0,3.0);
299 pPhiZCResidualsHistBnil->SetXTitle("\\phi [rad]");
300 pPhiZCResidualsHistBnil->SetYTitle("\\deltaz [cm]");
301 TH2F* pPtYAResidualsHistBnil = new TH2F("fPtYAResidualsHistBnil","Pt-\\deltar\\phi side A (z>0), Bnil",20,.3,8.,80,-1.5,1.5);
302 pPtYAResidualsHistBnil->SetXTitle("Pt [rad]");
303 pPtYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
304 TH2F* pPtYCResidualsHistBnil = new TH2F("fPtYCResidualsHistBnil","Pt-\\deltar\\phi side C (z<0), Bnil",20,.3,8.,80,-1.5,1.5);
305 pPtYCResidualsHistBnil->SetXTitle("Pt [rad]");
306 pPtYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
307 TH2F* pPtZAResidualsHistBnil = new TH2F("fPtZAResidualsHistBnil","Pt-\\deltaz side A (z>0), Bnil",20,.3,8.,80,-3.0,3.0);
308 pPtZAResidualsHistBnil->SetXTitle("Pt");
309 pPtZAResidualsHistBnil->SetYTitle("\\deltaz [cm]");
310 TH2F* pPtZCResidualsHistBnil = new TH2F("fPtZCResidualsHistBnil","Pt-\\deltaz side C (z<0), Bnil",20,.3,8.,80,-3.0,3.0);
311 pPtZCResidualsHistBnil->SetXTitle("Pt");
312 pPtZCResidualsHistBnil->SetYTitle("\\deltaz [cm]");
313 TH2F* pLowPtYAResidualsHistBnil = new TH2F("fLowPtYAResidualsHistBnil","Pt-\\deltar\\phi side A (z>0), Bnil",100,0.0,1.5,80,-1.5,1.5);
314 pLowPtYAResidualsHistBnil->SetXTitle("Pt [rad]");
315 pLowPtYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
316 TH2F* pLowPtYCResidualsHistBnil = new TH2F("fLowPtYCResidualsHistBnil","Pt-\\deltar\\phi side C (z<0), Bnil",100,0.0,1.5,80,-1.5,1.5);
317 pLowPtYCResidualsHistBnil->SetXTitle("Pt [rad]");
318 pLowPtYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
319 TH2F* pLowPtZAResidualsHistBnil = new TH2F("fLowPtZAResidualsHistBnil","Pt-\\deltaz side A (z>0), Bnil",100,0.0,1.5,80,-3.0,3.0);
320 pLowPtZAResidualsHistBnil->SetXTitle("Pt");
321 pLowPtZAResidualsHistBnil->SetYTitle("\\deltaz [cm]");
322 TH2F* pLowPtZCResidualsHistBnil = new TH2F("fLowPtZCResidualsHistBnil","Pt-\\deltaz side C (z<0), Bnil",100,0.0,1.5,80,-3.0,3.0);
323 pLowPtZCResidualsHistBnil->SetXTitle("Pt");
324 pLowPtZCResidualsHistBnil->SetYTitle("\\deltaz [cm]");
325 TList* listBnil = new TList();
326 fList->Add(listBnil); //2
327 listBnil->Add(pZYAResidualsHistBnil);
328 listBnil->Add(pZYCResidualsHistBnil);
329 listBnil->Add(pLPAResidualsHistBnil);
330 listBnil->Add(pLPCResidualsHistBnil);
331 listBnil->Add(pPhiYAResidualsHistBnil);
332 listBnil->Add(pPhiYCResidualsHistBnil);
333 listBnil->Add(pPhiZAResidualsHistBnil);
334 listBnil->Add(pPhiZCResidualsHistBnil);
335 listBnil->Add(pPtYAResidualsHistBnil);
336 listBnil->Add(pPtYCResidualsHistBnil);
337 listBnil->Add(pPtZAResidualsHistBnil);
338 listBnil->Add(pPtZCResidualsHistBnil);
339 listBnil->Add(pLowPtYAResidualsHistBnil);
340 listBnil->Add(pLowPtYCResidualsHistBnil);
341 listBnil->Add(pLowPtZAResidualsHistBnil);
342 listBnil->Add(pLowPtZCResidualsHistBnil);
343 TH1F* pNmatchingEff=new TH1F("pNmatchingEff","matching efficiency",50,0.,1.);
344 fList->Add(pNmatchingEff); //3
346 TH1I* pChecks = new TH1I("pChecks","various checks",10,0,10);
347 pChecks->GetXaxis()->SetBinLabel(1+kNoESD,"no ESD");
348 pChecks->GetXaxis()->SetBinLabel(1+kNoESDfriend,"no ESDfriend");
349 pChecks->GetXaxis()->SetBinLabel(1+kNoFriendTrack,"no friendtrack");
350 pChecks->GetXaxis()->SetBinLabel(1+kNoITSoutParams,"no params");
351 pChecks->GetXaxis()->SetBinLabel(1+kESDfriend,"ESD friends");
352 pChecks->GetXaxis()->SetBinLabel(1+kFriendsSkipBit,"friends+skipbit");
353 pChecks->GetXaxis()->SetBinLabel(1+kFriendTrack,"friendtrack");
354 pChecks->GetXaxis()->SetBinLabel(1+kITSoutParams,"params");
355 fList->Add(pChecks); //4
357 fAligner = new AliRelAlignerKalman();
358 fAligner->SetRejectOutliers(fRejectOutliers);
359 fAligner->SetOutRejSigma(fOutRejSigma);
360 fAligner->SetRejectOutliersSigma2Median(fRejectOutliersSigma2Median);
361 fAligner->SetOutRejSigma2Median(fOutRejSigma2Median);
362 fList->Add(fAligner);
364 fDebugTree = new TTree("debugTree","tree with debug info");
365 fDebugTree->Branch("aligner","AliRelAlignerKalman",&fAligner);
368 //________________________________________________________________________
369 void AliAnalysisTaskITSTPCalignment::UserExec(Option_t *)
372 // Called for each event
374 TH1I* pChecks = static_cast<TH1I*>(fList->At(4));
375 //check for ESD and friends
376 AliESDEvent* event = dynamic_cast<AliESDEvent*>(InputEvent());
377 AliESDfriend* eventFriend = dynamic_cast<AliESDfriend*>(ESDfriend());
381 pChecks->Fill(kNoESD);
387 AliError("no ESD friend");
388 pChecks->Fill(kNoESDfriend);
393 pChecks->Fill(kESDfriend);
396 if (eventFriend->TestSkipBit()) pChecks->Fill(kFriendsSkipBit);
398 //Update the parmeters
399 AnalyzeESDevent(event);
402 PostData(0, fDebugTree);
404 PostData(2, fArrayITSglobal);
405 PostData(3, fArrayITSsa);
408 //________________________________________________________________________
409 void AliAnalysisTaskITSTPCalignment::AnalyzeESDevent(AliESDEvent* event)
411 //analyze an ESD event with track matching
412 Int_t ntracks = event->GetNumberOfTracks();
413 if (ntracks==0) return;
415 if (fUseITSoutGlobalTrack)
417 AliRelAlignerKalman* alignerGlobal = fArrayITSglobal->GetAligner(event);
419 alignerGlobal->AddESDevent(event);
422 if (fUseITSoutITSSAtrack)
424 //get the aligner for the event
425 //do it with matching
426 TObjArray arrayParamITS(ntracks);
427 TObjArray arrayParamTPC(ntracks);
429 Int_t n = FindMatchingTracks(arrayParamITS, arrayParamTPC, event);
430 AliInfo(Form("n matched: %i\n",n));
432 TH1F* nMatchingEff=static_cast<TH1F*>(fList->At(3));
433 Float_t ratio = (float)n/(float)ntracks;
434 nMatchingEff->Fill(ratio);
436 AliRelAlignerKalman* alignerSA = NULL;
437 if (n>0) alignerSA = fArrayITSsa->GetAligner(event);
439 for (Int_t i=0;i<n;i++)
441 AliExternalTrackParam* paramsITS=static_cast<AliExternalTrackParam*>(arrayParamITS[i]);
442 AliExternalTrackParam* paramsTPC=static_cast<AliExternalTrackParam*>(arrayParamTPC[i]);
444 if (!(paramsITS&¶msTPC)) continue;
447 if (fDoQA) DoQA(paramsITS,paramsTPC);
450 if (fAligner->AddTrackParams(paramsITS, paramsTPC))
452 fAligner->SetRunNumber(event->GetRunNumber());
453 fAligner->SetMagField(event->GetMagneticField());
454 fAligner->SetTimeStamp(event->GetTimeStamp());
458 if (alignerSA) alignerSA->AddTrackParams(paramsITS, paramsTPC);
460 arrayParamITS.Delete();
461 arrayParamTPC.Delete();
465 //________________________________________________________________________
466 void AliAnalysisTaskITSTPCalignment::Terminate(Option_t *)
468 // Called once at the end of the query
471 //________________________________________________________________________
472 Int_t AliAnalysisTaskITSTPCalignment::FindMatchingTracks(TObjArray& arrITS, TObjArray& arrTPC, AliESDEvent* pESD)
474 //find matching tracks and return tobjarrays with the params
475 //fUniqueID of param object contains tracknumber in event.
477 TH1I* pChecks = static_cast<TH1I*>(fList->At(4));
478 Int_t ntracks = pESD->GetNumberOfTracks();
479 Double_t magfield = pESD->GetMagneticField();
481 Int_t* matchedArr = new Int_t[ntracks]; //storage for index of ITS track for which a match was found
482 for (Int_t i=0;i<ntracks;i++)
488 for (Int_t i=0; i<ntracks; i++)
490 //get track1 and friend
491 AliESDtrack* track1 = pESD->GetTrack(i);
492 if (!track1) continue;
494 if (track1->GetNcls(0) < fMinPointsVol1) continue;
496 if (!( ( track1->IsOn(AliESDtrack::kITSrefit)) &&
497 (!track1->IsOn(AliESDtrack::kTPCin)) )) continue;
499 const AliESDfriendTrack* constfriendtrack1 = track1->GetFriendTrack();
500 if (!constfriendtrack1) {AliInfo(Form("no friendTrack1\n"));pChecks->Fill(kNoFriendTrack);continue;}
501 pChecks->Fill(kFriendTrack);
502 AliESDfriendTrack friendtrack1(*constfriendtrack1);
504 if (!friendtrack1.GetITSOut()) {pChecks->Fill(kNoITSoutParams);continue;}
505 pChecks->Fill(kITSoutParams);
506 AliExternalTrackParam params1(*(friendtrack1.GetITSOut()));
508 Double_t bestd = 1000.; //best distance
509 Bool_t newi = kTRUE; //whether we start with a new i
510 for (Int_t j=0; j<ntracks; j++)
512 if (matchedArr[j]>0 && matchedArr[j]!=i) continue; //already matched, everything tried
513 //get track2 and friend
514 AliESDtrack* track2 = pESD->GetTrack(j);
515 if (!track2) continue;
516 if (track1==track2) continue;
517 if (!(track2->IsOn(AliESDtrack::kTPCin))) continue; //must be TPC
518 if (!(track2->IsOn(AliESDtrack::kITSout))) continue; //must have ITS
520 //if (track2->GetNcls(0) != track1->GetNcls(0)) continue;
521 //if (track2->GetITSClusterMap() != track1->GetITSClusterMap()) continue;
522 if (track2->GetNcls(1) < fMinPointsVol2) continue; //min 80 clusters in TPC
523 if (track2->GetTgl() > 1.) continue; //acceptance
524 //cut crossing tracks
525 if (!track2->GetInnerParam()) continue;
526 if (!track2->GetOuterParam()) continue;
527 if (track2->GetOuterParam()->GetZ()*track2->GetInnerParam()->GetZ()<0) continue;
528 if (track2->GetInnerParam()->GetX()>90) continue;
529 if (TMath::Abs(track2->GetInnerParam()->GetZ())<5.) continue; //too close to membrane?
531 AliExternalTrackParam params2(*(track2->GetInnerParam()));
533 //bring to same reference plane
534 if (!params2.Rotate(params1.GetAlpha())) continue;
535 if (!params2.PropagateTo(params1.GetX(), magfield)) continue;
537 //pt cut, only for data with magfield
538 if ((params2.Pt()<fMinPt)&&(TMath::Abs(magfield)>1.)) continue;
542 // AliStack* stack = fMC->Stack();
543 // Int_t label = track2->GetLabel();
544 // if (label<0) continue;
545 // TParticle* particle = stack->Particle(label);
546 // if (particle->Pt()<fMinPt) continue;
549 const Double32_t* p1 = params1.GetParameter();
550 const Double32_t* p2 = params2.GetParameter();
553 Double_t dy = TMath::Abs(p2[0]-p1[0]);
554 Double_t dz = TMath::Abs(p2[1]-p1[1]);
555 Double_t dphi = TMath::Abs(p2[2]-p1[2]);
556 Double_t dlam = TMath::Abs(p2[3]-p1[3]);
557 if (dy > 2.0) continue;
558 if (dz > 10.0) continue;
559 if (dphi > 0.1 ) continue;
560 if (dlam > 0.1 ) continue;
563 Double_t d = TMath::Sqrt(dy*dy+dz*dz+dphi*dphi+dlam*dlam);
564 if ( d >= bestd) continue;
566 matchedArr[j]=i; //j-th track matches i-th (ITS) track
567 if (newi) iMatched++; newi=kFALSE; //increment at most once per i
568 params1.SetUniqueID(i); //store tracknummer
569 params2.SetUniqueID(j);
570 if (arrITS[iMatched] && arrTPC[iMatched])
572 *(arrITS[iMatched]) = params1;
573 *(arrTPC[iMatched]) = params2;
577 arrITS[iMatched] = new AliExternalTrackParam(params1);
578 arrTPC[iMatched] = new AliExternalTrackParam(params2);
582 delete [] matchedArr;
586 //________________________________________________________________________
587 void AliAnalysisTaskITSTPCalignment::DoQA(AliExternalTrackParam* paramsITS,
588 AliExternalTrackParam* paramsTPC)
590 //fill qa histograms in a given list, per field direction (5,-5,0)
591 Float_t resy = paramsTPC->GetY() - paramsITS->GetY();
592 Float_t resz = paramsTPC->GetZ() - paramsITS->GetZ();
593 Float_t ressnp = paramsTPC->GetSnp() - paramsITS->GetSnp();
594 Float_t restgl = paramsTPC->GetTgl() - paramsITS->GetTgl();
597 Double_t field = fInputEvent->GetMagneticField();
598 if (field >= 1.) pList = static_cast<TList*>(fList->At(0));
599 if (field <= -1.) pList = static_cast<TList*>(fList->At(1));
600 if (field < 1. && field > -1.) pList = static_cast<TList*>(fList->At(2));
603 TH2F* pZYAResidualsHist = static_cast<TH2F*>(pList->At(0));
604 TH2F* pZYCResidualsHist = static_cast<TH2F*>(pList->At(1));
605 TH2F* pLPAResidualsHist = static_cast<TH2F*>(pList->At(2));
606 TH2F* pLPCResidualsHist = static_cast<TH2F*>(pList->At(3));
607 TH2F* pPhiYAResidualsHist = static_cast<TH2F*>(pList->At(4));
608 TH2F* pPhiYCResidualsHist = static_cast<TH2F*>(pList->At(5));
609 TH2F* pPhiZAResidualsHist = static_cast<TH2F*>(pList->At(6));
610 TH2F* pPhiZCResidualsHist = static_cast<TH2F*>(pList->At(7));
611 TH2F* pPtYAResidualsHist = static_cast<TH2F*>(pList->At(8));
612 TH2F* pPtYCResidualsHist = static_cast<TH2F*>(pList->At(9));
613 TH2F* pPtZAResidualsHist = static_cast<TH2F*>(pList->At(10));
614 TH2F* pPtZCResidualsHist = static_cast<TH2F*>(pList->At(11));
615 TH2F* pLowPtYAResidualsHist = static_cast<TH2F*>(pList->At(12));
616 TH2F* pLowPtYCResidualsHist = static_cast<TH2F*>(pList->At(13));
617 TH2F* pLowPtZAResidualsHist = static_cast<TH2F*>(pList->At(14));
618 TH2F* pLowPtZCResidualsHist = static_cast<TH2F*>(pList->At(15));
620 if (paramsITS->GetZ() > 0.0)
622 pZYAResidualsHist->Fill(resz,resy);
623 pLPAResidualsHist->Fill(restgl,ressnp);
624 pPhiYAResidualsHist->Fill(paramsTPC->Phi(),resy);
625 pPhiZAResidualsHist->Fill(paramsTPC->Phi(),resz);
626 pPtYAResidualsHist->Fill(paramsTPC->Pt(),resy);
627 pPtZAResidualsHist->Fill(paramsTPC->Pt(),resz);
628 pLowPtYAResidualsHist->Fill(paramsTPC->Pt(),resy);
629 pLowPtZAResidualsHist->Fill(paramsTPC->Pt(),resz);
633 pZYCResidualsHist->Fill(resz,resy);
634 pLPCResidualsHist->Fill(restgl,ressnp);
635 pPhiYCResidualsHist->Fill(paramsTPC->Phi(),resy);
636 pPhiZCResidualsHist->Fill(paramsTPC->Phi(),resz);
637 pPtYCResidualsHist->Fill(paramsTPC->Pt(),resy);
638 pPtZCResidualsHist->Fill(paramsTPC->Pt(),resz);
639 pLowPtYCResidualsHist->Fill(paramsTPC->Pt(),resy);
640 pLowPtZCResidualsHist->Fill(paramsTPC->Pt(),resz);