Using trigger selection as an optional argument for the
[u/mrichter/AliRoot.git] / PWGPP / 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
63f4d30f 119
51cfc50d 120 fArrayITSglobal = new AliRelAlignerKalmanArray();
63f4d30f 121 fArrayITSglobal->SetName("array ITS global");
51cfc50d 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();
63f4d30f 131 fArrayITSsa->SetName("array ITS SA");
51cfc50d 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();
63f4d30f 142 fList->SetName("QA");
143 fList->SetOwner(kTRUE);
4a84c20d 144
d6d927ac 145 TH2F* pZYAResidualsHistBpos = new TH2F("fZYAResidualsHistBpos","z-r\\phi residuals side A (z>0), Bpos", 100, -4., 4., 100, -1.5, 1.5 );
51cfc50d 146 pZYAResidualsHistBpos->SetXTitle("\\deltaz [cm]");
147 pZYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 148 TH2F* pZYCResidualsHistBpos = new TH2F("fZYCResidualsHistBpos","z-r\\phi residuals side C (z<0), Bpos", 100, -4., 4., 100, -1.5, 1.5 );
51cfc50d 149 pZYCResidualsHistBpos->SetXTitle("\\deltaz [cm]");
150 pZYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 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 );
51cfc50d 152 pLPAResidualsHistBpos->SetXTitle("\\deltasin(\\phi)");
153 pLPAResidualsHistBpos->SetYTitle("\\deltatan(\\lambda)");
d6d927ac 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 );
51cfc50d 155 pLPCResidualsHistBpos->SetXTitle("\\deltasin(\\phi)");
156 pLPCResidualsHistBpos->SetYTitle("\\deltatan(\\lambda)");
d6d927ac 157 TH2F* pPhiYAResidualsHistBpos = new TH2F("fPhiYAResidualsHistBpos","\\phi-\\deltar\\phi side A (z>0), Bpos",36,0.0,6.3,100,-1.5,1.5);
51cfc50d 158 pPhiYAResidualsHistBpos->SetXTitle("\\phi [rad]");
159 pPhiYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 160 TH2F* pPhiYCResidualsHistBpos = new TH2F("fPhiYCResidualsHistBpos","\\phi-\\deltar\\phi side C (z<0), Bpos",36,0.0,6.3,100,-1.5,1.5);
51cfc50d 161 pPhiYCResidualsHistBpos->SetXTitle("\\phi [rad]");
162 pPhiYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 163 TH2F* pPhiZAResidualsHistBpos = new TH2F("fPhiZAResidualsHistBpos","\\phi-\\deltaz side A (z>0), Bpos",36,0.0,6.3,100,-3.0,3.0);
51cfc50d 164 pPhiZAResidualsHistBpos->SetXTitle("\\phi [rad]");
165 pPhiZAResidualsHistBpos->SetYTitle("\\deltaz [cm]");
d6d927ac 166 TH2F* pPhiZCResidualsHistBpos = new TH2F("fPhiZCResidualsHistBpos","\\phi-\\deltaz side C (z<0), Bpos",36,0.0,6.3,100,-3.0,3.0);
51cfc50d 167 pPhiZCResidualsHistBpos->SetXTitle("\\phi [rad]");
168 pPhiZCResidualsHistBpos->SetYTitle("\\deltaz [cm]");
d6d927ac 169 TH2F* pPtYAResidualsHistBpos = new TH2F("fPtYAResidualsHistBpos","Pt-\\deltar\\phi side A (z>0), Bpos",20,.3,8.,80,-1.5,1.5);
51cfc50d 170 pPtYAResidualsHistBpos->SetXTitle("Pt [rad]");
171 pPtYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 172 TH2F* pPtYCResidualsHistBpos = new TH2F("fPtYCResidualsHistBpos","Pt-\\deltar\\phi side C (z<0), Bpos",20,.3,8.,80,-1.5,1.5);
51cfc50d 173 pPtYCResidualsHistBpos->SetXTitle("Pt [rad]");
174 pPtYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 175 TH2F* pPtZAResidualsHistBpos = new TH2F("fPtZAResidualsHistBpos","Pt-\\deltaz side A (z>0), Bpos",20,.3,8.,80,-3.0,3.0);
51cfc50d 176 pPtZAResidualsHistBpos->SetXTitle("Pt");
177 pPtZAResidualsHistBpos->SetYTitle("\\deltaz [cm]");
d6d927ac 178 TH2F* pPtZCResidualsHistBpos = new TH2F("fPtZCResidualsHistBpos","Pt-\\deltaz side C (z<0), Bpos",20,.3,8.,80,-3.0,3.0);
51cfc50d 179 pPtZCResidualsHistBpos->SetXTitle("Pt");
180 pPtZCResidualsHistBpos->SetYTitle("\\deltaz [cm]");
d6d927ac 181 TH2F* pLowPtYAResidualsHistBpos = new TH2F("fLowPtYAResidualsHistBpos","Pt-\\deltar\\phi side A (z>0), Bpos",100,0.0,1.5,80,-1.5,1.5);
51cfc50d 182 pLowPtYAResidualsHistBpos->SetXTitle("Pt [rad]");
183 pLowPtYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 184 TH2F* pLowPtYCResidualsHistBpos = new TH2F("fLowPtYCResidualsHistBpos","Pt-\\deltar\\phi side C (z<0), Bpos",100,0.0,1.5,80,-1.5,1.5);
51cfc50d 185 pLowPtYCResidualsHistBpos->SetXTitle("Pt [rad]");
186 pLowPtYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 187 TH2F* pLowPtZAResidualsHistBpos = new TH2F("fLowPtZAResidualsHistBpos","Pt-\\deltaz side A (z>0), Bpos",100,0.0,1.5,80,-3.0,3.0);
51cfc50d 188 pLowPtZAResidualsHistBpos->SetXTitle("Pt");
189 pLowPtZAResidualsHistBpos->SetYTitle("\\deltaz [cm]");
d6d927ac 190 TH2F* pLowPtZCResidualsHistBpos = new TH2F("fLowPtZCResidualsHistBpos","Pt-\\deltaz side C (z<0), Bpos",100,0.0,1.5,80,-3.0,3.0);
51cfc50d 191 pLowPtZCResidualsHistBpos->SetXTitle("Pt");
192 pLowPtZCResidualsHistBpos->SetYTitle("\\deltaz [cm]");
193 TList* listBpos = new TList();
63f4d30f 194 listBpos->SetName("B+");
195 listBpos->SetOwner(kTRUE);
d6d927ac 196 fList->Add(listBpos); //0
51cfc50d 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);
4a84c20d 213
d6d927ac 214 TH2F* pZYAResidualsHistBneg = new TH2F("fZYAResidualsHistBneg","z-r\\phi residuals side A (z>0), Bneg", 100, -4., 4., 100, -1.5, 1.5 );
51cfc50d 215 pZYAResidualsHistBneg->SetXTitle("\\deltaz [cm]");
216 pZYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 217 TH2F* pZYCResidualsHistBneg = new TH2F("fZYCResidualsHistBneg","z-r\\phi residuals side C (z<0), Bneg", 100, -4., 4., 100, -1.5, 1.5 );
51cfc50d 218 pZYCResidualsHistBneg->SetXTitle("\\deltaz [cm]");
219 pZYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 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 );
51cfc50d 221 pLPAResidualsHistBneg->SetXTitle("\\deltasin(\\phi)");
222 pLPAResidualsHistBneg->SetYTitle("\\deltatan(\\lambda)");
d6d927ac 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 );
51cfc50d 224 pLPCResidualsHistBneg->SetXTitle("\\deltasin(\\phi)");
225 pLPCResidualsHistBneg->SetYTitle("\\deltatan(\\lambda)");
d6d927ac 226 TH2F* pPhiYAResidualsHistBneg = new TH2F("fPhiYAResidualsHistBneg","\\phi-\\deltar\\phi side A (z>0), Bneg",36,0.0,6.3,100,-1.5,1.5);
51cfc50d 227 pPhiYAResidualsHistBneg->SetXTitle("\\phi [rad]");
228 pPhiYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 229 TH2F* pPhiYCResidualsHistBneg = new TH2F("fPhiYCResidualsHistBneg","\\phi-\\deltar\\phi side C (z<0), Bneg",36,0.0,6.3,100,-1.5,1.5);
51cfc50d 230 pPhiYCResidualsHistBneg->SetXTitle("\\phi [rad]");
231 pPhiYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 232 TH2F* pPhiZAResidualsHistBneg = new TH2F("fPhiZAResidualsHistBneg","\\phi-\\deltaz side A (z>0), Bneg",36,0.0,6.3,100,-3.0,3.0);
51cfc50d 233 pPhiZAResidualsHistBneg->SetXTitle("\\phi [rad]");
234 pPhiZAResidualsHistBneg->SetYTitle("\\deltaz [cm]");
d6d927ac 235 TH2F* pPhiZCResidualsHistBneg = new TH2F("fPhiZCResidualsHistBneg","\\Pt-\\deltaz side C (z<0), Bneg",36,0.0,6.3,100,-3.0,3.0);
51cfc50d 236 pPhiZCResidualsHistBneg->SetXTitle("\\phi [rad]");
237 pPhiZCResidualsHistBneg->SetYTitle("\\deltaz [cm]");
d6d927ac 238 TH2F* pPtYAResidualsHistBneg = new TH2F("fPtYAResidualsHistBneg","Pt-\\deltar\\phi side A (z>0), Bneg",20,.3,8.,80,-1.5,1.5);
51cfc50d 239 pPtYAResidualsHistBneg->SetXTitle("Pt [rad]");
240 pPtYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 241 TH2F* pPtYCResidualsHistBneg = new TH2F("fPtYCResidualsHistBneg","Pt-\\deltar\\phi side C (z<0), Bneg",20,.3,8.,80,-1.5,1.5);
51cfc50d 242 pPtYCResidualsHistBneg->SetXTitle("Pt [rad]");
243 pPtYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 244 TH2F* pPtZAResidualsHistBneg = new TH2F("fPtZAResidualsHistBneg","Pt-\\deltaz side A (z>0), Bneg",20,.3,8.,80,-3.0,3.0);
51cfc50d 245 pPtZAResidualsHistBneg->SetXTitle("Pt");
246 pPtZAResidualsHistBneg->SetYTitle("\\deltaz [cm]");
d6d927ac 247 TH2F* pPtZCResidualsHistBneg = new TH2F("fPtZCResidualsHistBneg","Pt-\\deltaz side C (z<0), Bneg",20,.3,8.,80,-3.0,3.0);
51cfc50d 248 pPtZCResidualsHistBneg->SetXTitle("Pt");
249 pPtZCResidualsHistBneg->SetYTitle("\\deltaz [cm]");
d6d927ac 250 TH2F* pLowPtYAResidualsHistBneg = new TH2F("fLowPtYAResidualsHistBneg","Pt-\\deltar\\phi side A (z>0), Bneg",100,0.0,1.5,80,-1.5,1.5);
51cfc50d 251 pLowPtYAResidualsHistBneg->SetXTitle("Pt [rad]");
252 pLowPtYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 253 TH2F* pLowPtYCResidualsHistBneg = new TH2F("fLowPtYCResidualsHistBneg","Pt-\\deltar\\phi side C (z<0), Bneg",100,0.0,1.5,80,-1.5,1.5);
51cfc50d 254 pLowPtYCResidualsHistBneg->SetXTitle("Pt [rad]");
255 pLowPtYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 256 TH2F* pLowPtZAResidualsHistBneg = new TH2F("fLowPtZAResidualsHistBneg","Pt-\\deltaz side A (z>0), Bneg",100,0.0,1.5,80,-3.0,3.0);
51cfc50d 257 pLowPtZAResidualsHistBneg->SetXTitle("Pt");
258 pLowPtZAResidualsHistBneg->SetYTitle("\\deltaz [cm]");
d6d927ac 259 TH2F* pLowPtZCResidualsHistBneg = new TH2F("fLowPtZCResidualsHistBneg","Pt-\\deltaz side C (z<0), Bneg",100,0.0,1.5,80,-3.0,3.0);
51cfc50d 260 pLowPtZCResidualsHistBneg->SetXTitle("Pt");
261 pLowPtZCResidualsHistBneg->SetYTitle("\\deltaz [cm]");
262 TList* listBneg = new TList();
63f4d30f 263 listBneg->SetName("B-");
264 listBneg->SetOwner(kTRUE);
d6d927ac 265 fList->Add(listBneg); //1
51cfc50d 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
d6d927ac 283 TH2F* pZYAResidualsHistBnil = new TH2F("fZYAResidualsHistBnil","z-r\\phi residuals side A (z>0), Bnil", 100, -3., 3., 100, -1.5, 1.5 );
51cfc50d 284 pZYAResidualsHistBnil->SetXTitle("\\deltaz [cm]");
285 pZYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 286 TH2F* pZYCResidualsHistBnil = new TH2F("fZYCResidualsHistBnil","z-r\\phi residuals side C (z<0), Bnil", 100, -3., 3., 100, -1.5, 1.5 );
51cfc50d 287 pZYCResidualsHistBnil->SetXTitle("\\deltaz [cm]");
288 pZYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 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 );
51cfc50d 290 pLPAResidualsHistBnil->SetXTitle("\\deltasin(\\phi)");
291 pLPAResidualsHistBnil->SetYTitle("\\deltatan(\\lambda)");
d6d927ac 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 );
51cfc50d 293 pLPCResidualsHistBnil->SetXTitle("\\deltasin(\\phi)");
294 pLPCResidualsHistBnil->SetYTitle("\\deltatan(\\lambda)");
d6d927ac 295 TH2F* pPhiYAResidualsHistBnil = new TH2F("fPhiYAResidualsHistBnil","\\phi-\\deltar\\phi side A (z>0), Bnil",36,0.0,6.3,100,-1.5,1.5);
51cfc50d 296 pPhiYAResidualsHistBnil->SetXTitle("\\phi [rad]");
297 pPhiYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 298 TH2F* pPhiYCResidualsHistBnil = new TH2F("fPhiYCResidualsHistBnil","\\phi-\\deltar\\phi side C (z<0), Bnil",36,0.0,6.3,100,-1.5,1.5);
51cfc50d 299 pPhiYCResidualsHistBnil->SetXTitle("\\phi [rad]");
300 pPhiYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 301 TH2F* pPhiZAResidualsHistBnil = new TH2F("fPhiZAResidualsHistBnil","\\phi-\\deltaz side A (z>0), Bnil",36,0.0,6.3,100,-3.0,3.0);
51cfc50d 302 pPhiZAResidualsHistBnil->SetXTitle("\\phi [rad]");
303 pPhiZAResidualsHistBnil->SetYTitle("\\deltaz [cm]");
d6d927ac 304 TH2F* pPhiZCResidualsHistBnil = new TH2F("fPhiZCResidualsHistBnil","\\Pt-\\deltaz side C (z<0), Bnil",36,0.0,6.3,100,-3.0,3.0);
51cfc50d 305 pPhiZCResidualsHistBnil->SetXTitle("\\phi [rad]");
306 pPhiZCResidualsHistBnil->SetYTitle("\\deltaz [cm]");
d6d927ac 307 TH2F* pPtYAResidualsHistBnil = new TH2F("fPtYAResidualsHistBnil","Pt-\\deltar\\phi side A (z>0), Bnil",20,.3,8.,80,-1.5,1.5);
51cfc50d 308 pPtYAResidualsHistBnil->SetXTitle("Pt [rad]");
309 pPtYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 310 TH2F* pPtYCResidualsHistBnil = new TH2F("fPtYCResidualsHistBnil","Pt-\\deltar\\phi side C (z<0), Bnil",20,.3,8.,80,-1.5,1.5);
51cfc50d 311 pPtYCResidualsHistBnil->SetXTitle("Pt [rad]");
312 pPtYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 313 TH2F* pPtZAResidualsHistBnil = new TH2F("fPtZAResidualsHistBnil","Pt-\\deltaz side A (z>0), Bnil",20,.3,8.,80,-3.0,3.0);
51cfc50d 314 pPtZAResidualsHistBnil->SetXTitle("Pt");
315 pPtZAResidualsHistBnil->SetYTitle("\\deltaz [cm]");
d6d927ac 316 TH2F* pPtZCResidualsHistBnil = new TH2F("fPtZCResidualsHistBnil","Pt-\\deltaz side C (z<0), Bnil",20,.3,8.,80,-3.0,3.0);
51cfc50d 317 pPtZCResidualsHistBnil->SetXTitle("Pt");
318 pPtZCResidualsHistBnil->SetYTitle("\\deltaz [cm]");
d6d927ac 319 TH2F* pLowPtYAResidualsHistBnil = new TH2F("fLowPtYAResidualsHistBnil","Pt-\\deltar\\phi side A (z>0), Bnil",100,0.0,1.5,80,-1.5,1.5);
51cfc50d 320 pLowPtYAResidualsHistBnil->SetXTitle("Pt [rad]");
321 pLowPtYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 322 TH2F* pLowPtYCResidualsHistBnil = new TH2F("fLowPtYCResidualsHistBnil","Pt-\\deltar\\phi side C (z<0), Bnil",100,0.0,1.5,80,-1.5,1.5);
51cfc50d 323 pLowPtYCResidualsHistBnil->SetXTitle("Pt [rad]");
324 pLowPtYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
d6d927ac 325 TH2F* pLowPtZAResidualsHistBnil = new TH2F("fLowPtZAResidualsHistBnil","Pt-\\deltaz side A (z>0), Bnil",100,0.0,1.5,80,-3.0,3.0);
51cfc50d 326 pLowPtZAResidualsHistBnil->SetXTitle("Pt");
327 pLowPtZAResidualsHistBnil->SetYTitle("\\deltaz [cm]");
d6d927ac 328 TH2F* pLowPtZCResidualsHistBnil = new TH2F("fLowPtZCResidualsHistBnil","Pt-\\deltaz side C (z<0), Bnil",100,0.0,1.5,80,-3.0,3.0);
51cfc50d 329 pLowPtZCResidualsHistBnil->SetXTitle("Pt");
330 pLowPtZCResidualsHistBnil->SetYTitle("\\deltaz [cm]");
331 TList* listBnil = new TList();
63f4d30f 332 listBnil->SetName("B0");
333 listBnil->SetOwner(kTRUE);
d6d927ac 334 fList->Add(listBnil); //2
51cfc50d 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);
63f4d30f 351
51cfc50d 352 TH1F* pNmatchingEff=new TH1F("pNmatchingEff","matching efficiency",50,0.,1.);
d6d927ac 353 fList->Add(pNmatchingEff); //3
354
63f4d30f 355 TH1I* pChecks = new TH1I("some counters","some counters",10,0,10);
d6d927ac 356 pChecks->GetXaxis()->SetBinLabel(1+kNoESD,"no ESD");
357 pChecks->GetXaxis()->SetBinLabel(1+kNoESDfriend,"no ESDfriend");
358 pChecks->GetXaxis()->SetBinLabel(1+kNoFriendTrack,"no friendtrack");
63f4d30f 359 pChecks->GetXaxis()->SetBinLabel(1+kNoITSoutParams,"no ITS out params");
360 pChecks->GetXaxis()->SetBinLabel(1+kESDfriend,"ESD friend");
d6d927ac 361 pChecks->GetXaxis()->SetBinLabel(1+kFriendsSkipBit,"friends+skipbit");
362 pChecks->GetXaxis()->SetBinLabel(1+kFriendTrack,"friendtrack");
63f4d30f 363 pChecks->GetXaxis()->SetBinLabel(1+kITSoutParams,"ITS out params");
d6d927ac 364 fList->Add(pChecks); //4
51cfc50d 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);
63f4d30f 375
376 // Post output data.
377 PostData(0, fDebugTree);
378 PostData(1, fList);
379 PostData(2, fArrayITSglobal);
380 PostData(3, fArrayITSsa);
4a84c20d 381}
382
383//________________________________________________________________________
d6d927ac 384void AliAnalysisTaskITSTPCalignment::UserExec(Option_t *)
4a84c20d 385{
386 // Main loop
387 // Called for each event
51cfc50d 388
d6d927ac 389 TH1I* pChecks = static_cast<TH1I*>(fList->At(4));
51cfc50d 390 //check for ESD and friends
d6d927ac 391 AliESDEvent* event = dynamic_cast<AliESDEvent*>(InputEvent());
392 AliESDfriend* eventFriend = dynamic_cast<AliESDfriend*>(ESDfriend());
393 if (!event)
4a84c20d 394 {
51cfc50d 395 AliError("no ESD");
d6d927ac 396 pChecks->Fill(kNoESD);
4a84c20d 397 return;
398 }
399
d6d927ac 400 if (!eventFriend)
4a84c20d 401 {
51cfc50d 402 AliError("no ESD friend");
d6d927ac 403 pChecks->Fill(kNoESDfriend);
51cfc50d 404 return;
4a84c20d 405 }
d6d927ac 406 else
407 {
408 pChecks->Fill(kESDfriend);
409 }
4a84c20d 410
63f4d30f 411 //event->SetESDfriend(eventFriend);
412
413 if (eventFriend->TestSkipBit())
414 {
415 pChecks->Fill(kFriendsSkipBit);
416 return;
417 }
51cfc50d 418
419 //Update the parmeters
d6d927ac 420 AnalyzeESDevent(event);
51cfc50d 421}
422
423//________________________________________________________________________
d6d927ac 424void AliAnalysisTaskITSTPCalignment::AnalyzeESDevent(AliESDEvent* event)
51cfc50d 425{
426 //analyze an ESD event with track matching
d6d927ac 427 Int_t ntracks = event->GetNumberOfTracks();
51cfc50d 428 if (ntracks==0) return;
429
430 if (fUseITSoutGlobalTrack)
4a84c20d 431 {
d6d927ac 432 AliRelAlignerKalman* alignerGlobal = fArrayITSglobal->GetAligner(event);
51cfc50d 433 if (alignerGlobal)
d6d927ac 434 alignerGlobal->AddESDevent(event);
4a84c20d 435 }
436
51cfc50d 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
d6d927ac 444 Int_t n = FindMatchingTracks(arrayParamITS, arrayParamTPC, event);
63f4d30f 445 AliInfo(Form("matched tracklet pairs: %i\n",n));
51cfc50d 446
9a136fa5 447 TH1F* nMatchingEff=static_cast<TH1F*>(fList->At(3));
51cfc50d 448 Float_t ratio = (float)n/(float)ntracks;
449 nMatchingEff->Fill(ratio);
450
451 AliRelAlignerKalman* alignerSA = NULL;
d6d927ac 452 if (n>0) alignerSA = fArrayITSsa->GetAligner(event);
51cfc50d 453
454 for (Int_t i=0;i<n;i++)
4a84c20d 455 {
9a136fa5 456 AliExternalTrackParam* paramsITS=static_cast<AliExternalTrackParam*>(arrayParamITS[i]);
457 AliExternalTrackParam* paramsTPC=static_cast<AliExternalTrackParam*>(arrayParamTPC[i]);
458
459 if (!(paramsITS&&paramsTPC)) continue;
51cfc50d 460
461 //QA
462 if (fDoQA) DoQA(paramsITS,paramsTPC);
463
464 //debug
465 if (fAligner->AddTrackParams(paramsITS, paramsTPC))
4a84c20d 466 {
d6d927ac 467 fAligner->SetRunNumber(event->GetRunNumber());
468 fAligner->SetMagField(event->GetMagneticField());
469 fAligner->SetTimeStamp(event->GetTimeStamp());
51cfc50d 470 }
4a84c20d 471
63f4d30f 472 //alignment
51cfc50d 473 if (alignerSA) alignerSA->AddTrackParams(paramsITS, paramsTPC);
474 }
475 arrayParamITS.Delete();
476 arrayParamTPC.Delete();
477 }
4a84c20d 478}
479
480//________________________________________________________________________
481void AliAnalysisTaskITSTPCalignment::Terminate(Option_t *)
482{
483 // Called once at the end of the query
63f4d30f 484 fArrayITSsa->Print("a");
485 fAligner->Print();
51cfc50d 486}
487
488//________________________________________________________________________
489Int_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
d6d927ac 494 TH1I* pChecks = static_cast<TH1I*>(fList->At(4));
51cfc50d 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();
d6d927ac 517 if (!constfriendtrack1) {AliInfo(Form("no friendTrack1\n"));pChecks->Fill(kNoFriendTrack);continue;}
518 pChecks->Fill(kFriendTrack);
51cfc50d 519 AliESDfriendTrack friendtrack1(*constfriendtrack1);
520
d6d927ac 521 if (!friendtrack1.GetITSOut()) {pChecks->Fill(kNoITSoutParams);continue;}
522 pChecks->Fill(kITSoutParams);
51cfc50d 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
9a136fa5 599 delete [] matchedArr;
51cfc50d 600 return iMatched+1;
601}
602
603//________________________________________________________________________
604void 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;
d6d927ac 614 Double_t field = fInputEvent->GetMagneticField();
9a136fa5 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));
51cfc50d 618 if (!pList) return;
619
9a136fa5 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));
51cfc50d 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 }
4a84c20d 659}