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(100.),
  82   fCutDCAZ(100.),
  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   fHistNEvents->Fill(kEvAll);
 343   //
 344   if (!AcceptCentrality(esd)) return;
 345   fHistNEvents->Fill(kEvCnt);
 346
 347   const AliESDVertex* vtx=0,*vtxSPD=0;
 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     CopyUserInfo();
 437   }
 438
 439   //
 440   PostData(1,fOutput);
 441
 442 }
 443
 444 //___________________________________________________________________________
 445 Bool_t AliAnalysisTaskITSAlignQA::AcceptTrack(const AliESDtrack * track, const AliESDVertex* vtx)
 446 {
 447   // track selection cuts
 448   Bool_t accept=kTRUE;
 449   if(fUseITSsaTracks){
 450     if(track->GetNcls(1)>0) accept=kFALSE;
 451   }else{
 452     if(track->GetNcls(1)<fMinTPCpts) accept=kFALSE;
 453   }
 454   if(track->GetNcls(0) < fMinITSpts) accept=kFALSE;
 455   Int_t trstatus=track->GetStatus();
 456   if(!(trstatus&AliESDtrack::kITSrefit)) accept=kFALSE;
 457   Float_t pt = 0;
 458   if (fUseTPCMomentum && track->IsOn(AliESDtrack::kTPCrefit)) pt = track->GetTPCInnerParam()->Pt();
 459   else pt = track->Pt();
 460   //
 461   if(pt<fMinPt) accept=kFALSE;
 462   //
 463   // if vertex constraint is used, apply soft DCA cut
 464   if (vtx) {
 465     Double_t dz[2],cov[3];
 466     AliExternalTrackParam trc = *track;
 467     if (!trc.PropagateToDCA(vtx, fFitter->GetBz(), 3.0, dz, cov)) accept=kFALSE;
 468     else {
 469       if (dz[0]*dz[0]/(1e-4+cov[0])>fCutDCAXY) accept=kFALSE;
 470       if (dz[1]*dz[1]/(4e-4+cov[2])>fCutDCAZ)  accept=kFALSE;
 471     }
 472   }
 473   //
 474   if(accept) fHistPtAccept->Fill(pt);
 475   return accept;
 476 }
 477
 478 //___________________________________________________________________________
 479 Bool_t AliAnalysisTaskITSAlignQA::AcceptVertex(const AliESDVertex * vtx, const AliESDVertex * vtxSPD) {
 480   // vertex selection cuts
 481   if (!vtx || vtx->GetStatus()<1) return kFALSE;
 482   if (!vtxSPD || vtxSPD->GetStatus()<1) return kFALSE;
 483   if (vtx->GetNContributors()<fMinVtxContributors) return kFALSE;
 484   if (TMath::Abs(vtx->GetZ()-vtxSPD->GetZ())>0.3) return kFALSE;
 485   return kTRUE;
 486 }
 487
 488 //___________________________________________________________________________
 489 void AliAnalysisTaskITSAlignQA::FitAndFillSPD(Int_t iLayer, const AliTrackPointArray *array, Int_t npts,AliESDtrack * track){
 490   // fit track and fills histos for SPD
 491   fFitter->AttachPoints(array,0, npts-1);
 492   Int_t iPtSPD[4],modIdSPD[4];
 493   Int_t nPtSPD=0;
 494   Double_t resGlo[3],resLoc[3];
 495   for(Int_t ipt=0; ipt<npts; ipt++) {
 496     AliTrackPoint point;
 497     Int_t modId;
 498     array->GetPoint(point,ipt);
 499     Int_t volId = point.GetVolumeID();
 500     if (volId == kVtxSensVID) continue; // this is a vertex constraint
 501     Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 502     if(layerId==iLayer){
 503       modId+=AliITSgeomTGeo::GetModuleIndex(layerId,1,1);
 504       iPtSPD[nPtSPD] = ipt;
 505       modIdSPD[nPtSPD] = modId;
 506       ++nPtSPD;
 507       fFitter->SetCovIScale(ipt,1e-4);
 508     }
 509   }
 510   if(nPtSPD>0){
 511     double pt = (fUseTPCMomentum && track->IsOn(AliESDtrack::kTPCrefit)) ? track->GetTPCInnerParam()->Pt() : track->Pt();
 512     fFitter->Fit(track->Charge(),pt,0.);
 513     Double_t chi2=fFitter->GetChi2NDF();
 514     if ( chi2<0 || chi2>1e4 ) return; // fit failed, abandon this track
 515     for (Int_t ip=0; ip<nPtSPD;ip++) {
 516       fFitter->GetResiduals(resGlo,iPtSPD[ip]);
 517       TGeoHMatrix *mcurr = AliITSgeomTGeo::GetMatrix(modIdSPD[ip]);
 518       mcurr->MasterToLocalVect(resGlo,resLoc);
 519       Int_t index=modIdSPD[ip];
 520       fHistSPDResidX[index]->Fill(pt,resLoc[0]);
 521       fHistSPDResidZ[index]->Fill(pt,resLoc[2]);
 522     }
 523   }
 524 }
 525 //___________________________________________________________________________
 526 void AliAnalysisTaskITSAlignQA::FitAndFillSDDrphi(const AliTrackPointArray *array, Int_t npts, AliESDtrack * track){
 527   // fit track and fills histos for SDD along rphi (drift coord.)
 528   Double_t dedx[4];
 529   track->GetITSdEdxSamples(dedx);
 530
 531   fFitter->AttachPoints(array,0, npts-1);
 532   Int_t iPtSDD[4],modIdSDD[4],modSide[4];
 533   Double_t xLocSDD[4],zLocSDD[4],drTime[4];
 534   Int_t nPtSDD=0;
 535   Int_t nPtSSDSPD=0;
 536   Double_t resGlo[3],resLoc[3];
 537   Float_t  posGloF[3];
 538   Double_t posGlo[3],posLoc[3];
 539
 540   for(Int_t ipt=0; ipt<npts; ipt++) {
 541     AliTrackPoint point;
 542     Int_t modId;
 543     array->GetPoint(point,ipt);
 544     Int_t volId = point.GetVolumeID();
 545     if (volId == kVtxSensVID) continue; // this is a vertex constraint
 546     Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 547     if(layerId==3 || layerId==4){
 548       drTime[nPtSDD] = point.GetDriftTime();
 549       modId+=AliITSgeomTGeo::GetModuleIndex(layerId,1,1);
 550       Int_t index=modId-kNSPDmods;
 551       if (fDoSDDDriftTime) {
 552         fHistSDDDrTimeAll[index]->Fill(drTime[nPtSDD]);
 553         if(point.IsExtra()) fHistSDDDrTimeExtra[index]->Fill(drTime[nPtSDD]);
 554         else fHistSDDDrTimeAttac[index]->Fill(drTime[nPtSDD]);
 555       }
 556       if (fDoSDDdEdxCalib) {
 557         Float_t dedxLay=dedx[layerId-3];
 558         if(dedxLay>1.) fHistSDDdEdxvsDrTime[index]->Fill(drTime[nPtSDD],dedxLay);
 559       }
 560       iPtSDD[nPtSDD] = ipt;
 561       modIdSDD[nPtSDD] = modId;
 562       modSide[nPtSDD] = point.GetClusterType()&BIT(16) ? 0:1;
 563       point.GetXYZ(posGloF);
 564       for(Int_t icoor=0;icoor<3;icoor++) posGlo[icoor]=posGloF[icoor];
 565       AliITSgeomTGeo::GlobalToLocal(modId,posGlo,posLoc);
 566       xLocSDD[nPtSDD]=posLoc[0];
 567       zLocSDD[nPtSDD]=posLoc[2];
 568       ++nPtSDD;
 569       fFitter->SetCovIScale(ipt,1e-4); // scaling for inverted errors of SDD
 570     }else{
 571       ++nPtSSDSPD;
 572     }
 573   }
 574   if(nPtSDD>0 && nPtSSDSPD>=2){
 575     double pt = (fUseTPCMomentum && track->IsOn(AliESDtrack::kTPCrefit)) ? track->GetTPCInnerParam()->Pt() : track->Pt();
 576     fFitter->Fit(track->Charge(),pt,0.);
 577     Double_t chi2=fFitter->GetChi2NDF();
 578     if ( chi2<0 || chi2>1e4 ) return; // fit failed, abandon this track
 579     for (Int_t ip=0; ip<nPtSDD;ip++) {
 580       fFitter->GetResiduals(resGlo,iPtSDD[ip]);
 581       TGeoHMatrix *mcurr = AliITSgeomTGeo::GetMatrix(modIdSDD[ip]);
 582       mcurr->MasterToLocalVect(resGlo,resLoc);
 583       Int_t index=modIdSDD[ip]-kNSPDmods;
 584       if (fDoSDDResiduals) {
 585         fHistSDDResidX[index]->Fill(pt,resLoc[0]);
 586         fHistSDDResidXvsX[index]->Fill(xLocSDD[ip],resLoc[0]);
 587         fHistSDDResidXvsZ[index]->Fill(zLocSDD[ip],resLoc[0]);
 588       }
 589       //
 590       if (fDoSDDVDriftCalib) {
 591         double cf = modSide[ip] ? 1.e4:-1.e4;
 592         double xMeas = cf*xLocSDD[ip];            // measured coordinate in microns
 593         double xRes  = cf*resLoc[0];             // X residual in microns
 594         double xDriftTrue  = 3.5085e4 - (xMeas + xRes);   // "true" drift distance
 595         //
 596         fHProfSDDResidXvsXD[index][modSide[ip]]->Fill(xDriftTrue, xRes);
 597         fHProfSDDResidXvsZ[index][modSide[ip]]->Fill(zLocSDD[ip]*1e4, xRes);
 598         fHProfSDDDrTimevsXD[index][modSide[ip]]->Fill(xDriftTrue, drTime[ip]);
 599         fHProfSDDDrTimevsZ[index][modSide[ip]]->Fill(zLocSDD[ip]*1e4, drTime[ip]);
 600       }
 601     }
 602   }
 603 }
 604 //___________________________________________________________________________
 605 void AliAnalysisTaskITSAlignQA::FitAndFillSDDz(Int_t iLayer, const AliTrackPointArray *array, Int_t npts, AliESDtrack * track){
 606   // fit track and fills histos for SDD along z
 607
 608   fFitter->AttachPoints(array,0, npts-1);
 609   Int_t iPtSDD[4],modIdSDD[4];
 610   Double_t xLocSDD[4],zLocSDD[4];
 611   Int_t nPtSDD=0;
 612   Double_t resGlo[3],resLoc[3];
 613   Float_t  posGloF[3];
 614   Double_t posGlo[3],posLoc[3];
 615   for(Int_t ipt=0; ipt<npts; ipt++) {
 616     AliTrackPoint point;
 617     Int_t modId;
 618     array->GetPoint(point,ipt);
 619     Int_t volId = point.GetVolumeID();
 620     if (volId == kVtxSensVID) continue; // this is a vertex constraint
 621     Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 622     if(layerId==iLayer){
 623       modId+=AliITSgeomTGeo::GetModuleIndex(layerId,1,1);
 624       iPtSDD[nPtSDD] = ipt;
 625       modIdSDD[nPtSDD] = modId;
 626       point.GetXYZ(posGloF);
 627       for(Int_t icoor=0;icoor<3;icoor++) posGlo[icoor]=posGloF[icoor];
 628       AliITSgeomTGeo::GlobalToLocal(modId,posGlo,posLoc);
 629       xLocSDD[nPtSDD]=posLoc[0];
 630       zLocSDD[nPtSDD]=posLoc[2];
 631       ++nPtSDD;
 632       fFitter->SetCovIScale(ipt,1e-4); // scaling for inverted errors of SDD
 633     }
 634   }
 635   if(nPtSDD>0){
 636     double pt = (fUseTPCMomentum && track->IsOn(AliESDtrack::kTPCrefit)) ? track->GetTPCInnerParam()->Pt() : track->Pt();
 637     fFitter->Fit(track->Charge(),pt,0.);
 638     Double_t chi2=fFitter->GetChi2NDF();
 639     if ( chi2<0 || chi2>1e4 ) return; // fit failed, abandon this track
 640     for (Int_t ip=0; ip<nPtSDD;ip++) {
 641       fFitter->GetResiduals(resGlo,iPtSDD[ip]);
 642       TGeoHMatrix *mcurr = AliITSgeomTGeo::GetMatrix(modIdSDD[ip]);
 643       mcurr->MasterToLocalVect(resGlo,resLoc);
 644       Int_t index=modIdSDD[ip]-kNSPDmods;
 645       fHistSDDResidZ[index]->Fill(pt,resLoc[2]);
 646       fHistSDDResidZvsX[index]->Fill(xLocSDD[ip],resLoc[2]);
 647       fHistSDDResidZvsZ[index]->Fill(zLocSDD[ip],resLoc[2]);
 648     }
 649   }
 650 }
 651 //___________________________________________________________________________
 652 void AliAnalysisTaskITSAlignQA::FitAndFillSSD(Int_t iLayer, const AliTrackPointArray *array, Int_t npts,AliESDtrack * track){
 653   // fit track and fills histos for SSD
 654   fFitter->AttachPoints(array,0, npts-1);
 655   Int_t iPtSSD[4],modIdSSD[4];
 656   Int_t nPtSSD=0;
 657   Double_t resGlo[3],resLoc[3];
 658   for(Int_t ipt=0; ipt<npts; ipt++) {
 659     AliTrackPoint point;
 660     Int_t modId;
 661     array->GetPoint(point,ipt);
 662     Int_t volId = point.GetVolumeID();
 663     if (volId == kVtxSensVID) continue; // this is a vertex constraint
 664     Int_t layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 665     if(layerId==iLayer){
 666       modId+=AliITSgeomTGeo::GetModuleIndex(layerId,1,1);
 667       iPtSSD[nPtSSD] = ipt;
 668       modIdSSD[nPtSSD] = modId;
 669       ++nPtSSD;
 670       fFitter->SetCovIScale(ipt,1e-4);
 671     }
 672   }
 673   if(nPtSSD>0){
 674     double pt = (fUseTPCMomentum && track->IsOn(AliESDtrack::kTPCrefit)) ? track->GetTPCInnerParam()->Pt() : track->Pt();
 675     fFitter->Fit(track->Charge(),pt,0.);
 676     Double_t chi2=fFitter->GetChi2NDF();
 677     if ( chi2<0 || chi2>1e4 ) return; // fit failed, abandon this track
 678     for (Int_t ip=0; ip<nPtSSD;ip++) {
 679       fFitter->GetResiduals(resGlo,iPtSSD[ip]);
 680       TGeoHMatrix *mcurr = AliITSgeomTGeo::GetMatrix(modIdSSD[ip]);
 681       mcurr->MasterToLocalVect(resGlo,resLoc);
 682       Int_t index=modIdSSD[ip]-kNSPDmods-kNSDDmods;
 683       fHistSSDResidX[index]->Fill(pt,resLoc[0]);
 684       fHistSSDResidZ[index]->Fill(pt,resLoc[2]);
 685     }
 686   }
 687 }
 688 //______________________________________________________________________________
 689 void AliAnalysisTaskITSAlignQA::Terminate(Option_t */*option*/)
 690 {
 691   // Terminate analysis
 692   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 693   if (!fOutput) {
 694     printf("ERROR: fOutput not available\n");
 695     return;
 696   }
 697
 698   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
 699   if(fHistNEvents){
 700     AliInfo(Form("Number of analyzed events = %d, %d tracks accepted",
 701                  (Int_t)fHistNEvents->GetBinContent(kEvAcc+1),(Int_t)fHistNEvents->GetBinContent(kNTracks+1)));
 702   }else{
 703     printf("Warning: pointer to fHistNEvents is NULL\n");
 704   }
 705   //
 706   if (fDoFillTPTree) CreateUserInfo();
 707   //
 708   return;
 709 }
 710
 711
 712 //___________________________________________________________________________
 713 void AliAnalysisTaskITSAlignQA::LoadGeometryFromOCDB(){
 714   //method to get the gGeomanager
 715   // it is called at the CreatedOutputObject stage
 716   // to comply with the CAF environment
 717   AliInfo("Loading geometry");
 718
 719   AliCDBManager *man = AliCDBManager::Instance();
 720   man->SetDefaultStorage(fOCDBLocation.Data());
 721   man->SetRun(fRunNb);
 722   AliCDBEntry* obj = man->Get(AliCDBPath("GRP", "Geometry", "Data"));
 723   if(obj){
 724     AliGeomManager::SetGeometry((TGeoManager*)obj->GetObject());
 725     AliGeomManager::GetNalignable("ITS");
 726     AliGeomManager::ApplyAlignObjsFromCDB("ITS");
 727   }
 728   else AliFatal("Geometry object not found in OCDB");
 729 }
 730
 731
 732 //______________________________________________________________________________________
 733 AliTrackPointArray* AliAnalysisTaskITSAlignQA::PrepareTrack(const AliTrackPointArray* inp, const AliESDVertex* vtx)
 734 {
 735   // Extract from the global TrackPointArray the ITS part and optionally add vertex as the last measured point
 736   //
 737   int npts = inp->GetNPoints();
 738   int modID=0,nptITS = 0;
 739   int itsRefs[24];
 740   const UShort_t *vids = inp->GetVolumeID();
 741   for(int ipt=0; ipt<npts; ipt++) { // count ITS points
 742     if (vids[ipt]<=0) continue;
 743     int layerId = AliGeomManager::VolUIDToLayer(vids[ipt],modID);
 744     if(layerId<1 || layerId>6) continue;
 745     itsRefs[nptITS++] = ipt;
 746   }
 747   //
 748   AliTrackPointArray *trackCopy = new AliTrackPointArray(nptITS + (vtx ? 1:0)); // reserve extra space if vertex provided
 749   AliTrackPoint point;
 750   for(int ipt=0; ipt<nptITS; ipt++) {
 751     inp->GetPoint(point,itsRefs[ipt]);
 752     trackCopy->AddPoint(ipt,&point);
 753   }
 754   //
 755   if (vtx) {
 756     PrepareVertexConstraint(vtx,point);
 757     trackCopy->AddPoint(nptITS,&point); // add vertex constraint as a last point
 758   }
 759   return trackCopy;
 760 }
 761
 762 //_______________________________________________________________________________________
 763 void AliAnalysisTaskITSAlignQA::PrepareVertexConstraint(const AliESDVertex* vtx, AliTrackPoint &point)
 764 {
 765   // convert vertex to measured point with dummy VID
 766   if (!vtx) return;
 767   //
 768   double cmat[6];
 769   float cmatF[6];
 770   point.SetVolumeID(kVtxSensVID);
 771   //
 772   vtx->GetCovMatrix(cmat);
 773   cmatF[0] = cmat[0]; // xx
 774   cmatF[1] = cmat[1]; // xy
 775   cmatF[2] = cmat[3]; // xz
 776   cmatF[3] = cmat[2]; // yy
 777   cmatF[4] = cmat[4]; // yz
 778   cmatF[5] = cmat[5]; // zz
 779   point.SetXYZ(vtx->GetX(),vtx->GetY(),vtx->GetZ(), cmatF);
 780 }
 781
 782
 783 //_______________________________________________________________________________________
 784 Bool_t AliAnalysisTaskITSAlignQA::AcceptCentrality(const AliESDEvent *esd) const
 785 {
 786   // check if events is in the required multiplicity range
 787   //
 788   const AliMultiplicity *alimult = esd->GetMultiplicity();
 789   Int_t nclsSPDouter=0;
 790   if(alimult) nclsSPDouter = alimult->GetNumberOfITSClusters(1);
 791   if(nclsSPDouter<fMinMult || nclsSPDouter>fMaxMult) return kFALSE;
 792   //
 793   return kTRUE;
 794 }
 795
 796 //_______________________________________________________________________________________
 797 void AliAnalysisTaskITSAlignQA::CreateUserInfo()
 798 {
 799   // if needed, set user info of the output tree
 800   if (!fTPTree) {
 801     AliError("TrackPoints summary tree does not exist");
 802     return;
 803   }
 804   TList* uInfo = fTPTree->GetUserInfo();
 805   TMap  *cdbMapCopy  = (TMap*)uInfo->FindObject("cdbMap");
 806   TList *cdbListCopy = (TList*)uInfo->FindObject("cdbList");
 807   TList *bzList      = (TList*)uInfo->FindObject("BzkGauss");
 808   if (cdbMapCopy && cdbListCopy && bzList) return; //already done
 809   //
 810   const TMap *cdbMap = AliCDBManager::Instance()->GetStorageMap();
 811   const TList *cdbList = AliCDBManager::Instance()->GetRetrievedIds();
 812   //
 813   cdbMapCopy = new TMap(cdbMap->GetEntries());
 814   cdbMapCopy->SetOwner(1);
 815   cdbMapCopy->SetName("cdbMap");
 816   TIter iter(cdbMap->GetTable());
 817   //
 818   TPair* pair = 0;
 819   while((pair = dynamic_cast<TPair*> (iter.Next()))){
 820     TObjString* keyStr = dynamic_cast<TObjString*> (pair->Key());
 821     TObjString* valStr = dynamic_cast<TObjString*> (pair->Value());
 822     if (keyStr && valStr) {
 823       cdbMapCopy->Add(new TObjString(keyStr->GetName()), new TObjString(valStr->GetName()));
 824       AliInfo(Form("Add %s : %s to cdbMap of ITSTPUserInfo",keyStr->GetName(),valStr->GetName()));
 825     }
 826   }
 827   //
 828   cdbListCopy = new TList();
 829   cdbListCopy->SetOwner(1);
 830   cdbListCopy->SetName("cdbList");
 831   //
 832   TIter iter2(cdbList);
 833   AliCDBId* id=0;
 834   while((id = dynamic_cast<AliCDBId*> (iter2.Next()))){
 835     cdbListCopy->Add(new TObjString(id->ToString().Data()));
 836     AliInfo(Form("Add %s to cdbList of ITSTPUserInfo",id->ToString().Data()));
 837   }
 838   //
 839   uInfo->Add(cdbMapCopy);
 840   uInfo->Add(cdbListCopy);
 841   //
 842   AliMagF *fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 843   Double_t bz = fld ? fld->SolenoidField() : 0;
 844   TString bzString; bzString+=bz;
 845   TObjString *bzObjString = new TObjString(bzString);
 846   bzList = new TList();
 847   bzList->SetOwner(1);
 848   bzList->SetName("BzkGauss");
 849   bzList->Add(bzObjString);
 850   uInfo->Add(bzList);
 851   //
 852 }
 853
 854 //_______________________________________________________________________________________
 855 void AliAnalysisTaskITSAlignQA::CopyUserInfo()
 856 {
 857   // if available, copy the UserInfo from the ESDtree to the output tree
 858   static Bool_t done = kFALSE;
 859   if (done) return;
 860   if (!fTPTree) {
 861     AliError("TrackPoints summary tree does not exist");
 862     return;
 863   }
 864   AliESDInputHandler *handler = (AliESDInputHandler*)AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 865   TTree* esdTree = 0;
 866   if (!handler || !(esdTree=handler->GetTree())) return;
 867   if (esdTree->InheritsFrom(TChain::Class())) esdTree = esdTree->GetTree();
 868   TList* uInfoSrc = esdTree->GetUserInfo();
 869   const TMap  *cdbMapSrc  = (TMap*)uInfoSrc->FindObject("cdbMap");
 870   const TList *cdbListSrc = (TList*)uInfoSrc->FindObject("cdbList");
 871   if (!cdbMapSrc || !cdbListSrc) return;
 872   //
 873   AliInfo("Create ITSTPUserInfo from esdTree");
 874   TList* uInfoDst = fTPTree->GetUserInfo();
 875   //
 876   TMap *cdbMapCopy = new TMap(cdbMapSrc->GetEntries());
 877   cdbMapCopy->SetOwner(1);
 878   cdbMapCopy->SetName("cdbMap");
 879   TIter iter(cdbMapSrc->GetTable());
 880   TPair* pair = 0;
 881   while((pair = dynamic_cast<TPair*> (iter.Next()))){
 882     TObjString* keyStr = dynamic_cast<TObjString*> (pair->Key());
 883     TObjString* valStr = dynamic_cast<TObjString*> (pair->Value());
 884     if (keyStr && valStr) {
 885       cdbMapCopy->Add(new TObjString(keyStr->GetName()), new TObjString(valStr->GetName()));
 886       AliInfo(Form("Add %s : %s to cdbMap of ITSTPUserInfo",keyStr->GetName(),valStr->GetName()));
 887     }
 888   }
 889   //
 890   TList *cdbListCopy = new TList();
 891   cdbListCopy->SetOwner(1);
 892   cdbListCopy->SetName("cdbList");
 893   //
 894   TIter iter2(cdbListSrc);
 895   TObjString* id=0;
 896   while((id = dynamic_cast<TObjString*> (iter2.Next()))){
 897     cdbListCopy->Add(new TObjString(*id));
 898     AliInfo(Form("Add %s to cdbList of ITSTPUserInfo",id->GetName()));
 899   }
 900   //
 901   uInfoDst->Add(cdbMapCopy);
 902   uInfoDst->Add(cdbListCopy);
 903   //
 904   AliMagF *fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 905   Double_t bz = fld ? fld->SolenoidField() : 0;
 906   TString bzString; bzString+=bz;
 907   TObjString *bzObjString = new TObjString(bzString);
 908   TList *bzList = new TList();
 909   bzList->SetOwner(1);
 910   bzList->SetName("BzkGauss");
 911   bzList->Add(bzObjString);
 912   uInfoDst->Add(bzList);
 913   //
 914   done = kTRUE;
 915 }