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