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 fESD(0),
  66 fOutput(0),
  67 fNtupleVertexESD(0),
  68 fhTrackRefs(0)
  69 {
  70   // Constructor
  71
  72   // Define input and output slots here
  73   // Input slot #0 works with a TChain
  74   DefineInput(0, TChain::Class());
  75   // Output slot #0 writes into a TList container
  76   DefineOutput(0, TList::Class());  //My private output
  77 }
  78 //________________________________________________________________________
  79 AliAnalysisTaskVertexESD::~AliAnalysisTaskVertexESD()
  80 {
  81   // Destructor
  82
  83   // histograms are in the output list and deleted when the output
  84   // list is deleted by the TSelector dtor
  85
  86   if (fOutput) {
  87     delete fOutput;
  88     fOutput = 0;
  89   }
  90 }
  91
  92 //________________________________________________________________________
  93 void AliAnalysisTaskVertexESD::ConnectInputData(Option_t *)
  94 {
  95   // Connect ESD or AOD here
  96   // Called once
  97
  98   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  99   if (!tree) {
 100     Printf("ERROR: Could not read chain from input slot 0");
 101   } else {
 102
 103     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 104
 105     if (!esdH) {
 106       Printf("ERROR: Could not get ESDInputHandler");
 107     } else {
 108       fESD = esdH->GetEvent();
 109     }
 110
 111   }
 112
 113   return;
 114 }
 115
 116 //________________________________________________________________________
 117 void AliAnalysisTaskVertexESD::CreateOutputObjects()
 118 {
 119   // Create histograms
 120   // Called once
 121
 122   // Several histograms are more conveniently managed in a TList
 123   fOutput = new TList;
 124   fOutput->SetOwner();
 125
 126   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");
 127
 128   fOutput->Add(fNtupleVertexESD);
 129
 130   fhTrackRefs = new TH2F("fhTrackRefs","Track references; x; y",1000,-4,4,1000,-4,4);
 131   fOutput->Add(fhTrackRefs);
 132
 133  return;
 134 }
 135
 136 //________________________________________________________________________
 137 void AliAnalysisTaskVertexESD::Exec(Option_t *)
 138 {
 139   // Main loop
 140   // Called for each event
 141
 142   if (!fESD) {
 143     Printf("ERROR: fESD not available");
 144     return;
 145   }
 146
 147
 148   TArrayF mcVertex(3);
 149   mcVertex[0]=9999.; mcVertex[1]=9999.; mcVertex[2]=9999.;
 150   Float_t dNchdy=-999.;
 151
 152   // ***********  MC info ***************
 153   if (fReadMC) {
 154     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 155     if (!eventHandler) {
 156       Printf("ERROR: Could not retrieve MC event handler");
 157       return;
 158     }
 159
 160     AliMCEvent* mcEvent = eventHandler->MCEvent();
 161     if (!mcEvent) {
 162       Printf("ERROR: Could not retrieve MC event");
 163       return;
 164     }
 165
 166     AliStack* stack = mcEvent->Stack();
 167     if (!stack) {
 168       AliDebug(AliLog::kError, "Stack not available");
 169       return;
 170     }
 171
 172     AliHeader* header = mcEvent->Header();
 173     if (!header) {
 174       AliDebug(AliLog::kError, "Header not available");
 175       return;
 176     }
 177     AliGenEventHeader* genHeader = header->GenEventHeader();
 178     genHeader->PrimaryVertex(mcVertex);
 179
 180     Int_t ngenpart = (Int_t)stack->GetNtrack();
 181     //printf("# generated particles = %d\n",ngenpart);
 182     dNchdy=0;
 183     for(Int_t ip=0; ip<ngenpart; ip++) {
 184       TParticle* part = (TParticle*)stack->Particle(ip);
 185       // keep only electorns, muons, pions, kaons and protons
 186       Int_t apdg = TMath::Abs(part->GetPdgCode());
 187       if(apdg!=11 && apdg!=13 && apdg!=211 && apdg!=321 && apdg!=2212) continue;
 188       // reject secondaries
 189       if(TMath::Sqrt((part->Vx()-mcVertex[0])*(part->Vx()-mcVertex[0])+(part->Vy()-mcVertex[1])*(part->Vy()-mcVertex[1]))>0.0010) continue;
 190       // reject incoming protons
 191       Double_t energy  = part->Energy();
 192       if(energy>900.) continue;
 193       Double_t pz = part->Pz();
 194       Double_t y = 0.5*TMath::Log((energy+pz+1.e-13)/(energy-pz+1.e-13));
 195       if(TMath::Abs(y)<1.0) dNchdy += 0.5; // count 1/2 of particles in |y|<1
 196       // tracks refs
 197       TClonesArray *trefs=0;
 198       Int_t ntrefs = mcEvent->GetParticleAndTR(ip,part,trefs);
 199       if(ntrefs<0) continue;
 200       for(Int_t iref=0; iref<ntrefs; iref++) {
 201         AliTrackReference *tref = (AliTrackReference*)trefs->At(iref);
 202         if(tref->R()>10.) continue;
 203         fhTrackRefs->Fill(tref->X(),tref->Y());
 204       }
 205     }
 206     //printf("# primary particles = %7.1f\n",dNchdy);
 207   }
 208   // ***********  MC info ***************
 209
 210
 211   // Trigger
 212   ULong64_t triggerMask;
 213   ULong64_t spdFO = (1 << 14);
 214   ULong64_t v0left = (1 << 11);
 215   ULong64_t v0right = (1 << 12);
 216
 217   triggerMask=fESD->GetTriggerMask();
 218   // MB1: SPDFO || V0L || V0R
 219   Bool_t eventTriggered = (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)));
 220   //MB2: GFO && V0R
 221   //triggerMask & spdFO && ((triggerMask&v0left) || (triggerMask&v0right))
 222
 223
 224   Int_t ntracks = fESD->GetNumberOfTracks();
 225   Int_t nITS5or6=0,nTPCin=0,nTPCinEta09=0;
 226   //printf("Tracks # = %d\n",fESD->GetNumberOfTracks());
 227   for(Int_t itr=0; itr<ntracks; itr++) {
 228     AliESDtrack *t = fESD->GetTrack(itr);
 229     if(t->GetNcls(0)>=5) nITS5or6++;
 230     Double_t z0; t->GetZAt(0,fESD->GetMagneticField(),z0);
 231     if(t->GetNcls(1)>0 && TMath::Abs(t->GetD(0,0,fESD->GetMagneticField()))<2.8 && TMath::Abs(z0)<20) {
 232       nTPCin++;
 233       if(TMath::Abs(t->GetTgl())<1.5) nTPCinEta09++;
 234     }
 235   }
 236
 237
 238   const AliESDVertex *spdv=fESD->GetPrimaryVertexSPD();
 239   const AliESDVertex *tpcv=fESD->GetPrimaryVertexTPC();
 240   const AliESDVertex *trkv=fESD->GetPrimaryVertex();
 241
 242   //Float_t tpccontrorig=tpcv->GetNContributors();
 243
 244   //tpcv = 0;
 245   //tpcv = ReconstructPrimaryVertexTPC();
 246
 247   const AliMultiplicity *alimult = fESD->GetMultiplicity();
 248   Int_t ntrklets=0,spd0cls=0;
 249   if(alimult) {
 250     ntrklets = alimult->GetNumberOfTracklets();
 251     for(Int_t l=0;l<alimult->GetNumberOfTracklets();l++){
 252       if(alimult->GetDeltaPhi(l)<-9998.) ntrklets--;
 253     }
 254     spd0cls = alimult->GetNumberOfSingleClusters()+ntrklets;
 255   }
 256
 257
 258   Int_t isize=37;
 259   Float_t xnt[37];
 260
 261   Int_t index=0;
 262
 263   xnt[index++]=(Float_t)fESD->GetRunNumber();
 264   xnt[index++]=(Float_t)fESD->GetTimeStamp();
 265
 266   xnt[index++]=mcVertex[0];
 267   xnt[index++]=mcVertex[1];
 268   xnt[index++]=mcVertex[2];
 269
 270   xnt[index++]=spdv->GetXv();
 271   xnt[index++]=spdv->GetXRes();
 272   xnt[index++]=spdv->GetYv();
 273   xnt[index++]=spdv->GetYRes();
 274   xnt[index++]=spdv->GetZv();
 275   xnt[index++]=spdv->GetZRes();
 276   xnt[index++]=spdv->GetNContributors();
 277
 278   xnt[index++]=tpcv->GetXv();
 279   xnt[index++]=tpcv->GetXRes();
 280   xnt[index++]=tpcv->GetYv();
 281   xnt[index++]=tpcv->GetYRes();
 282   xnt[index++]=tpcv->GetZv();
 283   xnt[index++]=tpcv->GetZRes();
 284   xnt[index++]=tpcv->GetNContributors();
 285
 286   xnt[index++]=trkv->GetXv();
 287   xnt[index++]=trkv->GetXRes();
 288   xnt[index++]=trkv->GetYv();
 289   xnt[index++]=trkv->GetYRes();
 290   xnt[index++]=trkv->GetZv();
 291   xnt[index++]=trkv->GetZRes();
 292   xnt[index++]=trkv->GetNContributors();// tpccontrorig;
 293
 294
 295   xnt[index++]=float(ntrklets);
 296   xnt[index++]=float(ntracks);
 297   xnt[index++]=float(nITS5or6);
 298   xnt[index++]=float(nTPCin);
 299   xnt[index++]=float(nTPCinEta09);
 300
 301   xnt[index++]=float(dNchdy);
 302
 303   xnt[index++]=(eventTriggered ? 1. : 0.);
 304
 305   TString spdtitle = spdv->GetTitle();
 306   xnt[index++]=(spdtitle.Contains("vertexer: 3D") ? 1. : 0.);
 307
 308   xnt[index++]=spd0cls;
 309
 310   TString trktitle = trkv->GetTitle();
 311   xnt[index++]=(trktitle.Contains("WithConstraint") ? 1. : 0.);
 312
 313   TString tpctitle = tpcv->GetTitle();
 314   xnt[index++]=(tpctitle.Contains("WithConstraint") ? 1. : 0.);
 315
 316
 317   if(index!=isize) printf("AliAnalysisTaskVertexESD: ERROR, index!=isize\n");
 318
 319   fNtupleVertexESD->Fill(xnt);
 320
 321   // Post the data already here
 322   PostData(0, fOutput);
 323
 324   return;
 325 }
 326
 327 //________________________________________________________________________
 328 void AliAnalysisTaskVertexESD::Terminate(Option_t *)
 329 {
 330   // Draw result to the screen
 331   // Called once at the end of the query
 332   fOutput = dynamic_cast<TList*> (GetOutputData(0));
 333   if (!fOutput) {
 334     Printf("ERROR: fOutput not available");
 335     return;
 336   }
 337
 338   fNtupleVertexESD = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleVertexESD"));
 339
 340   return;
 341 }
 342
 343 //_________________________________________________________________________
 344 AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexTPC() const {
 345   // On the fly reco of TPC vertex from ESD
 346
 347   AliVertexerTracks vertexer(fESD->GetMagneticField());
 348   vertexer.SetTPCMode(); // defaults
 349   //vertexer.SetTPCMode(0.1,1.0,5.,0,1,3.,0.1,1.5);
 350   Double_t pos[3]={+0.0220,-0.0340,+0.270};
 351   Double_t err[3]={0.0200,0.0200,7.5};
 352   AliESDVertex *initVertex = new AliESDVertex(pos,err);
 353   vertexer.SetVtxStart(initVertex);
 354   delete initVertex;
 355   vertexer.SetConstraintOff();
 356
 357   return vertexer.FindPrimaryVertex(fESD);
 358 }