AliPhysics  ced2227 (ced2227)
AliAnalysisTaskITSTPCalignment.cxx
Go to the documentation of this file.
1 // AliAnalysisTaskITSTPCalignment
3 // Runs the relative ITS TPC alignment procedure and TPC vdrift calib
4 // Origin: Mikolaj Krzewicki, mikolaj.krzewicki@cern.ch
6 
7 #include <TChain.h>
8 #include <TTree.h>
9 #include <TH2.h>
10 #include <TProfile.h>
11 #include <TCanvas.h>
12 #include <TArrayI.h>
13 #include <TObjArray.h>
14 #include <TGraphErrors.h>
15 #include <AliMagF.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>
22 #include <AliStack.h>
23 #include <AliExternalTrackParam.h>
24 #include <AliESDtrack.h>
25 #include <AliESDfriendTrack.h>
26 #include <AliESDInputHandler.h>
27 #include "AliRelAlignerKalman.h"
30 #include "AliLog.h"
31 
33 
34 //________________________________________________________________________
37  fArrayITSglobal(0),
38  fArrayITSsa(0),
39  fDebugTree(0),
40  fAligner(0),
41  fList(0),
42  fFillDebugTree(kFALSE),
43  fDoQA(kTRUE),
44  fT0(0),
45  fTend(0),
46  fSlotWidth(0),
47  fMinPt(0.4),
48  fMinPointsVol1(3),
49  fMinPointsVol2(80),
50  fRejectOutliers(kTRUE),
51  fOutRejSigma(1.),
52  fRejectOutliersSigma2Median(kFALSE),
53  fOutRejSigma2Median(5.),
54  fOutRejSigmaOnMerge(10.),
55  fUseITSoutGlobalTrack(kFALSE),
56  fUseITSoutITSSAtrack(kTRUE)
57 {
58  //dummy ctor
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());
65 }
66 
67 //________________________________________________________________________
69  AliAnalysisTaskSE(name),
70  fArrayITSglobal(0),
71  fArrayITSsa(0),
72  fDebugTree(0),
73  fAligner(0),
74  fList(0),
75  fFillDebugTree(kFALSE),
76  fDoQA(kTRUE),
77  fT0(0),
78  fTend(0),
79  fSlotWidth(0),
80  fMinPt(0.4),
81  fMinPointsVol1(3),
82  fMinPointsVol2(80),
83  fRejectOutliers(kTRUE),
84  fOutRejSigma(1.),
85  fRejectOutliersSigma2Median(kFALSE),
86  fOutRejSigma2Median(5.),
87  fOutRejSigmaOnMerge(10.),
88  fUseITSoutGlobalTrack(kFALSE),
89  fUseITSoutITSSAtrack(kTRUE)
90 {
91  // Constructor
92 
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());
99 }
100 
101 //________________________________________________________________________
103 {
104  //ON INPUT FILE/TREE CHANGE
105 
106  //fArray->Print();
107  if (fFillDebugTree)
108  {
109  if (fAligner->GetNUpdates()>0) fDebugTree->Fill();
110  }
111 
112  return kTRUE;
113 }
114 
115 //________________________________________________________________________
117 {
118  // Create output objects
119 
121  fArrayITSglobal->SetName("array ITS global");
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);
131  fArrayITSsa->SetName("array ITS SA");
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);
140 
141  fList = new TList();
142  fList->SetName("QA");
143  fList->SetOwner(kTRUE);
144 
145  TH2F* pZYAResidualsHistBpos = new TH2F("fZYAResidualsHistBpos","z-r\\phi residuals side A (z>0), Bpos", 100, -4., 4., 100, -1.5, 1.5 );
146  pZYAResidualsHistBpos->SetXTitle("\\deltaz [cm]");
147  pZYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
148  TH2F* pZYCResidualsHistBpos = new TH2F("fZYCResidualsHistBpos","z-r\\phi residuals side C (z<0), Bpos", 100, -4., 4., 100, -1.5, 1.5 );
149  pZYCResidualsHistBpos->SetXTitle("\\deltaz [cm]");
150  pZYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
151  TH2F* pLPAResidualsHistBpos = new TH2F("fLPAResidualsHistBpos", "sin(\\phi) tan(\\lambda) residuals side A (z>0), Bpos", 100, -.05, 0.05, 100, -0.05, 0.05 );
152  pLPAResidualsHistBpos->SetXTitle("\\deltasin(\\phi)");
153  pLPAResidualsHistBpos->SetYTitle("\\deltatan(\\lambda)");
154  TH2F* pLPCResidualsHistBpos = new TH2F("fLPCResidualsHistBpos", "sin(\\phi) tan(\\lambda) residuals side C (z<0), Bpos", 100, -.05, 0.05, 100, -0.05, 0.05 );
155  pLPCResidualsHistBpos->SetXTitle("\\deltasin(\\phi)");
156  pLPCResidualsHistBpos->SetYTitle("\\deltatan(\\lambda)");
157  TH2F* pPhiYAResidualsHistBpos = new TH2F("fPhiYAResidualsHistBpos","\\phi-\\deltar\\phi side A (z>0), Bpos",36,0.0,6.3,100,-1.5,1.5);
158  pPhiYAResidualsHistBpos->SetXTitle("\\phi [rad]");
159  pPhiYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
160  TH2F* pPhiYCResidualsHistBpos = new TH2F("fPhiYCResidualsHistBpos","\\phi-\\deltar\\phi side C (z<0), Bpos",36,0.0,6.3,100,-1.5,1.5);
161  pPhiYCResidualsHistBpos->SetXTitle("\\phi [rad]");
162  pPhiYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
163  TH2F* pPhiZAResidualsHistBpos = new TH2F("fPhiZAResidualsHistBpos","\\phi-\\deltaz side A (z>0), Bpos",36,0.0,6.3,100,-3.0,3.0);
164  pPhiZAResidualsHistBpos->SetXTitle("\\phi [rad]");
165  pPhiZAResidualsHistBpos->SetYTitle("\\deltaz [cm]");
166  TH2F* pPhiZCResidualsHistBpos = new TH2F("fPhiZCResidualsHistBpos","\\phi-\\deltaz side C (z<0), Bpos",36,0.0,6.3,100,-3.0,3.0);
167  pPhiZCResidualsHistBpos->SetXTitle("\\phi [rad]");
168  pPhiZCResidualsHistBpos->SetYTitle("\\deltaz [cm]");
169  TH2F* pPtYAResidualsHistBpos = new TH2F("fPtYAResidualsHistBpos","Pt-\\deltar\\phi side A (z>0), Bpos",20,.3,8.,80,-1.5,1.5);
170  pPtYAResidualsHistBpos->SetXTitle("Pt [rad]");
171  pPtYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
172  TH2F* pPtYCResidualsHistBpos = new TH2F("fPtYCResidualsHistBpos","Pt-\\deltar\\phi side C (z<0), Bpos",20,.3,8.,80,-1.5,1.5);
173  pPtYCResidualsHistBpos->SetXTitle("Pt [rad]");
174  pPtYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
175  TH2F* pPtZAResidualsHistBpos = new TH2F("fPtZAResidualsHistBpos","Pt-\\deltaz side A (z>0), Bpos",20,.3,8.,80,-3.0,3.0);
176  pPtZAResidualsHistBpos->SetXTitle("Pt");
177  pPtZAResidualsHistBpos->SetYTitle("\\deltaz [cm]");
178  TH2F* pPtZCResidualsHistBpos = new TH2F("fPtZCResidualsHistBpos","Pt-\\deltaz side C (z<0), Bpos",20,.3,8.,80,-3.0,3.0);
179  pPtZCResidualsHistBpos->SetXTitle("Pt");
180  pPtZCResidualsHistBpos->SetYTitle("\\deltaz [cm]");
181  TH2F* pLowPtYAResidualsHistBpos = new TH2F("fLowPtYAResidualsHistBpos","Pt-\\deltar\\phi side A (z>0), Bpos",100,0.0,1.5,80,-1.5,1.5);
182  pLowPtYAResidualsHistBpos->SetXTitle("Pt [rad]");
183  pLowPtYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
184  TH2F* pLowPtYCResidualsHistBpos = new TH2F("fLowPtYCResidualsHistBpos","Pt-\\deltar\\phi side C (z<0), Bpos",100,0.0,1.5,80,-1.5,1.5);
185  pLowPtYCResidualsHistBpos->SetXTitle("Pt [rad]");
186  pLowPtYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
187  TH2F* pLowPtZAResidualsHistBpos = new TH2F("fLowPtZAResidualsHistBpos","Pt-\\deltaz side A (z>0), Bpos",100,0.0,1.5,80,-3.0,3.0);
188  pLowPtZAResidualsHistBpos->SetXTitle("Pt");
189  pLowPtZAResidualsHistBpos->SetYTitle("\\deltaz [cm]");
190  TH2F* pLowPtZCResidualsHistBpos = new TH2F("fLowPtZCResidualsHistBpos","Pt-\\deltaz side C (z<0), Bpos",100,0.0,1.5,80,-3.0,3.0);
191  pLowPtZCResidualsHistBpos->SetXTitle("Pt");
192  pLowPtZCResidualsHistBpos->SetYTitle("\\deltaz [cm]");
193  TList* listBpos = new TList();
194  listBpos->SetName("B+");
195  listBpos->SetOwner(kTRUE);
196  fList->Add(listBpos); //0
197  listBpos->Add(pZYAResidualsHistBpos);
198  listBpos->Add(pZYCResidualsHistBpos);
199  listBpos->Add(pLPAResidualsHistBpos);
200  listBpos->Add(pLPCResidualsHistBpos);
201  listBpos->Add(pPhiYAResidualsHistBpos);
202  listBpos->Add(pPhiYCResidualsHistBpos);
203  listBpos->Add(pPhiZAResidualsHistBpos);
204  listBpos->Add(pPhiZCResidualsHistBpos);
205  listBpos->Add(pPtYAResidualsHistBpos);
206  listBpos->Add(pPtYCResidualsHistBpos);
207  listBpos->Add(pPtZAResidualsHistBpos);
208  listBpos->Add(pPtZCResidualsHistBpos);
209  listBpos->Add(pLowPtYAResidualsHistBpos);
210  listBpos->Add(pLowPtYCResidualsHistBpos);
211  listBpos->Add(pLowPtZAResidualsHistBpos);
212  listBpos->Add(pLowPtZCResidualsHistBpos);
213 
214  TH2F* pZYAResidualsHistBneg = new TH2F("fZYAResidualsHistBneg","z-r\\phi residuals side A (z>0), Bneg", 100, -4., 4., 100, -1.5, 1.5 );
215  pZYAResidualsHistBneg->SetXTitle("\\deltaz [cm]");
216  pZYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
217  TH2F* pZYCResidualsHistBneg = new TH2F("fZYCResidualsHistBneg","z-r\\phi residuals side C (z<0), Bneg", 100, -4., 4., 100, -1.5, 1.5 );
218  pZYCResidualsHistBneg->SetXTitle("\\deltaz [cm]");
219  pZYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
220  TH2F* pLPAResidualsHistBneg = new TH2F("fLPAResidualsHistBneg", "sin(\\phi) tan(\\lambda) residuals side A (z>0), Bneg", 100, -.05, 0.05, 100, -0.05, 0.05 );
221  pLPAResidualsHistBneg->SetXTitle("\\deltasin(\\phi)");
222  pLPAResidualsHistBneg->SetYTitle("\\deltatan(\\lambda)");
223  TH2F* pLPCResidualsHistBneg = new TH2F("fLPCResidualsHistBneg", "sin(\\phi) tan(\\lambda) residuals side C (z<0), Bneg", 100, -.05, 0.05, 100, -0.05, 0.05 );
224  pLPCResidualsHistBneg->SetXTitle("\\deltasin(\\phi)");
225  pLPCResidualsHistBneg->SetYTitle("\\deltatan(\\lambda)");
226  TH2F* pPhiYAResidualsHistBneg = new TH2F("fPhiYAResidualsHistBneg","\\phi-\\deltar\\phi side A (z>0), Bneg",36,0.0,6.3,100,-1.5,1.5);
227  pPhiYAResidualsHistBneg->SetXTitle("\\phi [rad]");
228  pPhiYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
229  TH2F* pPhiYCResidualsHistBneg = new TH2F("fPhiYCResidualsHistBneg","\\phi-\\deltar\\phi side C (z<0), Bneg",36,0.0,6.3,100,-1.5,1.5);
230  pPhiYCResidualsHistBneg->SetXTitle("\\phi [rad]");
231  pPhiYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
232  TH2F* pPhiZAResidualsHistBneg = new TH2F("fPhiZAResidualsHistBneg","\\phi-\\deltaz side A (z>0), Bneg",36,0.0,6.3,100,-3.0,3.0);
233  pPhiZAResidualsHistBneg->SetXTitle("\\phi [rad]");
234  pPhiZAResidualsHistBneg->SetYTitle("\\deltaz [cm]");
235  TH2F* pPhiZCResidualsHistBneg = new TH2F("fPhiZCResidualsHistBneg","\\Pt-\\deltaz side C (z<0), Bneg",36,0.0,6.3,100,-3.0,3.0);
236  pPhiZCResidualsHistBneg->SetXTitle("\\phi [rad]");
237  pPhiZCResidualsHistBneg->SetYTitle("\\deltaz [cm]");
238  TH2F* pPtYAResidualsHistBneg = new TH2F("fPtYAResidualsHistBneg","Pt-\\deltar\\phi side A (z>0), Bneg",20,.3,8.,80,-1.5,1.5);
239  pPtYAResidualsHistBneg->SetXTitle("Pt [rad]");
240  pPtYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
241  TH2F* pPtYCResidualsHistBneg = new TH2F("fPtYCResidualsHistBneg","Pt-\\deltar\\phi side C (z<0), Bneg",20,.3,8.,80,-1.5,1.5);
242  pPtYCResidualsHistBneg->SetXTitle("Pt [rad]");
243  pPtYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
244  TH2F* pPtZAResidualsHistBneg = new TH2F("fPtZAResidualsHistBneg","Pt-\\deltaz side A (z>0), Bneg",20,.3,8.,80,-3.0,3.0);
245  pPtZAResidualsHistBneg->SetXTitle("Pt");
246  pPtZAResidualsHistBneg->SetYTitle("\\deltaz [cm]");
247  TH2F* pPtZCResidualsHistBneg = new TH2F("fPtZCResidualsHistBneg","Pt-\\deltaz side C (z<0), Bneg",20,.3,8.,80,-3.0,3.0);
248  pPtZCResidualsHistBneg->SetXTitle("Pt");
249  pPtZCResidualsHistBneg->SetYTitle("\\deltaz [cm]");
250  TH2F* pLowPtYAResidualsHistBneg = new TH2F("fLowPtYAResidualsHistBneg","Pt-\\deltar\\phi side A (z>0), Bneg",100,0.0,1.5,80,-1.5,1.5);
251  pLowPtYAResidualsHistBneg->SetXTitle("Pt [rad]");
252  pLowPtYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
253  TH2F* pLowPtYCResidualsHistBneg = new TH2F("fLowPtYCResidualsHistBneg","Pt-\\deltar\\phi side C (z<0), Bneg",100,0.0,1.5,80,-1.5,1.5);
254  pLowPtYCResidualsHistBneg->SetXTitle("Pt [rad]");
255  pLowPtYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
256  TH2F* pLowPtZAResidualsHistBneg = new TH2F("fLowPtZAResidualsHistBneg","Pt-\\deltaz side A (z>0), Bneg",100,0.0,1.5,80,-3.0,3.0);
257  pLowPtZAResidualsHistBneg->SetXTitle("Pt");
258  pLowPtZAResidualsHistBneg->SetYTitle("\\deltaz [cm]");
259  TH2F* pLowPtZCResidualsHistBneg = new TH2F("fLowPtZCResidualsHistBneg","Pt-\\deltaz side C (z<0), Bneg",100,0.0,1.5,80,-3.0,3.0);
260  pLowPtZCResidualsHistBneg->SetXTitle("Pt");
261  pLowPtZCResidualsHistBneg->SetYTitle("\\deltaz [cm]");
262  TList* listBneg = new TList();
263  listBneg->SetName("B-");
264  listBneg->SetOwner(kTRUE);
265  fList->Add(listBneg); //1
266  listBneg->Add(pZYAResidualsHistBneg);
267  listBneg->Add(pZYCResidualsHistBneg);
268  listBneg->Add(pLPAResidualsHistBneg);
269  listBneg->Add(pLPCResidualsHistBneg);
270  listBneg->Add(pPhiYAResidualsHistBneg);
271  listBneg->Add(pPhiYCResidualsHistBneg);
272  listBneg->Add(pPhiZAResidualsHistBneg);
273  listBneg->Add(pPhiZCResidualsHistBneg);
274  listBneg->Add(pPtYAResidualsHistBneg);
275  listBneg->Add(pPtYCResidualsHistBneg);
276  listBneg->Add(pPtZAResidualsHistBneg);
277  listBneg->Add(pPtZCResidualsHistBneg);
278  listBneg->Add(pLowPtYAResidualsHistBneg);
279  listBneg->Add(pLowPtYCResidualsHistBneg);
280  listBneg->Add(pLowPtZAResidualsHistBneg);
281  listBneg->Add(pLowPtZCResidualsHistBneg);
282 
283  TH2F* pZYAResidualsHistBnil = new TH2F("fZYAResidualsHistBnil","z-r\\phi residuals side A (z>0), Bnil", 100, -3., 3., 100, -1.5, 1.5 );
284  pZYAResidualsHistBnil->SetXTitle("\\deltaz [cm]");
285  pZYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
286  TH2F* pZYCResidualsHistBnil = new TH2F("fZYCResidualsHistBnil","z-r\\phi residuals side C (z<0), Bnil", 100, -3., 3., 100, -1.5, 1.5 );
287  pZYCResidualsHistBnil->SetXTitle("\\deltaz [cm]");
288  pZYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
289  TH2F* pLPAResidualsHistBnil = new TH2F("fLPAResidualsHistBnil", "sin(\\phi) tan(\\lambda) residuals side A (z>0), Bnil", 100, -.05, 0.05, 100, -0.05, 0.05 );
290  pLPAResidualsHistBnil->SetXTitle("\\deltasin(\\phi)");
291  pLPAResidualsHistBnil->SetYTitle("\\deltatan(\\lambda)");
292  TH2F* pLPCResidualsHistBnil = new TH2F("fLPCResidualsHistBnil", "sin(\\phi) tan(\\lambda) residuals side C (z<0), Bnil", 100, -.05, 0.05, 100, -0.05, 0.05 );
293  pLPCResidualsHistBnil->SetXTitle("\\deltasin(\\phi)");
294  pLPCResidualsHistBnil->SetYTitle("\\deltatan(\\lambda)");
295  TH2F* pPhiYAResidualsHistBnil = new TH2F("fPhiYAResidualsHistBnil","\\phi-\\deltar\\phi side A (z>0), Bnil",36,0.0,6.3,100,-1.5,1.5);
296  pPhiYAResidualsHistBnil->SetXTitle("\\phi [rad]");
297  pPhiYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
298  TH2F* pPhiYCResidualsHistBnil = new TH2F("fPhiYCResidualsHistBnil","\\phi-\\deltar\\phi side C (z<0), Bnil",36,0.0,6.3,100,-1.5,1.5);
299  pPhiYCResidualsHistBnil->SetXTitle("\\phi [rad]");
300  pPhiYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
301  TH2F* pPhiZAResidualsHistBnil = new TH2F("fPhiZAResidualsHistBnil","\\phi-\\deltaz side A (z>0), Bnil",36,0.0,6.3,100,-3.0,3.0);
302  pPhiZAResidualsHistBnil->SetXTitle("\\phi [rad]");
303  pPhiZAResidualsHistBnil->SetYTitle("\\deltaz [cm]");
304  TH2F* pPhiZCResidualsHistBnil = new TH2F("fPhiZCResidualsHistBnil","\\Pt-\\deltaz side C (z<0), Bnil",36,0.0,6.3,100,-3.0,3.0);
305  pPhiZCResidualsHistBnil->SetXTitle("\\phi [rad]");
306  pPhiZCResidualsHistBnil->SetYTitle("\\deltaz [cm]");
307  TH2F* pPtYAResidualsHistBnil = new TH2F("fPtYAResidualsHistBnil","Pt-\\deltar\\phi side A (z>0), Bnil",20,.3,8.,80,-1.5,1.5);
308  pPtYAResidualsHistBnil->SetXTitle("Pt [rad]");
309  pPtYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
310  TH2F* pPtYCResidualsHistBnil = new TH2F("fPtYCResidualsHistBnil","Pt-\\deltar\\phi side C (z<0), Bnil",20,.3,8.,80,-1.5,1.5);
311  pPtYCResidualsHistBnil->SetXTitle("Pt [rad]");
312  pPtYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
313  TH2F* pPtZAResidualsHistBnil = new TH2F("fPtZAResidualsHistBnil","Pt-\\deltaz side A (z>0), Bnil",20,.3,8.,80,-3.0,3.0);
314  pPtZAResidualsHistBnil->SetXTitle("Pt");
315  pPtZAResidualsHistBnil->SetYTitle("\\deltaz [cm]");
316  TH2F* pPtZCResidualsHistBnil = new TH2F("fPtZCResidualsHistBnil","Pt-\\deltaz side C (z<0), Bnil",20,.3,8.,80,-3.0,3.0);
317  pPtZCResidualsHistBnil->SetXTitle("Pt");
318  pPtZCResidualsHistBnil->SetYTitle("\\deltaz [cm]");
319  TH2F* pLowPtYAResidualsHistBnil = new TH2F("fLowPtYAResidualsHistBnil","Pt-\\deltar\\phi side A (z>0), Bnil",100,0.0,1.5,80,-1.5,1.5);
320  pLowPtYAResidualsHistBnil->SetXTitle("Pt [rad]");
321  pLowPtYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
322  TH2F* pLowPtYCResidualsHistBnil = new TH2F("fLowPtYCResidualsHistBnil","Pt-\\deltar\\phi side C (z<0), Bnil",100,0.0,1.5,80,-1.5,1.5);
323  pLowPtYCResidualsHistBnil->SetXTitle("Pt [rad]");
324  pLowPtYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
325  TH2F* pLowPtZAResidualsHistBnil = new TH2F("fLowPtZAResidualsHistBnil","Pt-\\deltaz side A (z>0), Bnil",100,0.0,1.5,80,-3.0,3.0);
326  pLowPtZAResidualsHistBnil->SetXTitle("Pt");
327  pLowPtZAResidualsHistBnil->SetYTitle("\\deltaz [cm]");
328  TH2F* pLowPtZCResidualsHistBnil = new TH2F("fLowPtZCResidualsHistBnil","Pt-\\deltaz side C (z<0), Bnil",100,0.0,1.5,80,-3.0,3.0);
329  pLowPtZCResidualsHistBnil->SetXTitle("Pt");
330  pLowPtZCResidualsHistBnil->SetYTitle("\\deltaz [cm]");
331  TList* listBnil = new TList();
332  listBnil->SetName("B0");
333  listBnil->SetOwner(kTRUE);
334  fList->Add(listBnil); //2
335  listBnil->Add(pZYAResidualsHistBnil);
336  listBnil->Add(pZYCResidualsHistBnil);
337  listBnil->Add(pLPAResidualsHistBnil);
338  listBnil->Add(pLPCResidualsHistBnil);
339  listBnil->Add(pPhiYAResidualsHistBnil);
340  listBnil->Add(pPhiYCResidualsHistBnil);
341  listBnil->Add(pPhiZAResidualsHistBnil);
342  listBnil->Add(pPhiZCResidualsHistBnil);
343  listBnil->Add(pPtYAResidualsHistBnil);
344  listBnil->Add(pPtYCResidualsHistBnil);
345  listBnil->Add(pPtZAResidualsHistBnil);
346  listBnil->Add(pPtZCResidualsHistBnil);
347  listBnil->Add(pLowPtYAResidualsHistBnil);
348  listBnil->Add(pLowPtYCResidualsHistBnil);
349  listBnil->Add(pLowPtZAResidualsHistBnil);
350  listBnil->Add(pLowPtZCResidualsHistBnil);
351 
352  TH1F* pNmatchingEff=new TH1F("pNmatchingEff","matching efficiency",50,0.,1.);
353  fList->Add(pNmatchingEff); //3
354 
355  TH1I* pChecks = new TH1I("some counters","some counters",10,0,10);
356  pChecks->GetXaxis()->SetBinLabel(1+kNoESD,"no ESD");
357  pChecks->GetXaxis()->SetBinLabel(1+kNoESDfriend,"no ESDfriend");
358  pChecks->GetXaxis()->SetBinLabel(1+kNoFriendTrack,"no friendtrack");
359  pChecks->GetXaxis()->SetBinLabel(1+kNoITSoutParams,"no ITS out params");
360  pChecks->GetXaxis()->SetBinLabel(1+kESDfriend,"ESD friend");
361  pChecks->GetXaxis()->SetBinLabel(1+kFriendsSkipBit,"friends+skipbit");
362  pChecks->GetXaxis()->SetBinLabel(1+kFriendTrack,"friendtrack");
363  pChecks->GetXaxis()->SetBinLabel(1+kITSoutParams,"ITS out params");
364  fList->Add(pChecks); //4
365 
366  fAligner = new AliRelAlignerKalman();
367  fAligner->SetRejectOutliers(fRejectOutliers);
368  fAligner->SetOutRejSigma(fOutRejSigma);
369  fAligner->SetRejectOutliersSigma2Median(fRejectOutliersSigma2Median);
370  fAligner->SetOutRejSigma2Median(fOutRejSigma2Median);
371  fList->Add(fAligner);
372 
373  fDebugTree = new TTree("debugTree","tree with debug info");
374  fDebugTree->Branch("aligner","AliRelAlignerKalman",&fAligner);
375 
376  // Post output data.
377  PostData(0, fDebugTree);
378  PostData(1, fList);
379  PostData(2, fArrayITSglobal);
380  PostData(3, fArrayITSsa);
381 }
382 
383 //________________________________________________________________________
385 {
386  // Main loop
387  // Called for each event
388 
389  TH1I* pChecks = static_cast<TH1I*>(fList->At(4));
390  //check for ESD and friends
391  AliESDEvent* event = dynamic_cast<AliESDEvent*>(InputEvent());
392  AliESDfriend* eventFriend = dynamic_cast<AliESDfriend*>(ESDfriend());
393  if (!event)
394  {
395  AliError("no ESD");
396  pChecks->Fill(kNoESD);
397  return;
398  }
399 
400  if (!eventFriend)
401  {
402  AliError("no ESD friend");
403  pChecks->Fill(kNoESDfriend);
404  return;
405  }
406  else
407  {
408  pChecks->Fill(kESDfriend);
409  }
410 
411  //event->SetESDfriend(eventFriend);
412 
413  if (eventFriend->TestSkipBit())
414  {
415  pChecks->Fill(kFriendsSkipBit);
416  return;
417  }
418 
419  //Update the parmeters
420  AnalyzeESDevent(event);
421 }
422 
423 //________________________________________________________________________
425 {
426  //analyze an ESD event with track matching
427  Int_t ntracks = event->GetNumberOfTracks();
428  if (ntracks==0) return;
429 
431  {
432  AliRelAlignerKalman* alignerGlobal = fArrayITSglobal->GetAligner(event);
433  if (alignerGlobal)
434  alignerGlobal->AddESDevent(event);
435  }
436 
438  {
439  //get the aligner for the event
440  //do it with matching
441  TObjArray arrayParamITS(ntracks);
442  TObjArray arrayParamTPC(ntracks);
443 
444  Int_t n = FindMatchingTracks(arrayParamITS, arrayParamTPC, event);
445  AliInfo(Form("matched tracklet pairs: %i\n",n));
446 
447  TH1F* nMatchingEff=static_cast<TH1F*>(fList->At(3));
448  Float_t ratio = (float)n/(float)ntracks;
449  nMatchingEff->Fill(ratio);
450 
451  AliRelAlignerKalman* alignerSA = NULL;
452  if (n>0) alignerSA = fArrayITSsa->GetAligner(event);
453 
454  for (Int_t i=0;i<n;i++)
455  {
456  AliExternalTrackParam* paramsITS=static_cast<AliExternalTrackParam*>(arrayParamITS[i]);
457  AliExternalTrackParam* paramsTPC=static_cast<AliExternalTrackParam*>(arrayParamTPC[i]);
458 
459  if (!(paramsITS&&paramsTPC)) continue;
460 
461  //QA
462  if (fDoQA) DoQA(paramsITS,paramsTPC);
463 
464  //debug
465  if (fAligner->AddTrackParams(paramsITS, paramsTPC))
466  {
467  fAligner->SetRunNumber(event->GetRunNumber());
468  fAligner->SetMagField(event->GetMagneticField());
469  fAligner->SetTimeStamp(event->GetTimeStamp());
470  }
471 
472  //alignment
473  if (alignerSA) alignerSA->AddTrackParams(paramsITS, paramsTPC);
474  }
475  arrayParamITS.Delete();
476  arrayParamTPC.Delete();
477  }
478 }
479 
480 //________________________________________________________________________
482 {
483  // Called once at the end of the query
484  fArrayITSsa->Print("a");
485  fAligner->Print();
486 }
487 
488 //________________________________________________________________________
490 {
491  //find matching tracks and return tobjarrays with the params
492  //fUniqueID of param object contains tracknumber in event.
493 
494  TH1I* pChecks = static_cast<TH1I*>(fList->At(4));
495  Int_t ntracks = pESD->GetNumberOfTracks();
496  Double_t magfield = pESD->GetMagneticField();
497 
498  Int_t* matchedArr = new Int_t[ntracks]; //storage for index of ITS track for which a match was found
499  for (Int_t i=0;i<ntracks;i++)
500  {
501  matchedArr[i]=-1;
502  }
503 
504  Int_t iMatched=-1;
505  for (Int_t i=0; i<ntracks; i++)
506  {
507  //get track1 and friend
508  AliESDtrack* track1 = pESD->GetTrack(i);
509  if (!track1) continue;
510 
511  if (track1->GetNcls(0) < fMinPointsVol1) continue;
512 
513  if (!( ( track1->IsOn(AliESDtrack::kITSrefit)) &&
514  (!track1->IsOn(AliESDtrack::kTPCin)) )) continue;
515 
516  const AliESDfriendTrack* constfriendtrack1 = track1->GetFriendTrack();
517  if (!constfriendtrack1) {AliInfo(Form("no friendTrack1\n"));pChecks->Fill(kNoFriendTrack);continue;}
518  pChecks->Fill(kFriendTrack);
519  AliESDfriendTrack friendtrack1(*constfriendtrack1);
520 
521  if (!friendtrack1.GetITSOut()) {pChecks->Fill(kNoITSoutParams);continue;}
522  pChecks->Fill(kITSoutParams);
523  AliExternalTrackParam params1(*(friendtrack1.GetITSOut()));
524 
525  Double_t bestd = 1000.; //best distance
526  Bool_t newi = kTRUE; //whether we start with a new i
527  for (Int_t j=0; j<ntracks; j++)
528  {
529  if (matchedArr[j]>0 && matchedArr[j]!=i) continue; //already matched, everything tried
530  //get track2 and friend
531  AliESDtrack* track2 = pESD->GetTrack(j);
532  if (!track2) continue;
533  if (track1==track2) continue;
534  if (!(track2->IsOn(AliESDtrack::kTPCin))) continue; //must be TPC
535  if (!(track2->IsOn(AliESDtrack::kITSout))) continue; //must have ITS
536 
537  //if (track2->GetNcls(0) != track1->GetNcls(0)) continue;
538  //if (track2->GetITSClusterMap() != track1->GetITSClusterMap()) continue;
539  if (track2->GetNcls(1) < fMinPointsVol2) continue; //min 80 clusters in TPC
540  if (track2->GetTgl() > 1.) continue; //acceptance
541  //cut crossing tracks
542  if (!track2->GetInnerParam()) continue;
543  if (!track2->GetOuterParam()) continue;
544  if (track2->GetOuterParam()->GetZ()*track2->GetInnerParam()->GetZ()<0) continue;
545  if (track2->GetInnerParam()->GetX()>90) continue;
546  if (TMath::Abs(track2->GetInnerParam()->GetZ())<5.) continue; //too close to membrane?
547 
548  AliExternalTrackParam params2(*(track2->GetInnerParam()));
549 
550  //bring to same reference plane
551  if (!params2.Rotate(params1.GetAlpha())) continue;
552  if (!params2.PropagateTo(params1.GetX(), magfield)) continue;
553 
554  //pt cut, only for data with magfield
555  if ((params2.Pt()<fMinPt)&&(TMath::Abs(magfield)>1.)) continue;
556  //else
557  //if (fMC)
558  // {
559  // AliStack* stack = fMC->Stack();
560  // Int_t label = track2->GetLabel();
561  // if (label<0) continue;
562  // TParticle* particle = stack->Particle(label);
563  // if (particle->Pt()<fMinPt) continue;
564  // }
565 
566  const Double32_t* p1 = params1.GetParameter();
567  const Double32_t* p2 = params2.GetParameter();
568 
569  //hard cuts
570  Double_t dy = TMath::Abs(p2[0]-p1[0]);
571  Double_t dz = TMath::Abs(p2[1]-p1[1]);
572  Double_t dphi = TMath::Abs(p2[2]-p1[2]);
573  Double_t dlam = TMath::Abs(p2[3]-p1[3]);
574  if (dy > 2.0) continue;
575  if (dz > 10.0) continue;
576  if (dphi > 0.1 ) continue;
577  if (dlam > 0.1 ) continue;
578 
579  //best match only
580  Double_t d = TMath::Sqrt(dy*dy+dz*dz+dphi*dphi+dlam*dlam);
581  if ( d >= bestd) continue;
582  bestd = d;
583  matchedArr[j]=i; //j-th track matches i-th (ITS) track
584  if (newi) iMatched++; newi=kFALSE; //increment at most once per i
585  params1.SetUniqueID(i); //store tracknummer
586  params2.SetUniqueID(j);
587  if (arrITS[iMatched] && arrTPC[iMatched])
588  {
589  *(arrITS[iMatched]) = params1;
590  *(arrTPC[iMatched]) = params2;
591  }
592  else
593  {
594  arrITS[iMatched] = new AliExternalTrackParam(params1);
595  arrTPC[iMatched] = new AliExternalTrackParam(params2);
596  }//else
597  }//for j
598  }//for i
599  delete [] matchedArr;
600  return iMatched+1;
601 }
602 
603 //________________________________________________________________________
604 void AliAnalysisTaskITSTPCalignment::DoQA(AliExternalTrackParam* paramsITS,
605  AliExternalTrackParam* paramsTPC)
606 {
607  //fill qa histograms in a given list, per field direction (5,-5,0)
608  Float_t resy = paramsTPC->GetY() - paramsITS->GetY();
609  Float_t resz = paramsTPC->GetZ() - paramsITS->GetZ();
610  Float_t ressnp = paramsTPC->GetSnp() - paramsITS->GetSnp();
611  Float_t restgl = paramsTPC->GetTgl() - paramsITS->GetTgl();
612 
613  TList* pList=NULL;
614  Double_t field = fInputEvent->GetMagneticField();
615  if (field >= 1.) pList = static_cast<TList*>(fList->At(0));
616  if (field <= -1.) pList = static_cast<TList*>(fList->At(1));
617  if (field < 1. && field > -1.) pList = static_cast<TList*>(fList->At(2));
618  if (!pList) return;
619 
620  TH2F* pZYAResidualsHist = static_cast<TH2F*>(pList->At(0));
621  TH2F* pZYCResidualsHist = static_cast<TH2F*>(pList->At(1));
622  TH2F* pLPAResidualsHist = static_cast<TH2F*>(pList->At(2));
623  TH2F* pLPCResidualsHist = static_cast<TH2F*>(pList->At(3));
624  TH2F* pPhiYAResidualsHist = static_cast<TH2F*>(pList->At(4));
625  TH2F* pPhiYCResidualsHist = static_cast<TH2F*>(pList->At(5));
626  TH2F* pPhiZAResidualsHist = static_cast<TH2F*>(pList->At(6));
627  TH2F* pPhiZCResidualsHist = static_cast<TH2F*>(pList->At(7));
628  TH2F* pPtYAResidualsHist = static_cast<TH2F*>(pList->At(8));
629  TH2F* pPtYCResidualsHist = static_cast<TH2F*>(pList->At(9));
630  TH2F* pPtZAResidualsHist = static_cast<TH2F*>(pList->At(10));
631  TH2F* pPtZCResidualsHist = static_cast<TH2F*>(pList->At(11));
632  TH2F* pLowPtYAResidualsHist = static_cast<TH2F*>(pList->At(12));
633  TH2F* pLowPtYCResidualsHist = static_cast<TH2F*>(pList->At(13));
634  TH2F* pLowPtZAResidualsHist = static_cast<TH2F*>(pList->At(14));
635  TH2F* pLowPtZCResidualsHist = static_cast<TH2F*>(pList->At(15));
636 
637  if (paramsITS->GetZ() > 0.0)
638  {
639  pZYAResidualsHist->Fill(resz,resy);
640  pLPAResidualsHist->Fill(restgl,ressnp);
641  pPhiYAResidualsHist->Fill(paramsTPC->Phi(),resy);
642  pPhiZAResidualsHist->Fill(paramsTPC->Phi(),resz);
643  pPtYAResidualsHist->Fill(paramsTPC->Pt(),resy);
644  pPtZAResidualsHist->Fill(paramsTPC->Pt(),resz);
645  pLowPtYAResidualsHist->Fill(paramsTPC->Pt(),resy);
646  pLowPtZAResidualsHist->Fill(paramsTPC->Pt(),resz);
647  }
648  else
649  {
650  pZYCResidualsHist->Fill(resz,resy);
651  pLPCResidualsHist->Fill(restgl,ressnp);
652  pPhiYCResidualsHist->Fill(paramsTPC->Phi(),resy);
653  pPhiZCResidualsHist->Fill(paramsTPC->Phi(),resz);
654  pPtYCResidualsHist->Fill(paramsTPC->Pt(),resy);
655  pPtZCResidualsHist->Fill(paramsTPC->Pt(),resz);
656  pLowPtYCResidualsHist->Fill(paramsTPC->Pt(),resy);
657  pLowPtZCResidualsHist->Fill(paramsTPC->Pt(),resz);
658  }
659 }
double Double_t
Definition: External.C:58
Definition: External.C:236
Bool_t fFillDebugTree
list with QA histograms
void DoQA(AliExternalTrackParam *paramsITS, AliExternalTrackParam *paramsTPC)
AliRelAlignerKalman * GetAligner(UInt_t timestamp)
void SetupArray(Int_t t0, Int_t tend, Int_t slotwidth)
int Int_t
Definition: External.C:63
Definition: External.C:204
virtual void Print(Option_t *option="") const
float Float_t
Definition: External.C:68
Int_t FindMatchingTracks(TObjArray &arrITS, TObjArray &arrTPC, AliESDEvent *pESD)
AliRelAlignerKalman * GetAlignerTemplate()
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
TTree * fDebugTree
array of aligners ITS standalone
AliRelAlignerKalmanArray * fArrayITSsa
array of aligners with ITS global