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