PWGPP/ITS/AliMeanVertexCalibTask.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //*************************************************************************
  17 // Class AliMeanVertexCalibTask
  18 // AliAnalysisTask to extract from ESD the information on primary vertex
  19 // reconstruction in order to compute the MeanVertex object
  20 //
  21 // Author:  D.Caffarri, davide.caffarri@pd.infn.it
  22 //          A.Dainese, andrea.dainese@pd.infn.it
  23 //*************************************************************************
  24
  25
  26 #include <TH1F.h>
  27 #include <TH2F.h>
  28 #include <string.h>
  29
  30 #include "AliMultiplicity.h"
  31 #include "AliESDVertex.h"
  32 #include "AliESDEvent.h"
  33 #include "AliVertexerTracks.h"
  34 #include "AliCDBManager.h"
  35 #include "AliCDBEntry.h"
  36 #include "AliGRPRecoParam.h"
  37
  38 #include "AliMeanVertexCalibTask.h"
  39
  40
  41 ClassImp(AliMeanVertexCalibTask)
  42
  43 //_____________________________________________________________________
  44 AliMeanVertexCalibTask::AliMeanVertexCalibTask(const char *name):
  45   AliAnalysisTaskSE(name),
  46   fESD(0),
  47   fOutput(0),
  48   fOnlyITSTPCTracks(kFALSE),
  49   fOnlyITSSATracks(kTRUE)
  50 {
  51
  52   // Constructor
  53
  54   // Define input and output slots here
  55   // Output slot #0 writes into a TList container
  56   DefineOutput(1, TList::Class());  //My private output
  57 }
  58
  59 //________________________________________________________________________
  60 AliMeanVertexCalibTask::~AliMeanVertexCalibTask()
  61 {
  62   // Destructor
  63
  64   if (fOutput) {
  65     delete fOutput;
  66     fOutput = 0;
  67   }
  68
  69 }
  70 //________________________________________________________________________
  71 void AliMeanVertexCalibTask::UserCreateOutputObjects()
  72 {
  73
  74   // Several histograms are more conveniently managed in a TList
  75   fOutput = new TList;
  76   fOutput->SetOwner();
  77
  78
  79   TH1F* hSPDVertexX = new TH1F("hSPDVertexX","SPDVertex x; x vertex [cm]; events",200,-1,1);
  80   fOutput->Add(hSPDVertexX);
  81   TH1F* hSPDVertexY = new TH1F("hSPDVertexY","SPDVertex y; y vertex [cm]; events",200,-1,1);
  82   fOutput->Add(hSPDVertexY);
  83   TH1F* hSPDVertexZ = new TH1F("hSPDVertexZ","SPDVertex z; z vertex [cm]; events",200,-20,20);
  84   fOutput->Add(hSPDVertexZ);
  85   TH1F* hTRKVertexX = new TH1F("hTRKVertexX","TRKVertex x; x vertex [cm]; events",200,-1,1);
  86   fOutput->Add(hTRKVertexX);
  87   TH1F* hTRKVertexY = new TH1F("hTRKVertexY","TRKVertex y; y vertex [cm]; events",200,-1,1);
  88   fOutput->Add(hTRKVertexY);
  89   TH1F* hTRKVertexZ = new TH1F("hTRKVertexZ","TRKVertex z; z vertex [cm]; events",200,-20,20);
  90   fOutput->Add(hTRKVertexZ);
  91
  92   TH2F *hTRKVertexXvsMult = new TH2F("hTRKVertexXvsMult", "TRKVertex X vs mult", 200, -1, 1, 300, 0, 3000);
  93   fOutput->Add(hTRKVertexXvsMult);
  94
  95   TH2F *hTRKVertexYvsMult = new TH2F("hTRKVertexYvsMult", "TRKVertex Y vs mult", 200, -1, 1, 300, 0, 3000);
  96   fOutput->Add(hTRKVertexYvsMult);
  97
  98   TH2F *hTRKVertexXZ = new TH2F("hTRKVertexXZ", "TRKVertex XZ corr", 200, -1, 1, 200, -20, 20);
  99   fOutput->Add(hTRKVertexXZ);
 100
 101   TH2F *hTRKVertexYZ = new TH2F("hTRKVertexYZ", "TRKVertex YZ corr", 200, -1, 1, 200, -20, 20);
 102   fOutput->Add(hTRKVertexYZ);
 103
 104   TH1F* hTRKVertexXdefMult = new TH1F("hTRKVertexXdefMult","TRKVertex x Mult; x vertex [cm] 30<Mult<45; events",500,-1,1);
 105   fOutput->Add(hTRKVertexXdefMult);
 106   TH1F* hTRKVertexYdefMult = new TH1F("hTRKVertexYdefMult","TRKVertex y Mult; y vertex [cm] 30<Mult<45; events",500,-1,1);
 107   fOutput->Add(hTRKVertexYdefMult);
 108
 109   TH1F* hTRKVertexXHighMult = new TH1F("hTRKVertexXHighMult","TRKVertex x High Mult; x vertex [cm] Mult>1500; events",500,-0.5,0.5);
 110   fOutput->Add(hTRKVertexXHighMult);
 111   TH1F* hTRKVertexYHighMult = new TH1F("hTRKVertexYHighMult","TRKVertex y High Mult; y vertex [cm] Mult>1500; events",500,-0.5,0.5);
 112   fOutput->Add(hTRKVertexYHighMult);
 113
 114   TH1F* hITSSAVertexX = new TH1F("hITSSAVertexX","ITSSAVertex x; x vertex [cm]; events",200,-1,1);
 115   fOutput->Add(hITSSAVertexX);
 116   TH1F* hITSSAVertexY = new TH1F("hITSSAVertexY","ITSSAVertex y; y vertex [cm]; events",200,-1,1);
 117   fOutput->Add(hITSSAVertexY);
 118   TH1F* hITSSAVertexZ = new TH1F("hITSSAVertexZ","ITSSAVertex z; z vertex [cm]; events",200,-20,20);
 119   fOutput->Add(hITSSAVertexZ);
 120
 121   TH2F *hITSSAVertexXZ = new TH2F("hITSSAVertexXZ", "ITSSAVertex XZ corr", 200, -1, 1, 200, -20, 20);
 122   fOutput->Add(hITSSAVertexXZ);
 123
 124   TH2F *hITSSAVertexYZ = new TH2F("hITSSAVertexYZ", "ITSSAVertex YZ corr", 200, -1, 1, 200, -20, 20);
 125   fOutput->Add(hITSSAVertexYZ);
 126
 127   TH2F *hITSSAVertexXvsMult = new TH2F("hITSSAVertexXvsMult", "ITSSAVertex X vs mult", 200, -1, 1, 300, 0, 3000);
 128   fOutput->Add(hITSSAVertexXvsMult);
 129
 130   TH2F *hITSSAVertexYvsMult = new TH2F("hITSSAVertexYvsMult", "ITSSAVertex Y vs mult", 200, -1, 1, 300, 0, 3000);
 131   fOutput->Add(hITSSAVertexYvsMult);
 132
 133   TH1F* hITSSAVertexXdefMult = new TH1F("hITSSAVertexXdefMult","ITSSAVertex x Mult; x vertex [cm] 30<Mult<45; events",500,-1,1);
 134   fOutput->Add(hITSSAVertexXdefMult);
 135   TH1F* hITSSAVertexYdefMult = new TH1F("hITSSAVertexYdefMult","ITSSAVertex y Mult; y vertex [cm] 30<Mult<45; events",500,-1,1);
 136   fOutput->Add(hITSSAVertexYdefMult);
 137
 138
 139   TH1F* hITSSAVertexXHighMult = new TH1F("hITSSAVertexXHighMult","ITSSAVertex x High Mult; x vertex [cm] Mult>1500; events",500,-0.5,0.5);
 140   fOutput->Add(hITSSAVertexXHighMult);
 141   TH1F* hITSSAVertexYHighMult = new TH1F("hITSSAVertexYHighMult","ITSSAVertex y High Mult; y vertex [cm] Mult>1500; events",500,-0.5,0.5);
 142   fOutput->Add(hITSSAVertexYHighMult);
 143
 144   PostData(1, fOutput);
 145
 146   return;
 147 }
 148
 149
 150 //________________________________________________________________________
 151 void AliMeanVertexCalibTask::UserExec(Option_t *)
 152 {
 153   // Main loop
 154   // Called for each event
 155
 156   if (!InputEvent()) {
 157     Printf("ERROR: fESD not available");
 158     return;
 159   }
 160
 161   AliESDEvent* esdE = (AliESDEvent*) InputEvent();
 162
 163   const AliMultiplicity *alimult = esdE->GetMultiplicity();
 164   Int_t ntrklets=0;
 165   if(alimult) ntrklets = alimult->GetNumberOfTracklets();
 166
 167   const char* beamType = esdE->GetBeamType();
 168   // Printf("beam type = %s", beamType);
 169
 170   Bool_t kLowFlux = kTRUE, kHighFlux = kFALSE;
 171   // TString pp= "p-p";
 172   //TString pA= "p-A";
 173   TString AA= "A-A";
 174
 175   if (beamType == AA){
 176     kHighFlux = kTRUE;
 177     kLowFlux = kFALSE;
 178     // Printf ("high flux setting");
 179     }
 180
 181   AliCDBManager* man = AliCDBManager::Instance();
 182   //man->SetDefaultStorage("raw://");
 183   Int_t runNb = esdE->GetRunNumber();
 184   if (runNb > 0) {
 185     man->SetRun(runNb);
 186     // Printf("runNb = %d", runNb);
 187   }
 188
 189   AliCDBEntry *entry = (AliCDBEntry*)man->Get("GRP/Calib/RecoParam/");
 190   // Printf("entry = %p", entry);
 191   TObjArray *arrayRecoParam=0x0;
 192   if (entry) {
 193     arrayRecoParam = (TObjArray*)entry->GetObject();
 194     // Printf("arrayRecoParam = %p", arrayRecoParam);
 195   }
 196   else {
 197     Printf("CDBEntry not found");
 198     return;
 199   }
 200   AliGRPRecoParam *grpRecoParam=0x0;
 201   if (kLowFlux) grpRecoParam= (AliGRPRecoParam*)arrayRecoParam->At(1);
 202   else if (kHighFlux) grpRecoParam = (AliGRPRecoParam*)arrayRecoParam->At(2);
 203
 204   AliVertexerTracks *vertexer= new AliVertexerTracks(esdE->GetMagneticField());
 205   vertexer->SetITSMode();
 206   vertexer->SetConstraintOff();
 207
 208   if (grpRecoParam) {
 209     Int_t nCutsVertexer = grpRecoParam->GetVertexerTracksNCuts();
 210     Double_t *cutsVertexer = new Double_t[nCutsVertexer];
 211     grpRecoParam->GetVertexerTracksCutsITS(cutsVertexer,nCutsVertexer);
 212     vertexer->SetCuts(cutsVertexer,nCutsVertexer);
 213     delete [] cutsVertexer; cutsVertexer = NULL;
 214    }
 215
 216   vertexer->SetConstraintOff();
 217
 218   //track vertex recomputed from the vertexer
 219   AliESDVertex *trkv = vertexer->FindPrimaryVertex(esdE);
 220
 221   //const AliESDVertex *trkv = esdE->GetPrimaryVertexTracks();
 222
 223   //SPD vertex taken from the ESD
 224   const AliESDVertex *spdv = esdE->GetPrimaryVertexSPD();
 225
 226   //ITSSA vertex computed from the vertexer
 227   vertexer->SetITSpureSA();
 228   AliESDVertex *itsSAv = vertexer->FindPrimaryVertex(esdE);
 229
 230   if(spdv) {
 231     if(spdv->GetNContributors()>0) {
 232       TString title=spdv->GetTitle();
 233       if(title.Contains("3D")) {
 234         ((TH1F*)fOutput->FindObject("hSPDVertexX"))->Fill(spdv->GetX());
 235         ((TH1F*)fOutput->FindObject("hSPDVertexY"))->Fill(spdv->GetY());
 236       }
 237       ((TH1F*)fOutput->FindObject("hSPDVertexZ"))->Fill(spdv->GetZ());
 238     }
 239   }
 240
 241
 242   if(trkv) {
 243     if(trkv->GetNContributors()>0) {
 244       ((TH1F*)fOutput->FindObject("hTRKVertexX"))->Fill(trkv->GetX());
 245       ((TH1F*)fOutput->FindObject("hTRKVertexY"))->Fill(trkv->GetY());
 246       ((TH1F*)fOutput->FindObject("hTRKVertexZ"))->Fill(trkv->GetZ());
 247
 248       ((TH2F*)fOutput->FindObject("hTRKVertexXvsMult"))->Fill(trkv->GetX(), ntrklets);
 249       ((TH2F*)fOutput->FindObject("hTRKVertexYvsMult"))->Fill(trkv->GetY(), ntrklets);
 250
 251       if (ntrklets>30 && ntrklets<45){
 252         ((TH1F*)fOutput->FindObject("hTRKVertexXdefMult"))->Fill(trkv->GetX());
 253         ((TH1F*)fOutput->FindObject("hTRKVertexYdefMult"))->Fill(trkv->GetY());
 254       }
 255
 256       if (ntrklets>1500){
 257         ((TH1F*)fOutput->FindObject("hTRKVertexXHighMult"))->Fill(trkv->GetX());
 258         ((TH1F*)fOutput->FindObject("hTRKVertexYHighMult"))->Fill(trkv->GetY());
 259       }
 260
 261       ((TH2F*)fOutput->FindObject("hTRKVertexXZ"))->Fill(trkv->GetX(),trkv->GetZ());
 262       ((TH2F*)fOutput->FindObject("hTRKVertexYZ"))->Fill(trkv->GetY(),trkv->GetZ());
 263
 264     }
 265   }
 266
 267   if (itsSAv){
 268     if (itsSAv->GetNContributors()>0){
 269
 270       ((TH1F*)fOutput->FindObject("hITSSAVertexX"))->Fill(itsSAv->GetX());
 271       ((TH1F*)fOutput->FindObject("hITSSAVertexY"))->Fill(itsSAv->GetY());
 272       ((TH1F*)fOutput->FindObject("hITSSAVertexZ"))->Fill(itsSAv->GetZ());
 273
 274       ((TH2F*)fOutput->FindObject("hITSSAVertexXvsMult"))->Fill(itsSAv->GetX(), ntrklets);
 275       ((TH2F*)fOutput->FindObject("hITSSAVertexYvsMult"))->Fill(itsSAv->GetY(), ntrklets);
 276
 277       if (ntrklets>30 && ntrklets<45){
 278         ((TH1F*)fOutput->FindObject("hITSSAVertexXdefMult"))->Fill(itsSAv->GetX());
 279         ((TH1F*)fOutput->FindObject("hITSSAVertexYdefMult"))->Fill(itsSAv->GetY());
 280       }
 281
 282       if (ntrklets>1500){
 283         ((TH1F*)fOutput->FindObject("hITSSAVertexXHighMult"))->Fill(itsSAv->GetX());
 284         ((TH1F*)fOutput->FindObject("hITSSAVertexYHighMult"))->Fill(itsSAv->GetY());
 285       }
 286
 287       ((TH2F*)fOutput->FindObject("hITSSAVertexXZ"))->Fill(itsSAv->GetX(),itsSAv->GetZ());
 288       ((TH2F*)fOutput->FindObject("hITSSAVertexYZ"))->Fill(itsSAv->GetY(),itsSAv->GetZ());
 289
 290     }
 291   }
 292
 293   delete itsSAv;
 294   delete vertexer;
 295
 296   PostData(1, fOutput);
 297
 298   return;
 299
 300 }
 301
 302
 303 //________________________________________________________________________
 304 void AliMeanVertexCalibTask::Terminate(Option_t *)
 305 {
 306   // Draw result to the screen
 307   // Called once at the end of the query
 308   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 309   if (!fOutput) {
 310     Printf("ERROR: fOutput not available");
 311     return;
 312   }
 313
 314
 315   return;
 316
 317 }
 318
 319
 320 //__________________________________________________________________________
 321 // AliESDVertex* AliMeanVertexCalibTask::ReconstructPrimaryVertex(Bool_t constr,Int_t mode) const {
 322 //   // On the fly reco of ITS+TPC vertex from ESD
 323 //   // mode=0 use all tracks
 324 //   // mode=1 use odd-number tracks
 325 //   // mode=2 use even-number tracks
 326
 327 //   AliESDEvent* evt = (AliESDEvent*) fInputEvent;
 328 //   AliVertexerTracks vertexer(evt->GetMagneticField());
 329 //   if(evt->GetNumberOfTracks()<500) {
 330 //     vertexer.SetITSMode(); // defaults
 331 //     vertexer.SetMinClusters(4); // default is 5
 332 //   } else {
 333 //     vertexer.SetITSMode(0.1,0.1,0.5,5,1,3.,100.,1000.,3.,30.,1,1);// PbPb
 334 //   }
 335 //   if (fOnlyITSSATracks) vertexer.SetITSpureSA();
 336
 337 //   Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy);
 338 //   Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0};
 339 //   Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
 340 //   AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1);
 341 //   vertexer.SetVtxStart(initVertex);
 342 //   delete initVertex;
 343 //   if(!constr) vertexer.SetConstraintOff();
 344
 345 //   if(fOnlyITSTPCTracks || fOnlyITSSATracks || mode>0) {
 346 //     Int_t iskip=0;
 347 //     Int_t ntracks = evt->GetNumberOfTracks();
 348 //     Int_t *skip = new Int_t[ntracks];
 349 //     for(Int_t i=0;i<ntracks;i++) skip[i]=-1;
 350 //     for(Int_t itr=0;itr<ntracks; itr++) {
 351 //       AliESDtrack* track = evt->GetTrack(itr);
 352 //       if(fOnlyITSTPCTracks && track->GetNcls(1)==0) { // skip ITSSA
 353 //      skip[iskip++]=itr;
 354 //      continue;
 355 //       }
 356 //       if(fOnlyITSSATracks && track->GetNcls(1)>0) { // skip ITSTPC
 357 //      skip[iskip++]=itr;
 358 //      continue;
 359 //       }
 360 //       if(mode==1 && itr%2==0) skip[iskip++]=itr;
 361 //       if(mode==2 && itr%2==1) skip[iskip++]=itr;
 362 //     }
 363 //     vertexer.SetSkipTracks(iskip,skip);
 364 //     delete [] skip; skip=NULL;
 365 //   }
 366
 367 //   return vertexer.FindPrimaryVertex(evt);
 368 // }