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