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