PWG1/global/AliAnalysisTaskVertexESD.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 AliAnalysisTaskVertexESD
  18 // AliAnalysisTask to extract from ESD the information for the analysis
  19 // of the primary vertex reconstruction efficiency and resolution
  20 // (for MC events) and distributions (for real data). Three vertices:
  21 // - SPD tracklets
  22 // - ITS+TPC tracks
  23 // - TPC-only tracks
  24 //
  25 // Author: A.Dainese, andrea.dainese@pd.infn.it
  26 //*************************************************************************
  27
  28 #include <TChain.h>
  29 #include <TTree.h>
  30 #include <TNtuple.h>
  31 #include <TBranch.h>
  32 #include <TClonesArray.h>
  33 #include <TObjArray.h>
  34 #include <TH1F.h>
  35 #include <TH2F.h>
  36 #include <TCanvas.h>
  37
  38 #include "AliAnalysisTask.h"
  39 #include "AliAnalysisManager.h"
  40
  41 #include "AliMultiplicity.h"
  42 #include "AliVertexerTracks.h"
  43 #include "AliESDtrack.h"
  44 #include "AliExternalTrackParam.h"
  45 #include "AliESDVertex.h"
  46 #include "AliESDEvent.h"
  47 #include "AliESDInputHandler.h"
  48 #include "AliTrackReference.h"
  49 //#include "AliTriggerAnalysis.h"
  50
  51 #include "AliMCEventHandler.h"
  52 #include "AliMCEvent.h"
  53 #include "AliStack.h"
  54 #include "AliLog.h"
  55
  56 #include "AliGenEventHeader.h"
  57 #include "AliAnalysisTaskVertexESD.h"
  58
  59
  60 ClassImp(AliAnalysisTaskVertexESD)
  61
  62 //________________________________________________________________________
  63 AliAnalysisTaskVertexESD::AliAnalysisTaskVertexESD(const char *name) :
  64 AliAnalysisTask(name, "VertexESDTask"),
  65 fCheckEventType(kTRUE),
  66 fReadMC(kFALSE),
  67 fRecoVtxTPC(kFALSE),
  68 fRecoVtxITSTPC(kFALSE),
  69 fOnlyITSTPCTracks(kFALSE),
  70 fOnlyITSSATracks(kFALSE),
  71 fFillNtuple(kFALSE),
  72 fESD(0),
  73 fOutput(0),
  74 fNtupleVertexESD(0),
  75 fhSPDVertexX(0),
  76 fhSPDVertexY(0),
  77 fhSPDVertexZ(0),
  78 fhTRKVertexX(0),
  79 fhTRKVertexY(0),
  80 fhTRKVertexZ(0),
  81 fhTPCVertexX(0),
  82 fhTPCVertexY(0),
  83 fhTPCVertexZ(0),
  84 fhTrackRefs(0)
  85 {
  86   // Constructor
  87
  88   // Define input and output slots here
  89   // Input slot #0 works with a TChain
  90   DefineInput(0, TChain::Class());
  91   // Output slot #0 writes into a TList container
  92   DefineOutput(0, TList::Class());  //My private output
  93 }
  94 //________________________________________________________________________
  95 AliAnalysisTaskVertexESD::~AliAnalysisTaskVertexESD()
  96 {
  97   // Destructor
  98
  99   // histograms are in the output list and deleted when the output
 100   // list is deleted by the TSelector dtor
 101
 102   if (fOutput) {
 103     delete fOutput;
 104     fOutput = 0;
 105   }
 106 }
 107
 108 //________________________________________________________________________
 109 void AliAnalysisTaskVertexESD::ConnectInputData(Option_t *)
 110 {
 111   // Connect ESD or AOD here
 112   // Called once
 113
 114   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 115   if (!tree) {
 116     Printf("ERROR: Could not read chain from input slot 0");
 117   } else {
 118
 119     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 120
 121     if (!esdH) {
 122       Printf("ERROR: Could not get ESDInputHandler");
 123     } else {
 124       fESD = esdH->GetEvent();
 125     }
 126
 127   }
 128
 129   return;
 130 }
 131
 132 //________________________________________________________________________
 133 void AliAnalysisTaskVertexESD::CreateOutputObjects()
 134 {
 135   // Create histograms
 136   // Called once
 137
 138   // Several histograms are more conveniently managed in a TList
 139   fOutput = new TList;
 140   fOutput->SetOwner();
 141
 142   fNtupleVertexESD = new TNtuple("fNtupleVertexESD","vertices","run:tstamp:xtrue:ytrue:ztrue:xSPD:xerrSPD:ySPD:yerrSPD:zSPD:zerrSPD:ntrksSPD:xTPC:xerrTPC:yTPC:yerrTPC:zTPC:zerrTPC:ntrksTPC:xTRK:xerrTRK:yTRK:yerrTRK:zTRK:zerrTRK:ntrksTRK:ntrklets:nESDtracks:nITSrefit5or6:nTPCin:nTPCinEta09:dndygen:triggered:SPD3D:SPD0cls:constrTRK:constrTPC");
 143
 144   fOutput->Add(fNtupleVertexESD);
 145
 146   fhSPDVertexX = new TH1F("fhSPDVertexX","SPDVertex x; x vertex [cm]; events",200,-1,1);
 147   fOutput->Add(fhSPDVertexX);
 148   fhSPDVertexY = new TH1F("fhSPDVertexY","SPDVertex y; y vertex [cm]; events",200,-1,1);
 149   fOutput->Add(fhSPDVertexY);
 150   fhSPDVertexZ = new TH1F("fhSPDVertexZ","SPDVertex z; z vertex [cm]; events",200,-20,20);
 151   fOutput->Add(fhSPDVertexZ);
 152   fhTRKVertexX = new TH1F("fhTRKVertexX","TRKVertex x; x vertex [cm]; events",200,-1,1);
 153   fOutput->Add(fhTRKVertexX);
 154   fhTRKVertexY = new TH1F("fhTRKVertexY","TRKVertex y; y vertex [cm]; events",200,-1,1);
 155   fOutput->Add(fhTRKVertexY);
 156   fhTRKVertexZ = new TH1F("fhTRKVertexZ","TRKVertex z; z vertex [cm]; events",200,-20,20);
 157   fOutput->Add(fhTRKVertexZ);
 158   fhTPCVertexX = new TH1F("fhTPCVertexX","TPCVertex x; x vertex [cm]; events",200,-3,3);
 159   fOutput->Add(fhTPCVertexX);
 160   fhTPCVertexY = new TH1F("fhTPCVertexY","TPCVertex y; y vertex [cm]; events",200,-3,3);
 161   fOutput->Add(fhTPCVertexY);
 162   fhTPCVertexZ = new TH1F("fhTPCVertexZ","TPCVertex z; z vertex [cm]; events",200,-20,20);
 163   fOutput->Add(fhTPCVertexZ);
 164
 165   fhTrackRefs = new TH2F("fhTrackRefs","Track references; x; y",1000,-4,4,1000,-4,4);
 166   fOutput->Add(fhTrackRefs);
 167
 168   return;
 169 }
 170
 171 //________________________________________________________________________
 172 void AliAnalysisTaskVertexESD::Exec(Option_t *)
 173 {
 174   // Main loop
 175   // Called for each event
 176
 177   if (!fESD) {
 178     Printf("ERROR: fESD not available");
 179     return;
 180   }
 181
 182   if(fCheckEventType && (fESD->GetEventType())!=7) return;
 183
 184
 185   TArrayF mcVertex(3);
 186   mcVertex[0]=9999.; mcVertex[1]=9999.; mcVertex[2]=9999.;
 187   Float_t dNchdy=-999.;
 188
 189   // ***********  MC info ***************
 190   if (fReadMC) {
 191     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 192     if (!eventHandler) {
 193       Printf("ERROR: Could not retrieve MC event handler");
 194       return;
 195     }
 196
 197     AliMCEvent* mcEvent = eventHandler->MCEvent();
 198     if (!mcEvent) {
 199       Printf("ERROR: Could not retrieve MC event");
 200       return;
 201     }
 202
 203     AliStack* stack = mcEvent->Stack();
 204     if (!stack) {
 205       AliDebug(AliLog::kError, "Stack not available");
 206       return;
 207     }
 208
 209     AliHeader* header = mcEvent->Header();
 210     if (!header) {
 211       AliDebug(AliLog::kError, "Header not available");
 212       return;
 213     }
 214     AliGenEventHeader* genHeader = header->GenEventHeader();
 215     genHeader->PrimaryVertex(mcVertex);
 216
 217     Int_t ngenpart = (Int_t)stack->GetNtrack();
 218     //printf("# generated particles = %d\n",ngenpart);
 219     dNchdy=0;
 220     for(Int_t ip=0; ip<ngenpart; ip++) {
 221       TParticle* part = (TParticle*)stack->Particle(ip);
 222       // keep only electorns, muons, pions, kaons and protons
 223       Int_t apdg = TMath::Abs(part->GetPdgCode());
 224       if(apdg!=11 && apdg!=13 && apdg!=211 && apdg!=321 && apdg!=2212) continue;
 225       // reject secondaries
 226       if(TMath::Sqrt((part->Vx()-mcVertex[0])*(part->Vx()-mcVertex[0])+(part->Vy()-mcVertex[1])*(part->Vy()-mcVertex[1]))>0.0010) continue;
 227       // reject incoming protons
 228       Double_t energy  = part->Energy();
 229       if(energy>900.) continue;
 230       Double_t pz = part->Pz();
 231       Double_t y = 0.5*TMath::Log((energy+pz+1.e-13)/(energy-pz+1.e-13));
 232       if(TMath::Abs(y)<1.0) dNchdy += 0.5; // count 1/2 of particles in |y|<1
 233       // tracks refs
 234       TClonesArray *trefs=0;
 235       Int_t ntrefs = mcEvent->GetParticleAndTR(ip,part,trefs);
 236       if(ntrefs<0) continue;
 237       for(Int_t iref=0; iref<ntrefs; iref++) {
 238         AliTrackReference *tref = (AliTrackReference*)trefs->At(iref);
 239         if(tref->R()>10.) continue;
 240         fhTrackRefs->Fill(tref->X(),tref->Y());
 241       }
 242     }
 243     //printf("# primary particles = %7.1f\n",dNchdy);
 244   }
 245   // ***********  MC info ***************
 246
 247
 248   // Trigger
 249   //ULong64_t triggerMask;
 250   //ULong64_t spdFO = (1 << 14);
 251   //ULong64_t v0left = (1 << 10);
 252   //ULong64_t v0right = (1 << 11);
 253
 254   //triggerMask=fESD->GetTriggerMask();
 255   // MB1: SPDFO || V0L || V0R
 256   //Bool_t eventTriggered = (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)));
 257   //MB2: GFO && V0R
 258   //triggerMask & spdFO && ((triggerMask&v0left) || (triggerMask&v0right))
 259   //Bool_t eventTriggered = (triggerMask & spdFO);
 260
 261   //static AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis();
 262   Bool_t eventTriggered = 0;//triggerAnalysis->IsTriggerFired(fESD, AliTriggerAnalysis::kSPDGFO /*| AliTriggerAnalysis::kOfflineFlag*/);
 263
 264   Int_t ntracks = fESD->GetNumberOfTracks();
 265   Int_t nITS5or6=0,nTPCin=0,nTPCinEta09=0;
 266   //printf("Tracks # = %d\n",fESD->GetNumberOfTracks());
 267   for(Int_t itr=0; itr<ntracks; itr++) {
 268     AliESDtrack *t = fESD->GetTrack(itr);
 269     if(t->GetNcls(0)>=5) nITS5or6++;
 270     Double_t z0; t->GetZAt(0,fESD->GetMagneticField(),z0);
 271     if(t->GetNcls(1)>0 && TMath::Abs(t->GetD(0,0,fESD->GetMagneticField()))<2.8 && TMath::Abs(z0)<20) {
 272       nTPCin++;
 273       if(TMath::Abs(t->GetTgl())<1.5) nTPCinEta09++;
 274     }
 275   }
 276
 277
 278   const AliESDVertex *spdv=fESD->GetPrimaryVertexSPD();
 279   const AliESDVertex *tpcv=fESD->GetPrimaryVertexTPC();
 280   const AliESDVertex *trkv=fESD->GetPrimaryVertexTracks();
 281
 282   //Float_t tpccontrorig=tpcv->GetNContributors();
 283
 284   if(fRecoVtxTPC) {
 285     tpcv = 0;
 286     tpcv = ReconstructPrimaryVertexTPC();
 287   }
 288   if(fRecoVtxITSTPC) {
 289     trkv = 0;
 290     trkv = ReconstructPrimaryVertexITSTPC();
 291   }
 292
 293   const AliMultiplicity *alimult = fESD->GetMultiplicity();
 294   Int_t ntrklets=0,spd0cls=0;
 295   if(alimult) {
 296     ntrklets = alimult->GetNumberOfTracklets();
 297     for(Int_t l=0;l<alimult->GetNumberOfTracklets();l++){
 298       if(alimult->GetDeltaPhi(l)<-9998.) ntrklets--;
 299     }
 300     spd0cls = alimult->GetNumberOfSingleClusters()+ntrklets;
 301   }
 302
 303   // fill histos
 304
 305   if(spdv) {
 306     if(spdv->GetNContributors()>0) {
 307       TString title=spdv->GetTitle();
 308       if(title.Contains("3D")) {
 309         fhSPDVertexX->Fill(spdv->GetXv());
 310         fhSPDVertexY->Fill(spdv->GetYv());
 311       }
 312       fhSPDVertexZ->Fill(spdv->GetZv());
 313     }
 314   }
 315
 316   if(trkv) {
 317     if(trkv->GetNContributors()>0) {
 318       fhTRKVertexX->Fill(trkv->GetXv());
 319       fhTRKVertexY->Fill(trkv->GetYv());
 320       fhTRKVertexZ->Fill(trkv->GetZv());
 321     }
 322   }
 323
 324   if(tpcv) {
 325     if(tpcv->GetNContributors()>0) {
 326       fhTPCVertexX->Fill(tpcv->GetXv());
 327       fhTPCVertexY->Fill(tpcv->GetYv());
 328       fhTPCVertexZ->Fill(tpcv->GetZv());
 329     }
 330   }
 331
 332
 333   // fill ntuple
 334   Int_t isize=37;
 335   Float_t xnt[37];
 336
 337   Int_t index=0;
 338
 339   xnt[index++]=(Float_t)fESD->GetRunNumber();
 340   xnt[index++]=(Float_t)fESD->GetTimeStamp();
 341
 342   xnt[index++]=mcVertex[0];
 343   xnt[index++]=mcVertex[1];
 344   xnt[index++]=mcVertex[2];
 345
 346   xnt[index++]=spdv->GetXv();
 347   xnt[index++]=spdv->GetXRes();
 348   xnt[index++]=spdv->GetYv();
 349   xnt[index++]=spdv->GetYRes();
 350   xnt[index++]=spdv->GetZv();
 351   xnt[index++]=spdv->GetZRes();
 352   xnt[index++]=spdv->GetNContributors();
 353
 354   xnt[index++]=tpcv->GetXv();
 355   xnt[index++]=tpcv->GetXRes();
 356   xnt[index++]=tpcv->GetYv();
 357   xnt[index++]=tpcv->GetYRes();
 358   xnt[index++]=tpcv->GetZv();
 359   xnt[index++]=tpcv->GetZRes();
 360   xnt[index++]=tpcv->GetNContributors();
 361
 362   xnt[index++]=trkv->GetXv();
 363   xnt[index++]=trkv->GetXRes();
 364   xnt[index++]=trkv->GetYv();
 365   xnt[index++]=trkv->GetYRes();
 366   xnt[index++]=trkv->GetZv();
 367   xnt[index++]=trkv->GetZRes();
 368   xnt[index++]=trkv->GetNContributors();// tpccontrorig;
 369
 370
 371   xnt[index++]=float(ntrklets);
 372   xnt[index++]=float(ntracks);
 373   xnt[index++]=float(nITS5or6);
 374   xnt[index++]=float(nTPCin);
 375   xnt[index++]=float(nTPCinEta09);
 376
 377   xnt[index++]=float(dNchdy);
 378
 379   xnt[index++]=(eventTriggered ? 1. : 0.);
 380
 381   TString spdtitle = spdv->GetTitle();
 382   xnt[index++]=(spdtitle.Contains("vertexer: 3D") ? 1. : 0.);
 383
 384   xnt[index++]=spd0cls;
 385
 386   TString trktitle = trkv->GetTitle();
 387   xnt[index++]=(trktitle.Contains("WithConstraint") ? 1. : 0.);
 388
 389   TString tpctitle = tpcv->GetTitle();
 390   xnt[index++]=(tpctitle.Contains("WithConstraint") ? 1. : 0.);
 391
 392
 393   if(index!=isize) printf("AliAnalysisTaskVertexESD: ERROR, index!=isize\n");
 394
 395   if(fFillNtuple) fNtupleVertexESD->Fill(xnt);
 396
 397   // Post the data already here
 398   PostData(0, fOutput);
 399
 400   return;
 401 }
 402
 403 //________________________________________________________________________
 404 void AliAnalysisTaskVertexESD::Terminate(Option_t *)
 405 {
 406   // Draw result to the screen
 407   // Called once at the end of the query
 408   fOutput = dynamic_cast<TList*> (GetOutputData(0));
 409   if (!fOutput) {
 410     Printf("ERROR: fOutput not available");
 411     return;
 412   }
 413
 414   fNtupleVertexESD = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleVertexESD"));
 415
 416   return;
 417 }
 418
 419 //_________________________________________________________________________
 420 AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexTPC() const {
 421   // On the fly reco of TPC vertex from ESD
 422
 423   AliVertexerTracks vertexer(fESD->GetMagneticField());
 424   vertexer.SetTPCMode(); // defaults
 425   //vertexer.SetTPCMode(0.1,1.0,5.,0,1,3.,0.1,1.5);
 426   Double_t pos[3]={+0.0220,-0.0340,+0.270};
 427   Double_t err[3]={0.0200,0.0200,7.5};
 428   AliESDVertex *initVertex = new AliESDVertex(pos,err);
 429   vertexer.SetVtxStart(initVertex);
 430   delete initVertex;
 431   vertexer.SetConstraintOff();
 432
 433   return vertexer.FindPrimaryVertex(fESD);
 434 }
 435
 436 //_________________________________________________________________________
 437 AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexITSTPC() const {
 438   // On the fly reco of ITS+TPC vertex from ESD
 439
 440   AliVertexerTracks vertexer(fESD->GetMagneticField());
 441   vertexer.SetITSMode(); // defaults
 442   //vertexer.SetTPCMode(0.1,1.0,5.,0,1,3.,0.1,1.5);
 443   Double_t pos[3]={+0.0220,-0.0340,+0.270};
 444   Double_t err[3]={0.0200,0.0200,7.5};
 445   AliESDVertex *initVertex = new AliESDVertex(pos,err);
 446   vertexer.SetVtxStart(initVertex);
 447   delete initVertex;
 448   vertexer.SetConstraintOff();
 449
 450   // use only ITS-TPC or only ITS-SA tracks
 451   if(fOnlyITSTPCTracks || fOnlyITSSATracks) {
 452     Int_t iskip=0;
 453     Int_t *skip = new Int_t[fESD->GetNumberOfTracks()];
 454     for(Int_t itr=0;itr<fESD->GetNumberOfTracks(); itr++) {
 455       AliESDtrack* track = fESD->GetTrack(itr);
 456       if(fOnlyITSTPCTracks && track->GetNcls(1)==0) { // skip ITSSA
 457         skip[iskip++]=itr;
 458       }
 459       if(fOnlyITSSATracks && track->GetNcls(1)>0) { // skip ITSTPC
 460         skip[iskip++]=itr;
 461       }
 462     }
 463     vertexer.SetSkipTracks(iskip,skip);
 464     delete [] skip; skip=NULL;
 465   }
 466
 467   return vertexer.FindPrimaryVertex(fESD);
 468 }