Fix for coverity (AdC)
[u/mrichter/AliRoot.git] / PWGPP / AliAnalysisTaskITSTPCalignment.cxx
1 ////////////////////////////////////////////////////////////////////////////////
2 //  AliAnalysisTaskITSTPCalignment
3 //  Runs the relative ITS TPC alignment procedure and TPC vdrift calib
4 //  Origin: Mikolaj Krzewicki, mikolaj.krzewicki@cern.ch
5 ////////////////////////////////////////////////////////////////////////////////
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 }