Updating ITS macros (Ruben)
[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 "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"
15 #include "AliMultiplicity.h"
16 #include <TSystem.h>
17 #include <TTree.h>
18 #include <TH1F.h>
19 #include <TH2F.h>
20 #include <TProfile.h>
21 #include <TChain.h>
22 #include <TGeoGlobalMagField.h>
23 #include "AliAODHandler.h"
24 #include "AliESDInputHandlerRP.h"
25 #include "AliITSSumTP.h"
26 #include "AliMagF.h"
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
54 ClassImp(AliAnalysisTaskITSAlignQA)
55 //______________________________________________________________________________
56 AliAnalysisTaskITSAlignQA::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),
64   fDoSDDVDriftCalib(kTRUE),
65   fDoSDDDriftTime(kTRUE),
66   fDoFillTPTree(kFALSE),
67   fUseITSsaTracks(kFALSE),
68   fLoadGeometry(kFALSE),
69   fUseVertex(kFALSE),
70   fUseVertexForZOnly(kFALSE),
71   fUseTPCMomentum(kFALSE),
72   fMinVtxContributors(5),
73   fRemovePileupWithSPD(kTRUE),
74   fMinITSpts(3),
75   fMinTPCpts(70),
76   fMinPt(0.5),
77   fNPtBins(8),
78   fMinMult(0),
79   fMaxMult(1e9),
80   fCutDCAXY(1.e10),
81   fCutDCAZ(1.e10),
82   fFitter(0),
83   fITSSumTP(),
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 //___________________________________________________________________________
96 AliAnalysisTaskITSAlignQA::~AliAnalysisTaskITSAlignQA(){
97   //
98   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutput;
99   delete fFitter;
100   delete fITSSumTP;
101   //
102 }
103 //___________________________________________________________________________
104 void AliAnalysisTaskITSAlignQA::UserCreateOutputObjects() {
105   //
106
107   if(fLoadGeometry) LoadGeometryFromOCDB();
108
109   fOutput = new TList();
110   fOutput->SetOwner();
111   fOutput->SetName("OutputHistos");
112
113   fHistNEvents = new TH1F("hNEvents", "Number of processed events",kNEvStatBins,-0.5,kNEvStatBins-0.5);
114   fHistNEvents->Sumw2();
115   fHistNEvents->SetMinimum(0);
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");
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();
129   if(fDoSDDResiduals || fDoSDDdEdxCalib || fDoSDDVDriftCalib || fDoSDDDriftTime) CreateSDDHistos();
130   if(fDoSSDResiduals) CreateSSDHistos();
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   //
142   PostData(1,fOutput);
143 }
144
145 //___________________________________________________________________________
146 void 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 //___________________________________________________________________________
169 void AliAnalysisTaskITSAlignQA::CreateSDDHistos(){
170   // Histos for SDD
171
172   for(Int_t iMod=0; iMod<kNSDDmods; iMod++){
173     if (fDoSDDResiduals) {
174       fHistSDDResidX[iMod] = new TH2F(Form("hSDDResidX%d",iMod+kNSPDmods),
175                                       Form("hSDDResidX%d",iMod+kNSPDmods),
176                                       fNPtBins,fPtBinLimits,
177                                       300,-0.15,0.15);
178       fHistSDDResidX[iMod]->Sumw2();
179       fOutput->Add(fHistSDDResidX[iMod]);
180       
181       fHistSDDResidZ[iMod] = new TH2F(Form("hSDDResidZ%d",iMod+kNSPDmods),
182                                       Form("hSDDResidZ%d",iMod+kNSPDmods),
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]);
211       //
212     }
213     //
214     if (fDoSDDVDriftCalib) {
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       }
240     }
241     
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     }
249     //
250     if (fDoSDDDriftTime) {
251       fHistSDDDrTimeAll[iMod]=new TH1F(Form("hSDDDrTimeAll%d",iMod+kNSPDmods),
252                                        Form("hSDDDrTimeAll%d",iMod+kNSPDmods),
253                                        3200,0.,6400.);
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     }
272   }
273   return;
274   //
275 }
276
277 //___________________________________________________________________________
278 void 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 //______________________________________________________________________________
298 void AliAnalysisTaskITSAlignQA::UserExec(Option_t *)
299 {
300   //
301   static AliTrackPointArray* arrayITS = 0;
302   AliTrackPointArray* arrayITSNoVtx = 0;
303   //
304   AliESDEvent *esd = (AliESDEvent*) (InputEvent());
305   if (fITSSumTP) fITSSumTP->Reset();
306
307   if(!esd) {
308     printf("AliAnalysisTaskITSAlignQA::Exec(): bad ESD\n");
309     return;
310   } 
311
312   if(!ESDfriend()) {
313     printf("AliAnalysisTaskITSAlignQA::Exec(): bad ESDfriend\n");
314     return;
315   }
316   //
317   if (!AcceptCentrality(esd)) return;
318   fHistNEvents->Fill(kEvCnt);
319
320   const AliESDVertex* vtx=0,*vtxSPD=0;
321   fHistNEvents->Fill(kEvAll);
322   vtx    = esd->GetPrimaryVertex();
323   vtxSPD = esd->GetPrimaryVertexSPD();
324   //
325   if (fUseVertex) {  // check the vertex if it is requested as an extra point
326     if (!AcceptVertex(vtx,vtxSPD)) return;
327   }
328
329   fHistNEvents->Fill(kEvVtx);
330   if (fRemovePileupWithSPD){
331     // skip events tagged by SPD as pileup
332     if(esd->IsPileupFromSPD()) return;
333   }
334   fHistNEvents->Fill(kEvPlp);
335
336   //
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++) {
343     //
344     if (arrayITS) {delete arrayITS; arrayITS = 0;}  // reset points from previous tracks 
345     arrayITSNoVtx = 0;
346     //
347     AliESDtrack * track = esd->GetTrack(itrack);
348     if(!track) continue;
349     if(!AcceptTrack(track, vtx)) continue;
350     array = track->GetTrackPointArray();
351     if(!array) continue;
352     arrayITS = PrepareTrack(array, vtx);
353     if (fITSSumTP) {
354       arrayITSNoVtx = PrepareTrack(array, 0);
355       arrayITSNoVtx->SetUniqueID(itrack);
356       fITSSumTP->AddTrack(arrayITSNoVtx);
357     }
358     //
359     fHistNEvents->Fill(kNTracks);
360     //
361     Int_t npts  = arrayITS->GetNPoints();
362     Int_t npts1 = fUseVertexForZOnly ? npts-1 : npts;
363     //
364     if(fDoSPDResiduals){ 
365       FitAndFillSPD(1,arrayITS,npts1,track);
366       FitAndFillSPD(2,arrayITS,npts1,track);
367     }
368     if(fDoSDDResiduals || fDoSDDdEdxCalib || fDoSDDVDriftCalib || fDoSDDDriftTime) {
369       FitAndFillSDDrphi(arrayITS,npts,track);
370       if (fDoSDDResiduals) {
371         FitAndFillSDDz(3,arrayITS,npts1,track);
372         FitAndFillSDDz(4,arrayITS,npts1,track);
373       }
374     }
375     if(fDoSSDResiduals){ 
376       FitAndFillSSD(5,arrayITS,npts1,track);
377       FitAndFillSSD(6,arrayITS,npts1,track);
378     }
379   }
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   }
406
407   //
408   PostData(1,fOutput);
409   
410 }
411
412 //___________________________________________________________________________
413 Bool_t AliAnalysisTaskITSAlignQA::AcceptTrack(const AliESDtrack * track, const AliESDVertex* vtx)
414 {
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;
425   Float_t pt = 0;
426   if (fUseTPCMomentum) {
427     if (track->GetTPCInnerParam()) pt = track->GetTPCInnerParam()->Pt();
428     else                           pt = track->Pt();
429   }
430   if(pt<fMinPt) accept=kFALSE;
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   //
443   if(accept) fHistPtAccept->Fill(pt);
444   return accept;
445 }
446
447 //___________________________________________________________________________
448 Bool_t AliAnalysisTaskITSAlignQA::AcceptVertex(const AliESDVertex * vtx, const AliESDVertex * vtxSPD) {
449   // vertex selection cuts
450   if (!vtx || vtx->GetStatus()<1) return kFALSE;
451   if (!vtxSPD || vtxSPD->GetStatus()<1) return kFALSE;
452   if (vtx->GetNContributors()<fMinVtxContributors) return kFALSE;
453   if (TMath::Abs(vtx->GetZ()-vtxSPD->GetZ())>0.3) return kFALSE;
454   return kTRUE;
455 }
456
457 //___________________________________________________________________________
458 void 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();
469     if (volId == kVtxSensVID) continue; // this is a vertex constraint
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){
480     double pt = fUseTPCMomentum ? track->GetTPCInnerParam()->Pt() : track->Pt();
481     fFitter->Fit(track->Charge(),pt,0.);
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];
489       fHistSPDResidX[index]->Fill(pt,resLoc[0]);
490       fHistSPDResidZ[index]->Fill(pt,resLoc[2]);
491     }
492   }    
493 }
494 //___________________________________________________________________________
495 void AliAnalysisTaskITSAlignQA::FitAndFillSDDrphi(const AliTrackPointArray *array, Int_t npts, AliESDtrack * track){
496   // fit track and fills histos for SDD along rphi (drift coord.)
497   Double_t dedx[4];
498   track->GetITSdEdxSamples(dedx);
499
500   fFitter->AttachPoints(array,0, npts-1); 
501   Int_t iPtSDD[4],modIdSDD[4],modSide[4];
502   Double_t xLocSDD[4],zLocSDD[4],drTime[4];
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];
508
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();
514     if (volId == kVtxSensVID) continue; // this is a vertex constraint
515     Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId);
516     if(layerId==3 || layerId==4){
517       drTime[nPtSDD] = point.GetDriftTime();
518       modId+=AliITSgeomTGeo::GetModuleIndex(layerId,1,1);
519       Int_t index=modId-kNSPDmods;
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       }
529       iPtSDD[nPtSDD] = ipt;
530       modIdSDD[nPtSDD] = modId;
531       modSide[nPtSDD] = point.GetClusterType()&BIT(16) ? 0:1; 
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){
544     double pt = fUseTPCMomentum ? track->GetTPCInnerParam()->Pt() : track->Pt();
545     fFitter->Fit(track->Charge(),pt,0.);
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;
553       if (fDoSDDResiduals) {
554         fHistSDDResidX[index]->Fill(pt,resLoc[0]);
555         fHistSDDResidXvsX[index]->Fill(xLocSDD[ip],resLoc[0]);
556         fHistSDDResidXvsZ[index]->Fill(zLocSDD[ip],resLoc[0]);
557       }
558       //
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       }
570     }
571   }
572 }
573 //___________________________________________________________________________
574 void 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();
589     if (volId == kVtxSensVID) continue; // this is a vertex constraint
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){
605     double pt = fUseTPCMomentum ? track->GetTPCInnerParam()->Pt() : track->Pt();
606     fFitter->Fit(track->Charge(),pt,0.);
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;
614       fHistSDDResidZ[index]->Fill(pt,resLoc[2]);
615       fHistSDDResidZvsX[index]->Fill(xLocSDD[ip],resLoc[2]);
616       fHistSDDResidZvsZ[index]->Fill(zLocSDD[ip],resLoc[2]);
617     }
618   }
619 }
620 //___________________________________________________________________________
621 void 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();
632     if (volId == kVtxSensVID) continue; // this is a vertex constraint
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){
643     double pt = fUseTPCMomentum ? track->GetTPCInnerParam()->Pt() : track->Pt();
644     fFitter->Fit(track->Charge(),pt,0.);
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;
652       fHistSSDResidX[index]->Fill(pt,resLoc[0]);
653       fHistSSDResidZ[index]->Fill(pt,resLoc[2]);
654     }
655   }
656 }
657 //______________________________________________________________________________
658 void 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"));
668   if(fHistNEvents){
669     printf("Number of analyzed events = %d, %d tracks accepted\n",
670            (Int_t)fHistNEvents->GetBinContent(kEvAcc+1),(Int_t)fHistNEvents->GetBinContent(kNTracks+1));
671   }else{
672     printf("Warning: pointer to fHistNEvents is NULL\n");
673   }
674   return;
675 }
676
677
678 //___________________________________________________________________________
679 void 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"));
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");
695 }
696
697
698 //______________________________________________________________________________________
699 AliTrackPointArray* 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 }
727
728 //_______________________________________________________________________________________
729 void 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 }
747
748
749 //_______________________________________________________________________________________
750 Bool_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 }
761
762 //_______________________________________________________________________________________
763 void 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 }