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