AliPhysics  958ad07 (958ad07)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
28 #include "AliRelAlignerKalmanArray.h"
29 #include "AliAnalysisTaskITSTPCalignment.h"
30 #include "AliLog.h"
31 
32 ClassImp(AliAnalysisTaskITSTPCalignment)
33 
34 //________________________________________________________________________
35 AliAnalysisTaskITSTPCalignment::AliAnalysisTaskITSTPCalignment():
36  AliAnalysisTaskSE(),
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 //________________________________________________________________________
68 AliAnalysisTaskITSTPCalignment::AliAnalysisTaskITSTPCalignment(const char *name):
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 //________________________________________________________________________
102 Bool_t AliAnalysisTaskITSTPCalignment::UserNotify()
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 //________________________________________________________________________
116 void AliAnalysisTaskITSTPCalignment::UserCreateOutputObjects()
117 {
118  // Create output objects
119 
120  fArrayITSglobal = new AliRelAlignerKalmanArray();
121  fArrayITSglobal->SetName("array ITS global");
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("array ITS SA");
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);
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 //________________________________________________________________________
384 void AliAnalysisTaskITSTPCalignment::UserExec(Option_t *)
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 //________________________________________________________________________
424 void AliAnalysisTaskITSTPCalignment::AnalyzeESDevent(AliESDEvent* event)
425 {
426  //analyze an ESD event with track matching
427  Int_t ntracks = event->GetNumberOfTracks();
428  if (ntracks==0) return;
429 
430  if (fUseITSoutGlobalTrack)
431  {
432  AliRelAlignerKalman* alignerGlobal = fArrayITSglobal->GetAligner(event);
433  if (alignerGlobal)
434  alignerGlobal->AddESDevent(event);
435  }
436 
437  if (fUseITSoutITSSAtrack)
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 //________________________________________________________________________
481 void AliAnalysisTaskITSTPCalignment::Terminate(Option_t *)
482 {
483  // Called once at the end of the query
484  fArrayITSsa->Print("a");
485  fAligner->Print();
486 }
487 
488 //________________________________________________________________________
489 Int_t AliAnalysisTaskITSTPCalignment::FindMatchingTracks(TObjArray& arrITS, TObjArray& arrTPC, AliESDEvent* pESD)
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_t
double Double_t
Definition: External.C:58
TH2F
Definition: External.C:236
AliAnalysisTaskITSTPCalignment::fFillDebugTree
Bool_t fFillDebugTree
list with QA histograms
Definition: AliAnalysisTaskITSTPCalignment.h:73
AliAnalysisTaskITSTPCalignment::DoQA
void DoQA(AliExternalTrackParam *paramsITS, AliExternalTrackParam *paramsTPC)
Definition: AliAnalysisTaskITSTPCalignment.cxx:604
AliRelAlignerKalmanArray::GetAligner
AliRelAlignerKalman * GetAligner(UInt_t timestamp)
Definition: AliRelAlignerKalmanArray.cxx:180
AliRelAlignerKalmanArray.h
ClassImp
ClassImp(AliAnalysisTaskITSTPCalignment) AliAnalysisTaskITSTPCalignment
Definition: AliAnalysisTaskITSTPCalignment.cxx:32
AliRelAlignerKalmanArray::SetupArray
void SetupArray(Int_t t0, Int_t tend, Int_t slotwidth)
Definition: AliRelAlignerKalmanArray.cxx:112
Int_t
int Int_t
Definition: External.C:63
TH1I
Definition: External.C:204
AliRelAlignerKalmanArray::Print
virtual void Print(Option_t *option="") const
Definition: AliRelAlignerKalmanArray.cxx:280
Float_t
float Float_t
Definition: External.C:68
AliAnalysisTaskITSTPCalignment::FindMatchingTracks
Int_t FindMatchingTracks(TObjArray &arrITS, TObjArray &arrTPC, AliESDEvent *pESD)
Definition: AliAnalysisTaskITSTPCalignment.cxx:489
AliRelAlignerKalmanArray::GetAlignerTemplate
AliRelAlignerKalman * GetAlignerTemplate()
Definition: AliRelAlignerKalmanArray.cxx:132
Option_t
const char Option_t
Definition: External.C:48
Bool_t
bool Bool_t
Definition: External.C:53
AliAnalysisTaskITSTPCalignment::fDebugTree
TTree * fDebugTree
array of aligners ITS standalone
Definition: AliAnalysisTaskITSTPCalignment.h:69
AliAnalysisTaskITSTPCalignment.h
AliAnalysisTaskITSTPCalignment::fArrayITSsa
AliRelAlignerKalmanArray * fArrayITSsa
array of aligners with ITS global
Definition: AliAnalysisTaskITSTPCalignment.h:68