PWGPP/ITS/AliAnalysisTaskITSAlignQA.cxx

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