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 }