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