Updating ITS macros (Ruben)
[u/mrichter/AliRoot.git] / PWGPP / ITS / AliAnalysisTaskITSAlignQA.cxx
CommitLineData
89b48e9e 1#include "AliAnalysisTask.h"
2#include "AliAnalysisManager.h"
3#include "AliAnalysisDataContainer.h"
4#include "AliITSRecPoint.h"
5#include "AliESDEvent.h"
6#include "AliTrackPointArray.h"
7#include "AliITSgeomTGeo.h"
8#include "AliITSTPArrayFit.h"
9#include "AliESDfriend.h"
10#include "AliCDBManager.h"
11#include "AliCDBEntry.h"
12#include "AliITSCalibrationSDD.h"
13#include "AliITSresponseSDD.h"
14#include "AliGeomManager.h"
00a2ed4c 15#include "AliMultiplicity.h"
89b48e9e 16#include <TSystem.h>
17#include <TTree.h>
18#include <TH1F.h>
19#include <TH2F.h>
f09774b4 20#include <TProfile.h>
89b48e9e 21#include <TChain.h>
22#include <TGeoGlobalMagField.h>
6a5f5ccc 23#include "AliAODHandler.h"
89b48e9e 24#include "AliESDInputHandlerRP.h"
6a5f5ccc 25#include "AliITSSumTP.h"
26#include "AliMagF.h"
89b48e9e 27
28/**************************************************************************
29 * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
30 * *
31 * Author: The ALICE Off-line Project. *
32 * Contributors are mentioned in the code where appropriate. *
33 * *
34 * Permission to use, copy, modify and distribute this software and its *
35 * documentation strictly for non-commercial purposes is hereby granted *
36 * without fee, provided that the above copyright notice appears in all *
37 * copies and that both the copyright notice and this permission notice *
38 * appear in the supporting documentation. The authors make no claims *
39 * about the suitability of this software for any purpose. It is *
40 * provided "as is" without express or implied warranty. *
41 **************************************************************************/
42
43//*************************************************************************
44// Implementation of class AliAnalysiTaskITSAlignQA
45// AliAnalysisTaskSE to extract from ESD + ESDfriends
46// the track-to-point residuals and dE/dx vs, time for SDD modules
47//
48// Author: F. Prino, prino@to.infn.it
49//*************************************************************************
50
51
52#include "AliAnalysisTaskITSAlignQA.h"
53
54ClassImp(AliAnalysisTaskITSAlignQA)
55//______________________________________________________________________________
56AliAnalysisTaskITSAlignQA::AliAnalysisTaskITSAlignQA() : AliAnalysisTaskSE("SDD Calibration"),
57 fOutput(0),
58 fHistNEvents(0),
59 fHistPtAccept(0),
60 fDoSPDResiduals(kTRUE),
61 fDoSDDResiduals(kTRUE),
62 fDoSSDResiduals(kTRUE),
63 fDoSDDdEdxCalib(kTRUE),
826883ac 64 fDoSDDVDriftCalib(kTRUE),
65 fDoSDDDriftTime(kTRUE),
6a5f5ccc 66 fDoFillTPTree(kFALSE),
89b48e9e 67 fUseITSsaTracks(kFALSE),
9ea163ca 68 fLoadGeometry(kFALSE),
df69bb07 69 fUseVertex(kFALSE),
70 fUseVertexForZOnly(kFALSE),
6a5f5ccc 71 fUseTPCMomentum(kFALSE),
df69bb07 72 fMinVtxContributors(5),
5104b23e 73 fRemovePileupWithSPD(kTRUE),
89b48e9e 74 fMinITSpts(3),
75 fMinTPCpts(70),
76 fMinPt(0.5),
77 fNPtBins(8),
00a2ed4c 78 fMinMult(0),
79 fMaxMult(1e9),
6a5f5ccc 80 fCutDCAXY(1.e10),
81 fCutDCAZ(1.e10),
89b48e9e 82 fFitter(0),
6a5f5ccc 83 fITSSumTP(),
89b48e9e 84 fRunNb(0),
85 fOCDBLocation("local://$ALICE_ROOT/OCDB")
86{
87 //
88 fFitter = new AliITSTPArrayFit(5);
89 Double_t xbins[9]={0.3,0.5,0.75,1.,1.5,2.,3.,5.,10.};
90 SetPtBinLimits(8,xbins);
91 DefineOutput(1, TList::Class());
92}
93
94
95//___________________________________________________________________________
96AliAnalysisTaskITSAlignQA::~AliAnalysisTaskITSAlignQA(){
97 //
6a5f5ccc 98 if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutput;
99 delete fFitter;
100 delete fITSSumTP;
101 //
89b48e9e 102}
103//___________________________________________________________________________
104void AliAnalysisTaskITSAlignQA::UserCreateOutputObjects() {
105 //
9ea163ca 106
107 if(fLoadGeometry) LoadGeometryFromOCDB();
89b48e9e 108
109 fOutput = new TList();
110 fOutput->SetOwner();
111 fOutput->SetName("OutputHistos");
112
00a2ed4c 113 fHistNEvents = new TH1F("hNEvents", "Number of processed events",kNEvStatBins,-0.5,kNEvStatBins-0.5);
89b48e9e 114 fHistNEvents->Sumw2();
115 fHistNEvents->SetMinimum(0);
00a2ed4c 116 fHistNEvents->GetXaxis()->SetBinLabel(kEvAll+1,"All Events");
117 fHistNEvents->GetXaxis()->SetBinLabel(kEvCnt+1,"After Centrality cut");
118 fHistNEvents->GetXaxis()->SetBinLabel(kEvVtx+1,"After Vertex cut");
119 fHistNEvents->GetXaxis()->SetBinLabel(kEvPlp+1,"After Pileup cut");
120 fHistNEvents->GetXaxis()->SetBinLabel(kNTracks+1,"Tracks Accepted");
89b48e9e 121 fOutput->Add(fHistNEvents);
122
123 fHistPtAccept = new TH1F("hPtAccept","Pt distrib of accepted tracks",50,0.,5.);
124 fHistPtAccept->Sumw2();
125 fHistPtAccept->SetMinimum(0);
126 fOutput->Add(fHistPtAccept);
127
128 if(fDoSPDResiduals) CreateSPDHistos();
826883ac 129 if(fDoSDDResiduals || fDoSDDdEdxCalib || fDoSDDVDriftCalib || fDoSDDDriftTime) CreateSDDHistos();
89b48e9e 130 if(fDoSSDResiduals) CreateSSDHistos();
6a5f5ccc 131 //
132 if (fDoFillTPTree) {
133 fITSSumTP = new AliITSSumTP();
134 AliAODHandler* handler = dynamic_cast<AliAODHandler*>( AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler() );
135 if (!handler) AliFatal("TP tree requested but AOD handler is not available");
136 handler->AddBranch("AliITSSumTP",&fITSSumTP);
137 handler->SetFillAOD(kFALSE); // manual fill
138 CreateUserInfo();
139 //
140 }
141 //
9ea163ca 142 PostData(1,fOutput);
89b48e9e 143}
144
145//___________________________________________________________________________
146void AliAnalysisTaskITSAlignQA::CreateSPDHistos(){
147 // Histos for SPD
148
149
150 for(Int_t iMod=0; iMod<kNSPDmods; iMod++){
151 fHistSPDResidX[iMod] = new TH2F(Form("hSPDResidX%d",iMod),
152 Form("hSPDResidX%d",iMod),
153 fNPtBins,fPtBinLimits,
154 250,-0.05,0.05);
155 fHistSPDResidX[iMod]->Sumw2();
156 fOutput->Add(fHistSPDResidX[iMod]);
157
158 fHistSPDResidZ[iMod] = new TH2F(Form("hSPDResidZ%d",iMod),
159 Form("hSPDResidZ%d",iMod),
160 fNPtBins,fPtBinLimits,
161 250,-0.1,0.1);
162 fHistSPDResidZ[iMod]->Sumw2();
163 fOutput->Add(fHistSPDResidZ[iMod]);
164 }
165 return;
166}
167
168//___________________________________________________________________________
169void AliAnalysisTaskITSAlignQA::CreateSDDHistos(){
170 // Histos for SDD
171
172 for(Int_t iMod=0; iMod<kNSDDmods; iMod++){
826883ac 173 if (fDoSDDResiduals) {
f09774b4 174 fHistSDDResidX[iMod] = new TH2F(Form("hSDDResidX%d",iMod+kNSPDmods),
175 Form("hSDDResidX%d",iMod+kNSPDmods),
89b48e9e 176 fNPtBins,fPtBinLimits,
177 300,-0.15,0.15);
178 fHistSDDResidX[iMod]->Sumw2();
179 fOutput->Add(fHistSDDResidX[iMod]);
180
f09774b4 181 fHistSDDResidZ[iMod] = new TH2F(Form("hSDDResidZ%d",iMod+kNSPDmods),
182 Form("hSDDResidZ%d",iMod+kNSPDmods),
89b48e9e 183 fNPtBins,fPtBinLimits,
184 200,-0.1,0.1);
185 fHistSDDResidZ[iMod]->Sumw2();
186 fOutput->Add(fHistSDDResidZ[iMod]);
187
188 fHistSDDResidXvsX[iMod] = new TH2F(Form("hSDDResidXvsX%d",iMod+kNSPDmods),
189 Form("hSDDResidXvsX%d",iMod+kNSPDmods),
190 40,-3.5,3.5,300,-0.15,0.15);
191 fHistSDDResidXvsX[iMod]->Sumw2();
192 fOutput->Add(fHistSDDResidXvsX[iMod]);
193
194 fHistSDDResidXvsZ[iMod] = new TH2F(Form("hSDDResidXvsZ%d",iMod+kNSPDmods),
195 Form("hSDDResidXvsZ%d",iMod+kNSPDmods),
196 10,-3.8,3.8,300,-0.15,0.15);
197 fHistSDDResidXvsZ[iMod]->Sumw2();
198 fOutput->Add(fHistSDDResidXvsZ[iMod]);
199
200 fHistSDDResidZvsX[iMod] = new TH2F(Form("hSDDResidZvsX%d",iMod+kNSPDmods),
201 Form("hSDDResidZvsX%d",iMod+kNSPDmods),
202 40,-3.5,3.5,200,-0.1,0.1);
203 fHistSDDResidZvsX[iMod]->Sumw2();
204 fOutput->Add(fHistSDDResidZvsX[iMod]);
205
206 fHistSDDResidZvsZ[iMod] = new TH2F(Form("hSDDResidZvsZ%d",iMod+kNSPDmods),
207 Form("hSDDResidZvsZ%d",iMod+kNSPDmods),
208 10,-3.8,3.8,200,-0.1,0.1);
209 fHistSDDResidZvsZ[iMod]->Sumw2();
210 fOutput->Add(fHistSDDResidZvsZ[iMod]);
f09774b4 211 //
826883ac 212 }
213 //
214 if (fDoSDDVDriftCalib) {
f09774b4 215 for (int ix=0;ix<2;ix++) { // profile histos per side
216 //
217 char* hnm = Form("hpSDDResXvsXD%d_%d",iMod+kNSPDmods,ix);
218 int nbX = 50, nbZ = 20, nbXOffs = 2, nbZOffs = 2;
219 double xRange = 3.5085e4, zRange = 7.5264e4, xOffs = nbXOffs*xRange/nbX, zOffs = nbZOffs*zRange/nbZ;
220 fHProfSDDResidXvsXD[iMod][ix] = new TProfile(hnm, hnm, nbX+2*nbXOffs, -xOffs, xRange + xOffs);
221 fHProfSDDResidXvsXD[iMod][ix]->Sumw2();
222 fOutput->Add(fHProfSDDResidXvsXD[iMod][ix]);
223 //
224 hnm = Form("hpSDDDrTimevsXD%d_%d",iMod+kNSPDmods,ix);
225 fHProfSDDDrTimevsXD[iMod][ix] = new TProfile(hnm, hnm, nbX+2*nbXOffs, -xOffs, xRange + xOffs);
226 fHProfSDDDrTimevsXD[iMod][ix]->Sumw2();
227 fOutput->Add(fHProfSDDDrTimevsXD[iMod][ix]);
228 //
229 hnm = Form("hpSDDResXvsZ%d_%d",iMod+kNSPDmods,ix);
230 fHProfSDDResidXvsZ[iMod][ix] = new TProfile(hnm, hnm, nbZ+2*nbZOffs, -(0.5*zRange+zOffs),(0.5*zRange+zOffs));
231 fHProfSDDResidXvsZ[iMod][ix]->Sumw2();
232 fOutput->Add(fHProfSDDResidXvsZ[iMod][ix]);
233 //
234 hnm = Form("hpSDDDrTimevsZ%d_%d",iMod+kNSPDmods,ix);
235 fHProfSDDDrTimevsZ[iMod][ix] = new TProfile(hnm, hnm, nbZ+2*nbZOffs, -(0.5*zRange+zOffs),(0.5*zRange+zOffs));
236 fHProfSDDDrTimevsZ[iMod][ix]->Sumw2();
237 fOutput->Add(fHProfSDDDrTimevsZ[iMod][ix]);
238 //
239 }
89b48e9e 240 }
826883ac 241
89b48e9e 242 if(fDoSDDdEdxCalib){
243 fHistSDDdEdxvsDrTime[iMod] = new TH2F(Form("hSDDdEdxvsDrTime%d",iMod+kNSPDmods),
244 Form("hSDDdEdxvsDrTime%d",iMod+kNSPDmods),
245 16,0.,6400.,100,0.,300.);
246 fHistSDDdEdxvsDrTime[iMod]->Sumw2();
247 fOutput->Add(fHistSDDdEdxvsDrTime[iMod]);
248 }
826883ac 249 //
250 if (fDoSDDDriftTime) {
251 fHistSDDDrTimeAll[iMod]=new TH1F(Form("hSDDDrTimeAll%d",iMod+kNSPDmods),
252 Form("hSDDDrTimeAll%d",iMod+kNSPDmods),
89b48e9e 253 3200,0.,6400.);
826883ac 254 fHistSDDDrTimeAll[iMod]->Sumw2();
255 fHistSDDDrTimeAll[iMod]->SetMinimum(0.);
256 fOutput->Add(fHistSDDDrTimeAll[iMod]);
257
258 fHistSDDDrTimeExtra[iMod]=new TH1F(Form("hSDDDrTimeExtra%d",iMod+kNSPDmods),
259 Form("hSDDDrTimeExtra%d",iMod+kNSPDmods),
260 3200,0.,6400.);
261 fHistSDDDrTimeExtra[iMod]->Sumw2();
262 fHistSDDDrTimeExtra[iMod]->SetMinimum(0.);
263 fOutput->Add(fHistSDDDrTimeExtra[iMod]);
264
265 fHistSDDDrTimeAttac[iMod]=new TH1F(Form("hSDDDrTimeAttac%d",iMod+kNSPDmods),
266 Form("hSDDDrTimeAttac%d",iMod+kNSPDmods),
267 3200,0.,6400.);
268 fHistSDDDrTimeAttac[iMod]->Sumw2();
269 fHistSDDDrTimeAttac[iMod]->SetMinimum(0.);
270 fOutput->Add(fHistSDDDrTimeAttac[iMod]);
271 }
89b48e9e 272 }
273 return;
826883ac 274 //
89b48e9e 275}
826883ac 276
89b48e9e 277//___________________________________________________________________________
278void AliAnalysisTaskITSAlignQA::CreateSSDHistos(){
279 // Histos for SSD
280 for(Int_t iMod=0; iMod<kNSSDmods; iMod++){
281 fHistSSDResidX[iMod] = new TH2F(Form("hSSDResidX%d",iMod+kNSPDmods+kNSDDmods),
282 Form("hSSDResidX%d",iMod+kNSPDmods+kNSDDmods),
283 fNPtBins,fPtBinLimits,
284 250,-0.1,0.1);
285 fHistSSDResidX[iMod]->Sumw2();
286 fOutput->Add(fHistSSDResidX[iMod]);
287
288 fHistSSDResidZ[iMod] = new TH2F(Form("hSSDResidZ%d",iMod+kNSPDmods+kNSDDmods),
289 Form("hSSDResidZ%d",iMod+kNSPDmods+kNSDDmods),
290 fNPtBins,fPtBinLimits,
291 250,-1.,1.);
292 fHistSSDResidZ[iMod]->Sumw2();
293 fOutput->Add(fHistSSDResidZ[iMod]);
294 }
295 return;
296}
297//______________________________________________________________________________
298void AliAnalysisTaskITSAlignQA::UserExec(Option_t *)
299{
df69bb07 300 //
301 static AliTrackPointArray* arrayITS = 0;
6a5f5ccc 302 AliTrackPointArray* arrayITSNoVtx = 0;
89b48e9e 303 //
304 AliESDEvent *esd = (AliESDEvent*) (InputEvent());
6a5f5ccc 305 if (fITSSumTP) fITSSumTP->Reset();
89b48e9e 306
307 if(!esd) {
308 printf("AliAnalysisTaskITSAlignQA::Exec(): bad ESD\n");
309 return;
310 }
311
89b48e9e 312 if(!ESDfriend()) {
313 printf("AliAnalysisTaskITSAlignQA::Exec(): bad ESDfriend\n");
314 return;
315 }
df69bb07 316 //
00a2ed4c 317 if (!AcceptCentrality(esd)) return;
318 fHistNEvents->Fill(kEvCnt);
319
5104b23e 320 const AliESDVertex* vtx=0,*vtxSPD=0;
00a2ed4c 321 fHistNEvents->Fill(kEvAll);
6a5f5ccc 322 vtx = esd->GetPrimaryVertex();
323 vtxSPD = esd->GetPrimaryVertexSPD();
324 //
df69bb07 325 if (fUseVertex) { // check the vertex if it is requested as an extra point
5104b23e 326 if (!AcceptVertex(vtx,vtxSPD)) return;
327 }
00a2ed4c 328
329 fHistNEvents->Fill(kEvVtx);
5104b23e 330 if (fRemovePileupWithSPD){
331 // skip events tagged by SPD as pileup
332 if(esd->IsPileupFromSPD()) return;
df69bb07 333 }
00a2ed4c 334 fHistNEvents->Fill(kEvPlp);
5104b23e 335
df69bb07 336 //
89b48e9e 337 fFitter->SetBz(esd->GetMagneticField());
338
339 const AliTrackPointArray *array = 0;
340 Int_t ntracks = esd->GetNumberOfTracks();
341
342 for (Int_t itrack=0; itrack < ntracks; itrack++) {
df69bb07 343 //
344 if (arrayITS) {delete arrayITS; arrayITS = 0;} // reset points from previous tracks
6a5f5ccc 345 arrayITSNoVtx = 0;
df69bb07 346 //
89b48e9e 347 AliESDtrack * track = esd->GetTrack(itrack);
348 if(!track) continue;
6a5f5ccc 349 if(!AcceptTrack(track, vtx)) continue;
89b48e9e 350 array = track->GetTrackPointArray();
351 if(!array) continue;
df69bb07 352 arrayITS = PrepareTrack(array, vtx);
6a5f5ccc 353 if (fITSSumTP) {
354 arrayITSNoVtx = PrepareTrack(array, 0);
355 arrayITSNoVtx->SetUniqueID(itrack);
356 fITSSumTP->AddTrack(arrayITSNoVtx);
357 }
df69bb07 358 //
00a2ed4c 359 fHistNEvents->Fill(kNTracks);
360 //
df69bb07 361 Int_t npts = arrayITS->GetNPoints();
362 Int_t npts1 = fUseVertexForZOnly ? npts-1 : npts;
363 //
89b48e9e 364 if(fDoSPDResiduals){
df69bb07 365 FitAndFillSPD(1,arrayITS,npts1,track);
366 FitAndFillSPD(2,arrayITS,npts1,track);
89b48e9e 367 }
826883ac 368 if(fDoSDDResiduals || fDoSDDdEdxCalib || fDoSDDVDriftCalib || fDoSDDDriftTime) {
df69bb07 369 FitAndFillSDDrphi(arrayITS,npts,track);
826883ac 370 if (fDoSDDResiduals) {
371 FitAndFillSDDz(3,arrayITS,npts1,track);
372 FitAndFillSDDz(4,arrayITS,npts1,track);
373 }
9ea163ca 374 }
89b48e9e 375 if(fDoSSDResiduals){
df69bb07 376 FitAndFillSSD(5,arrayITS,npts1,track);
377 FitAndFillSSD(6,arrayITS,npts1,track);
89b48e9e 378 }
379 }
6a5f5ccc 380 //
381 if (fITSSumTP) { // store vertex and mometum info
382 fITSSumTP->SetVertex(vtx);
383 TObjArray& tps = fITSSumTP->GetTracks();
384 int ntp = tps.GetEntriesFast();
385 fITSSumTP->BookNTracks(ntp);
386 for (int it=ntp;it--;) {
387 AliTrackPointArray* tp = (AliTrackPointArray*)tps[it];
388 if (!tp) continue;
389 AliESDtrack* esdTr = esd->GetTrack(tp->GetUniqueID());
390 double crv = esdTr->GetC(esd->GetMagneticField());
391 double crve = TMath::Sqrt(esdTr->GetSigma1Pt2()) * esd->GetMagneticField()*kB2C;
392 fITSSumTP->SetCrvGlo(it,crv);
393 fITSSumTP->SetCrvGloErr(it,crve);
394 const AliExternalTrackParam* inTPC = esdTr->GetTPCInnerParam();
395 if (inTPC) {
396 crv = inTPC->GetC(esd->GetMagneticField());
397 crve = TMath::Sqrt(inTPC->GetSigma1Pt2()) * esd->GetMagneticField()*kB2C;
398 fITSSumTP->SetCrvTPC(it,crv);
399 fITSSumTP->SetCrvTPCErr(it,crve);
400 }
401 }
402 fITSSumTP->SetUniqueID(fCurrentRunNumber);
403 AliAODHandler* handler = dynamic_cast<AliAODHandler*>( AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler() );
404 if (!ntp) handler->FillTree();
405 }
89b48e9e 406
6a5f5ccc 407 //
89b48e9e 408 PostData(1,fOutput);
409
410}
df69bb07 411
89b48e9e 412//___________________________________________________________________________
6a5f5ccc 413Bool_t AliAnalysisTaskITSAlignQA::AcceptTrack(const AliESDtrack * track, const AliESDVertex* vtx)
414{
89b48e9e 415 // track selection cuts
416 Bool_t accept=kTRUE;
417 if(fUseITSsaTracks){
418 if(track->GetNcls(1)>0) accept=kFALSE;
419 }else{
420 if(track->GetNcls(1)<fMinTPCpts) accept=kFALSE;
421 }
422 if(track->GetNcls(0) < fMinITSpts) accept=kFALSE;
423 Int_t trstatus=track->GetStatus();
424 if(!(trstatus&AliESDtrack::kITSrefit)) accept=kFALSE;
6a5f5ccc 425 Float_t pt = 0;
426 if (fUseTPCMomentum) {
427 if (track->GetTPCInnerParam()) pt = track->GetTPCInnerParam()->Pt();
428 else pt = track->Pt();
429 }
89b48e9e 430 if(pt<fMinPt) accept=kFALSE;
6a5f5ccc 431 //
432 // if vertex constraint is used, apply soft DCA cut
433 if (vtx) {
434 Double_t dz[2],cov[3];
435 AliExternalTrackParam trc = *track;
436 if (!trc.PropagateToDCA(vtx, fFitter->GetBz(), 3.0, dz, cov)) accept = kFALSE;
437 else {
438 if (dz[0]*dz[0]/(1e-4+cov[0])>fCutDCAXY) accept = kFALSE;
439 if (dz[1]*dz[1]/(4e-4+cov[2])>fCutDCAZ) accept = kFALSE;
440 }
441 }
442 //
89b48e9e 443 if(accept) fHistPtAccept->Fill(pt);
444 return accept;
445}
df69bb07 446
447//___________________________________________________________________________
5104b23e 448Bool_t AliAnalysisTaskITSAlignQA::AcceptVertex(const AliESDVertex * vtx, const AliESDVertex * vtxSPD) {
df69bb07 449 // vertex selection cuts
5104b23e 450 if (!vtx || vtx->GetStatus()<1) return kFALSE;
451 if (!vtxSPD || vtxSPD->GetStatus()<1) return kFALSE;
df69bb07 452 if (vtx->GetNContributors()<fMinVtxContributors) return kFALSE;
5104b23e 453 if (TMath::Abs(vtx->GetZ()-vtxSPD->GetZ())>0.3) return kFALSE;
df69bb07 454 return kTRUE;
455}
456
89b48e9e 457//___________________________________________________________________________
458void AliAnalysisTaskITSAlignQA::FitAndFillSPD(Int_t iLayer, const AliTrackPointArray *array, Int_t npts,AliESDtrack * track){
459 // fit track and fills histos for SPD
460 fFitter->AttachPoints(array,0, npts-1);
461 Int_t iPtSPD[4],modIdSPD[4];
462 Int_t nPtSPD=0;
463 Double_t resGlo[3],resLoc[3];
464 for(Int_t ipt=0; ipt<npts; ipt++) {
465 AliTrackPoint point;
466 Int_t modId;
467 array->GetPoint(point,ipt);
468 Int_t volId = point.GetVolumeID();
df69bb07 469 if (volId == kVtxSensVID) continue; // this is a vertex constraint
89b48e9e 470 Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId);
471 if(layerId==iLayer){
472 modId+=AliITSgeomTGeo::GetModuleIndex(layerId,1,1);
473 iPtSPD[nPtSPD] = ipt;
474 modIdSPD[nPtSPD] = modId;
475 ++nPtSPD;
476 fFitter->SetCovIScale(ipt,1e-4);
477 }
478 }
479 if(nPtSPD>0){
6a5f5ccc 480 double pt = fUseTPCMomentum ? track->GetTPCInnerParam()->Pt() : track->Pt();
481 fFitter->Fit(track->Charge(),pt,0.);
89b48e9e 482 Double_t chi2=fFitter->GetChi2NDF();
483 if ( chi2<0 || chi2>1e4 ) return; // fit failed, abandon this track
484 for (Int_t ip=0; ip<nPtSPD;ip++) {
485 fFitter->GetResiduals(resGlo,iPtSPD[ip]);
486 TGeoHMatrix *mcurr = AliITSgeomTGeo::GetMatrix(modIdSPD[ip]);
487 mcurr->MasterToLocalVect(resGlo,resLoc);
488 Int_t index=modIdSPD[ip];
6a5f5ccc 489 fHistSPDResidX[index]->Fill(pt,resLoc[0]);
490 fHistSPDResidZ[index]->Fill(pt,resLoc[2]);
89b48e9e 491 }
492 }
493}
494//___________________________________________________________________________
9ea163ca 495void AliAnalysisTaskITSAlignQA::FitAndFillSDDrphi(const AliTrackPointArray *array, Int_t npts, AliESDtrack * track){
496 // fit track and fills histos for SDD along rphi (drift coord.)
89b48e9e 497 Double_t dedx[4];
498 track->GetITSdEdxSamples(dedx);
499
500 fFitter->AttachPoints(array,0, npts-1);
f09774b4 501 Int_t iPtSDD[4],modIdSDD[4],modSide[4];
502 Double_t xLocSDD[4],zLocSDD[4],drTime[4];
89b48e9e 503 Int_t nPtSDD=0;
504 Int_t nPtSSDSPD=0;
505 Double_t resGlo[3],resLoc[3];
506 Float_t posGloF[3];
507 Double_t posGlo[3],posLoc[3];
f09774b4 508
89b48e9e 509 for(Int_t ipt=0; ipt<npts; ipt++) {
510 AliTrackPoint point;
511 Int_t modId;
512 array->GetPoint(point,ipt);
513 Int_t volId = point.GetVolumeID();
df69bb07 514 if (volId == kVtxSensVID) continue; // this is a vertex constraint
89b48e9e 515 Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId);
516 if(layerId==3 || layerId==4){
f09774b4 517 drTime[nPtSDD] = point.GetDriftTime();
89b48e9e 518 modId+=AliITSgeomTGeo::GetModuleIndex(layerId,1,1);
519 Int_t index=modId-kNSPDmods;
826883ac 520 if (fDoSDDDriftTime) {
521 fHistSDDDrTimeAll[index]->Fill(drTime[nPtSDD]);
522 if(point.IsExtra()) fHistSDDDrTimeExtra[index]->Fill(drTime[nPtSDD]);
523 else fHistSDDDrTimeAttac[index]->Fill(drTime[nPtSDD]);
524 }
525 if (fDoSDDdEdxCalib) {
526 Float_t dedxLay=dedx[layerId-3];
527 if(dedxLay>1.) fHistSDDdEdxvsDrTime[index]->Fill(drTime[nPtSDD],dedxLay);
528 }
89b48e9e 529 iPtSDD[nPtSDD] = ipt;
530 modIdSDD[nPtSDD] = modId;
f09774b4 531 modSide[nPtSDD] = point.GetClusterType()&BIT(16) ? 0:1;
89b48e9e 532 point.GetXYZ(posGloF);
533 for(Int_t icoor=0;icoor<3;icoor++) posGlo[icoor]=posGloF[icoor];
534 AliITSgeomTGeo::GlobalToLocal(modId,posGlo,posLoc);
535 xLocSDD[nPtSDD]=posLoc[0];
536 zLocSDD[nPtSDD]=posLoc[2];
537 ++nPtSDD;
538 fFitter->SetCovIScale(ipt,1e-4); // scaling for inverted errors of SDD
539 }else{
540 ++nPtSSDSPD;
541 }
542 }
543 if(nPtSDD>0 && nPtSSDSPD>=2){
6a5f5ccc 544 double pt = fUseTPCMomentum ? track->GetTPCInnerParam()->Pt() : track->Pt();
545 fFitter->Fit(track->Charge(),pt,0.);
89b48e9e 546 Double_t chi2=fFitter->GetChi2NDF();
547 if ( chi2<0 || chi2>1e4 ) return; // fit failed, abandon this track
548 for (Int_t ip=0; ip<nPtSDD;ip++) {
549 fFitter->GetResiduals(resGlo,iPtSDD[ip]);
550 TGeoHMatrix *mcurr = AliITSgeomTGeo::GetMatrix(modIdSDD[ip]);
551 mcurr->MasterToLocalVect(resGlo,resLoc);
552 Int_t index=modIdSDD[ip]-kNSPDmods;
826883ac 553 if (fDoSDDResiduals) {
6a5f5ccc 554 fHistSDDResidX[index]->Fill(pt,resLoc[0]);
826883ac 555 fHistSDDResidXvsX[index]->Fill(xLocSDD[ip],resLoc[0]);
556 fHistSDDResidXvsZ[index]->Fill(zLocSDD[ip],resLoc[0]);
557 }
f09774b4 558 //
826883ac 559 if (fDoSDDVDriftCalib) {
560 double cf = modSide[ip] ? 1.e4:-1.e4;
561 double xMeas = cf*xLocSDD[ip]; // measured coordinate in microns
562 double xRes = cf*resLoc[0]; // X residual in microns
563 double xDriftTrue = 3.5085e4 - (xMeas + xRes); // "true" drift distance
564 //
565 fHProfSDDResidXvsXD[index][modSide[ip]]->Fill(xDriftTrue, xRes);
566 fHProfSDDResidXvsZ[index][modSide[ip]]->Fill(zLocSDD[ip]*1e4, xRes);
567 fHProfSDDDrTimevsXD[index][modSide[ip]]->Fill(xDriftTrue, drTime[ip]);
568 fHProfSDDDrTimevsZ[index][modSide[ip]]->Fill(zLocSDD[ip]*1e4, drTime[ip]);
569 }
9ea163ca 570 }
571 }
572}
573//___________________________________________________________________________
574void AliAnalysisTaskITSAlignQA::FitAndFillSDDz(Int_t iLayer, const AliTrackPointArray *array, Int_t npts, AliESDtrack * track){
575 // fit track and fills histos for SDD along z
576
577 fFitter->AttachPoints(array,0, npts-1);
578 Int_t iPtSDD[4],modIdSDD[4];
579 Double_t xLocSDD[4],zLocSDD[4];
580 Int_t nPtSDD=0;
581 Double_t resGlo[3],resLoc[3];
582 Float_t posGloF[3];
583 Double_t posGlo[3],posLoc[3];
584 for(Int_t ipt=0; ipt<npts; ipt++) {
585 AliTrackPoint point;
586 Int_t modId;
587 array->GetPoint(point,ipt);
588 Int_t volId = point.GetVolumeID();
df69bb07 589 if (volId == kVtxSensVID) continue; // this is a vertex constraint
9ea163ca 590 Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId);
591 if(layerId==iLayer){
592 modId+=AliITSgeomTGeo::GetModuleIndex(layerId,1,1);
593 iPtSDD[nPtSDD] = ipt;
594 modIdSDD[nPtSDD] = modId;
595 point.GetXYZ(posGloF);
596 for(Int_t icoor=0;icoor<3;icoor++) posGlo[icoor]=posGloF[icoor];
597 AliITSgeomTGeo::GlobalToLocal(modId,posGlo,posLoc);
598 xLocSDD[nPtSDD]=posLoc[0];
599 zLocSDD[nPtSDD]=posLoc[2];
600 ++nPtSDD;
601 fFitter->SetCovIScale(ipt,1e-4); // scaling for inverted errors of SDD
602 }
603 }
604 if(nPtSDD>0){
6a5f5ccc 605 double pt = fUseTPCMomentum ? track->GetTPCInnerParam()->Pt() : track->Pt();
606 fFitter->Fit(track->Charge(),pt,0.);
9ea163ca 607 Double_t chi2=fFitter->GetChi2NDF();
608 if ( chi2<0 || chi2>1e4 ) return; // fit failed, abandon this track
609 for (Int_t ip=0; ip<nPtSDD;ip++) {
610 fFitter->GetResiduals(resGlo,iPtSDD[ip]);
611 TGeoHMatrix *mcurr = AliITSgeomTGeo::GetMatrix(modIdSDD[ip]);
612 mcurr->MasterToLocalVect(resGlo,resLoc);
613 Int_t index=modIdSDD[ip]-kNSPDmods;
6a5f5ccc 614 fHistSDDResidZ[index]->Fill(pt,resLoc[2]);
89b48e9e 615 fHistSDDResidZvsX[index]->Fill(xLocSDD[ip],resLoc[2]);
616 fHistSDDResidZvsZ[index]->Fill(zLocSDD[ip],resLoc[2]);
617 }
618 }
619}
620//___________________________________________________________________________
621void AliAnalysisTaskITSAlignQA::FitAndFillSSD(Int_t iLayer, const AliTrackPointArray *array, Int_t npts,AliESDtrack * track){
622 // fit track and fills histos for SSD
623 fFitter->AttachPoints(array,0, npts-1);
624 Int_t iPtSSD[4],modIdSSD[4];
625 Int_t nPtSSD=0;
626 Double_t resGlo[3],resLoc[3];
627 for(Int_t ipt=0; ipt<npts; ipt++) {
628 AliTrackPoint point;
629 Int_t modId;
630 array->GetPoint(point,ipt);
631 Int_t volId = point.GetVolumeID();
df69bb07 632 if (volId == kVtxSensVID) continue; // this is a vertex constraint
89b48e9e 633 Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId);
634 if(layerId==iLayer){
635 modId+=AliITSgeomTGeo::GetModuleIndex(layerId,1,1);
636 iPtSSD[nPtSSD] = ipt;
637 modIdSSD[nPtSSD] = modId;
638 ++nPtSSD;
639 fFitter->SetCovIScale(ipt,1e-4);
640 }
641 }
642 if(nPtSSD>0){
6a5f5ccc 643 double pt = fUseTPCMomentum ? track->GetTPCInnerParam()->Pt() : track->Pt();
644 fFitter->Fit(track->Charge(),pt,0.);
89b48e9e 645 Double_t chi2=fFitter->GetChi2NDF();
646 if ( chi2<0 || chi2>1e4 ) return; // fit failed, abandon this track
647 for (Int_t ip=0; ip<nPtSSD;ip++) {
648 fFitter->GetResiduals(resGlo,iPtSSD[ip]);
649 TGeoHMatrix *mcurr = AliITSgeomTGeo::GetMatrix(modIdSSD[ip]);
650 mcurr->MasterToLocalVect(resGlo,resLoc);
651 Int_t index=modIdSSD[ip]-kNSPDmods-kNSDDmods;
6a5f5ccc 652 fHistSSDResidX[index]->Fill(pt,resLoc[0]);
653 fHistSSDResidZ[index]->Fill(pt,resLoc[2]);
89b48e9e 654 }
655 }
656}
657//______________________________________________________________________________
658void AliAnalysisTaskITSAlignQA::Terminate(Option_t */*option*/)
659{
660 // Terminate analysis
661 fOutput = dynamic_cast<TList*> (GetOutputData(1));
662 if (!fOutput) {
663 printf("ERROR: fOutput not available\n");
664 return;
665 }
666
667 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
96e2c861 668 if(fHistNEvents){
00a2ed4c 669 printf("Number of analyzed events = %d, %d tracks accepted\n",
670 (Int_t)fHistNEvents->GetBinContent(kEvAcc+1),(Int_t)fHistNEvents->GetBinContent(kNTracks+1));
96e2c861 671 }else{
672 printf("Warning: pointer to fHistNEvents is NULL\n");
673 }
89b48e9e 674 return;
675}
676
677
678//___________________________________________________________________________
679void AliAnalysisTaskITSAlignQA::LoadGeometryFromOCDB(){
680 //method to get the gGeomanager
681 // it is called at the CreatedOutputObject stage
682 // to comply with the CAF environment
683 AliInfo("Loading geometry");
684
685 AliCDBManager *man = AliCDBManager::Instance();
686 man->SetDefaultStorage(fOCDBLocation.Data());
687 man->SetRun(fRunNb);
688 AliCDBEntry* obj = man->Get(AliCDBPath("GRP", "Geometry", "Data"));
96e2c861 689 if(obj){
690 AliGeomManager::SetGeometry((TGeoManager*)obj->GetObject());
691 AliGeomManager::GetNalignable("ITS");
692 AliGeomManager::ApplyAlignObjsFromCDB("ITS");
693 }
694 else AliFatal("Geometry object not found in OCDB");
89b48e9e 695}
696
697
df69bb07 698//______________________________________________________________________________________
699AliTrackPointArray* AliAnalysisTaskITSAlignQA::PrepareTrack(const AliTrackPointArray* inp, const AliESDVertex* vtx)
700{
701 // Extract from the global TrackPointArray the ITS part and optionally add vertex as the last measured point
702 //
703 int npts = inp->GetNPoints();
704 int modID=0,nptITS = 0;
705 int itsRefs[24];
706 const UShort_t *vids = inp->GetVolumeID();
707 for(int ipt=0; ipt<npts; ipt++) { // count ITS points
708 if (vids[ipt]<=0) continue;
709 int layerId = AliGeomManager::VolUIDToLayer(vids[ipt],modID);
710 if(layerId<1 || layerId>6) continue;
711 itsRefs[nptITS++] = ipt;
712 }
713 //
714 AliTrackPointArray *trackCopy = new AliTrackPointArray(nptITS + (vtx ? 1:0)); // reserve extra space if vertex provided
715 AliTrackPoint point;
716 for(int ipt=0; ipt<nptITS; ipt++) {
717 inp->GetPoint(point,itsRefs[ipt]);
718 trackCopy->AddPoint(ipt,&point);
719 }
720 //
721 if (vtx) {
722 PrepareVertexConstraint(vtx,point);
723 trackCopy->AddPoint(nptITS,&point); // add vertex constraint as a last point
724 }
725 return trackCopy;
726}
89b48e9e 727
df69bb07 728//_______________________________________________________________________________________
729void AliAnalysisTaskITSAlignQA::PrepareVertexConstraint(const AliESDVertex* vtx, AliTrackPoint &point)
730{
731 // convert vertex to measured point with dummy VID
732 if (!vtx) return;
733 //
734 double cmat[6];
735 float cmatF[6];
736 point.SetVolumeID(kVtxSensVID);
737 //
738 vtx->GetCovMatrix(cmat);
739 cmatF[0] = cmat[0]; // xx
740 cmatF[1] = cmat[1]; // xy
741 cmatF[2] = cmat[3]; // xz
742 cmatF[3] = cmat[2]; // yy
743 cmatF[4] = cmat[4]; // yz
744 cmatF[5] = cmat[5]; // zz
745 point.SetXYZ(vtx->GetX(),vtx->GetY(),vtx->GetZ(), cmatF);
746}
89b48e9e 747
5104b23e 748
00a2ed4c 749//_______________________________________________________________________________________
750Bool_t AliAnalysisTaskITSAlignQA::AcceptCentrality(const AliESDEvent *esd) const
751{
752 // check if events is in the required multiplicity range
753 //
754 const AliMultiplicity *alimult = esd->GetMultiplicity();
755 Int_t nclsSPDouter=0;
756 if(alimult) nclsSPDouter = alimult->GetNumberOfITSClusters(1);
757 if(nclsSPDouter<fMinMult || nclsSPDouter>fMaxMult) return kFALSE;
758 //
759 return kTRUE;
760}
6a5f5ccc 761
762//_______________________________________________________________________________________
763void AliAnalysisTaskITSAlignQA::CreateUserInfo()
764{
765 // if needed, set user info of the output tree
766 AliAODHandler* handler = dynamic_cast<AliAODHandler*>( AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler() );
767 if (!handler) return;
768 TTree* outTree = handler->GetTree();
769 //
770 const TMap *cdbMap = AliCDBManager::Instance()->GetStorageMap();
771 const TList *cdbList = AliCDBManager::Instance()->GetRetrievedIds();
772 //
773 TMap *cdbMapCopy = new TMap(cdbMap->GetEntries());
774 cdbMapCopy->SetOwner(1);
775 cdbMapCopy->SetName("cdbMap");
776 TIter iter(cdbMap->GetTable());
777 //
778 TPair* pair = 0;
779 while((pair = dynamic_cast<TPair*> (iter.Next()))){
780 TObjString* keyStr = dynamic_cast<TObjString*> (pair->Key());
781 TObjString* valStr = dynamic_cast<TObjString*> (pair->Value());
782 if (keyStr && valStr) cdbMapCopy->Add(new TObjString(keyStr->GetName()), new TObjString(valStr->GetName()));
783 }
784
785 TList *cdbListCopy = new TList();
786 cdbListCopy->SetOwner(1);
787 cdbListCopy->SetName("cdbList");
788 //
789 TIter iter2(cdbList);
790 AliCDBId* id=0;
791 while((id = dynamic_cast<AliCDBId*> (iter2.Next()))){
792 cdbListCopy->Add(new TObjString(id->ToString().Data()));
793 }
794 //
795 outTree->GetUserInfo()->Add(cdbMapCopy);
796 outTree->GetUserInfo()->Add(cdbListCopy);
797 //
798 AliMagF *fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
799 Double_t bz = fld ? fld->SolenoidField() : 0;
800 TString bzString; bzString+=bz;
801 TObjString *bzObjString = new TObjString(bzString);
802 TList *bzList = new TList();
803 bzList->SetOwner(1);
804 bzList->SetName("BzkGauss");
805 bzList->Add(bzObjString);
806 outTree->GetUserInfo()->Add(bzList);
807 //
808}