PWGPP/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 #include <TGraphAsymmErrors.h>
  38 #include <TFile.h>
  39
  40 #include "AliAnalysisTask.h"
  41 #include "AliAnalysisManager.h"
  42
  43 #include "AliMultiplicity.h"
  44 #include "AliVertexerTracks.h"
  45 #include "AliESDtrack.h"
  46 #include "AliExternalTrackParam.h"
  47 #include "AliESDVertex.h"
  48 #include "AliVEvent.h"
  49 #include "AliESDInputHandler.h"
  50 #include "AliTrackReference.h"
  51 //#include "AliTriggerAnalysis.h"
  52
  53 #include "AliMCEventHandler.h"
  54 #include "AliMCEvent.h"
  55 #include "AliStack.h"
  56 #include "AliLog.h"
  57
  58 #include "AliGenEventHeader.h"
  59 #include "AliAnalysisTaskVertexESD.h"
  60
  61
  62 ClassImp(AliAnalysisTaskVertexESD)
  63
  64 //________________________________________________________________________
  65 AliAnalysisTaskVertexESD::AliAnalysisTaskVertexESD(const char *name) :
  66 AliAnalysisTaskSE(name),
  67   fCheckEventType(kTRUE),
  68   fReadMC(kFALSE),
  69   fRecoVtxTPC(kTRUE),
  70   fRecoVtxITSTPC(kTRUE),
  71   fRecoVtxITSTPCHalfEvent(kFALSE),
  72   fOnlyITSTPCTracks(kFALSE),
  73   fOnlyITSSATracks(kFALSE),
  74   fFillNtuple(kFALSE),
  75   fFillTreeBeamSpot(kFALSE),
  76   fESD(0),
  77   fOutput(0),
  78   fNtupleVertexESD(0),
  79   fhSPDVertexX(0),
  80   fhSPDVertexY(0),
  81   fhSPDVertexZ(0),
  82   fhSPDVertexZonly(0),
  83   fhTRKVertexX(0),
  84   fhTRKVertexY(0),
  85   fhTRKVertexZ(0),
  86   fhTPCVertexX(0),
  87   fhTPCVertexY(0),
  88   fhTPCVertexZ(0),
  89   fhTrackRefs(0),
  90   fTreeBeamSpot(0),
  91   fhTriggeredTrklets(0),
  92   fhSPD3DTrklets(0),
  93   fhSPDZTrklets(0),
  94   fhTRKTrklets(0),
  95   fhTRKcTrklets(0),
  96   fhTRKncTrklets(0),
  97   fhTPCTrklets(0),
  98   fhTPCcTrklets(0),
  99   fhTPCncTrklets(0),
 100   fhSPD3DZreco(0),
 101   fhSPDZZreco(0),
 102   fhSPDVertexXPile(0),
 103   fhSPDVertexYPile(0),
 104   fhSPDVertexZPile(0),
 105   fhSPDVertexDiffZPileContr2(0),
 106   fhSPDVertexDiffZPileContr3(0),
 107   fhSPDVertexDiffZPileContr4(0),
 108   fhSPDVertexDiffZPileContr5(0),
 109   fhSPDVertexDiffZPileDefault(0),
 110   fhSPDContributorsPile(0),
 111   fhSPDDispContributors(0),
 112   fTriggerType(AliVEvent::kMB),
 113   fhntrksSPDvsSPDcls(0),
 114   fhntrksZvsSPDcls(0)
 115 {
 116   // Constructor
 117
 118   // Define input and output slots here
 119   // Output slot #0 writes into a TList container
 120   DefineOutput(1, TList::Class());  //My private output
 121 }
 122 //________________________________________________________________________
 123 AliAnalysisTaskVertexESD::~AliAnalysisTaskVertexESD()
 124 {
 125   // Destructor
 126
 127   // histograms are in the output list and deleted when the output
 128   // list is deleted by the TSelector dtor
 129
 130   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 131     delete fOutput;
 132     fOutput = 0;
 133   }
 134 }
 135
 136
 137 //________________________________________________________________________
 138 void AliAnalysisTaskVertexESD::UserCreateOutputObjects()
 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:bunchcross:triggered:dndygen:xdiam:ydiam:zdiam:xerrdiam:yerrdiam:zerrdiam:xtrue:ytrue:ztrue:xSPD:xerrSPD:ySPD:yerrSPD:zSPD:zerrSPD:ntrksSPD:SPD3D:dphiSPD:xTPC:xerrTPC:yTPC:yerrTPC:zTPC:zerrTPC:ntrksTPC:constrTPC:xTRK:xerrTRK:yTRK:yerrTRK:zTRK:zerrTRK:ntrksTRK:constrTRK:ntrklets:nESDtracks:nITSrefit5or6:nTPCin:nTPCinEta09:SPD0cls:ispileupfromSPDdefault:xSPDp:xerrSPDp:ySPDp:yerrSPDp:zSPDp:zerrSPDp:ntrksSPDp:xTPCnc:xerrTPCnc:yTPCnc:yerrTPCnc:zTPCnc:zerrTPCnc:ntrksTPCnc:xTPCc:xerrTPCc:yTPCc:yerrTPCc:zTPCc:zerrTPCc:ntrksTPCc:xTRKnc:xerrTRKnc:yTRKnc:yerrTRKnc:zTRKnc:zerrTRKnc:ntrksTRKnc:xTRKc:xerrTRKc:yTRKc:yerrTRKc:zTRKc:zerrTRKc:ntrksTRKc:deltaxTRKnc:deltayTRKnc:deltazTRKnc:deltaxerrTRKnc:deltayerrTRKnc:deltazerrTRKnc:ntrksEvenTRKnc:ntrksOddTRKnc");
 148
 149   fOutput->Add(fNtupleVertexESD);
 150
 151   fhSPDVertexX = new TH1F("fhSPDVertexX","SPDVertex x; x vertex [cm]; events",200,-1,1);
 152   fOutput->Add(fhSPDVertexX);
 153   fhSPDVertexY = new TH1F("fhSPDVertexY","SPDVertex y; y vertex [cm]; events",200,-1,1);
 154   fOutput->Add(fhSPDVertexY);
 155   fhSPDVertexZ = new TH1F("fhSPDVertexZ","SPDVertex z; z vertex [cm]; events",250,-25,25);
 156   fOutput->Add(fhSPDVertexZ);
 157
 158   fhSPDVertexZonly = new TH1F("fhSPDVertexZonly","SPDVertexer z; z vertex [cm]; events",250,-25,25);
 159   fOutput->Add(fhSPDVertexZonly);
 160
 161   fhTRKVertexX = new TH1F("fhTRKVertexX","TRKVertex x; x vertex [cm]; events",2000,-1,1);
 162   fOutput->Add(fhTRKVertexX);
 163   fhTRKVertexY = new TH1F("fhTRKVertexY","TRKVertex y; y vertex [cm]; events",2000,-1,1);
 164   fOutput->Add(fhTRKVertexY);
 165   fhTRKVertexZ = new TH1F("fhTRKVertexZ","TRKVertex z; z vertex [cm]; events",250,-25,25);
 166   fOutput->Add(fhTRKVertexZ);
 167   fhTPCVertexX = new TH1F("fhTPCVertexX","TPCVertex x; x vertex [cm]; events",200,-3,3);
 168   fOutput->Add(fhTPCVertexX);
 169   fhTPCVertexY = new TH1F("fhTPCVertexY","TPCVertex y; y vertex [cm]; events",200,-3,3);
 170   fOutput->Add(fhTPCVertexY);
 171   fhTPCVertexZ = new TH1F("fhTPCVertexZ","TPCVertex z; z vertex [cm]; events",250,-25,25);
 172   fOutput->Add(fhTPCVertexZ);
 173
 174
 175   fhSPDVertexXPile = new TH1F("fhSPDVertexXPile","SPDVertexPile x; x vertex [cm]; events",200,-20,20);
 176   fOutput->Add(fhSPDVertexXPile);
 177   fhSPDVertexYPile = new TH1F("fhSPDVertexYPile","SPDVertexPile y; y vertex [cm]; events",200,-20,20);
 178   fOutput->Add(fhSPDVertexYPile);
 179   fhSPDVertexZPile = new TH1F("fhSPDVertexZPile","SPDVertexPile z; z vertex [cm]; events",200,-40,40);
 180   fOutput->Add(fhSPDVertexZPile);
 181
 182   fhSPDVertexDiffZPileContr2 = new TH1F("fhSPDVertexDiffZPileContr2","SPDVertexDiff z; zmain - zpile [cm]; events",2000,-100,100);
 183   fOutput->Add(fhSPDVertexDiffZPileContr2);
 184   fhSPDVertexDiffZPileContr3 = new TH1F("fhSPDVertexDiffZPileContr3","SPDVertexDiff z; zmain - zpile [cm]; events",2000,-100,100);
 185   fOutput->Add(fhSPDVertexDiffZPileContr3);
 186   fhSPDVertexDiffZPileContr4 = new TH1F("fhSPDVertexDiffZPileContr4","SPDVertexDiff z; zmain - zpile [cm]; events",2000,-100,100);
 187   fOutput->Add(fhSPDVertexDiffZPileContr4);
 188   fhSPDVertexDiffZPileContr5 = new TH1F("fhSPDVertexDiffZPileContr5","SPDVertexDiff z; zmain - zpile [cm]; events",2000,-100,100);
 189   fOutput->Add(fhSPDVertexDiffZPileContr5);
 190   fhSPDVertexDiffZPileDefault = new TH1F("fhSPDVertexDiffZPileDefault","SPDVertexDiff z; zmain - zpile [cm]; events",2000,-100,100);
 191   fOutput->Add(fhSPDVertexDiffZPileDefault);
 192
 193   fhTrackRefs = new TH2F("fhTrackRefs","Track references; x; y",1000,-4,4,1000,-4,4);
 194   fOutput->Add(fhTrackRefs);
 195
 196   fTreeBeamSpot = new TTree("fTreeBeamSpot", "BeamSpotTree");
 197   UShort_t triggered, ntrkletsS, ntrksTRKnc;
 198   UInt_t run, bx;
 199   Float_t cetTimeLHC,xTRKnc, yTRKnc, zTRKnc;
 200   fTreeBeamSpot->Branch("run", &run, "run/i");
 201   fTreeBeamSpot->Branch("cetTimeLHC", &cetTimeLHC, "cetTimeLHC/F");
 202   fTreeBeamSpot->Branch("bx", &bx, "bx/i");
 203   fTreeBeamSpot->Branch("triggered", &triggered, "triggered/s");
 204   fTreeBeamSpot->Branch("ntrkletsS", &ntrkletsS, "ntrkletsS/s");
 205   fTreeBeamSpot->Branch("xTRKnc", &xTRKnc, "xTRKnc/F");
 206   fTreeBeamSpot->Branch("yTRKnc", &yTRKnc, "yTRKnc/F");
 207   fTreeBeamSpot->Branch("zTRKnc", &zTRKnc, "zTRKnc/F");
 208   fTreeBeamSpot->Branch("ntrksTRKnc", &ntrksTRKnc, "ntrksTRKnc/s");
 209   fOutput->Add(fTreeBeamSpot);
 210
 211   Int_t nbinTrklets=29;
 212   Float_t lowTrklets[30]={-0.5,0.5,1.5,2.5,3.5,4.5,5.5,7.5,10.5,25.5,50.5,75.5,100.5,150.5,200.5,300.5,400.5,500.5,600.5,800.5,1000.5,1500.5,2000.5,2500.5,3000.5,4000.5,5000.5,6000.5,8000.5,10000.5};
 213   fhTriggeredTrklets = new TH1F("fhTriggeredTrklets","trklets dist for triggered ev.; ntrklets; entries",nbinTrklets,lowTrklets);
 214   fOutput->Add(fhTriggeredTrklets);
 215   fhSPD3DTrklets = new TH1F("fhSPD3DTrklets","trklets dist for SPD3D ev.; ntrklets; entries",nbinTrklets,lowTrklets);
 216   fOutput->Add(fhSPD3DTrklets);
 217   fhSPDZTrklets = new TH1F("fhSPDZTrklets","trklets dist for SPDZ ev.; ntrklets; entries",nbinTrklets,lowTrklets);
 218   fOutput->Add(fhSPDZTrklets);
 219   fhTRKTrklets = new TH1F("fhTRKTrklets","trklets dist for TRK ev.; ntrklets; entries",nbinTrklets,lowTrklets);
 220   fOutput->Add(fhTRKTrklets);
 221   fhTRKcTrklets = new TH1F("fhTRKcTrklets","trklets dist for TRKc ev.; ntrklets; entries",nbinTrklets,lowTrklets);
 222   fOutput->Add(fhTRKcTrklets);
 223   fhTRKncTrklets = new TH1F("fhTRKncTrklets","trklets dist for TRKnc ev.; ntrklets; entries",nbinTrklets,lowTrklets);
 224   fOutput->Add(fhTRKncTrklets);
 225   fhTPCTrklets = new TH1F("fhTPCTrklets","trklets dist for TPC ev.; ntrklets; entries",nbinTrklets,lowTrklets);
 226   fOutput->Add(fhTPCTrklets);
 227   fhTPCcTrklets = new TH1F("fhTPCcTrklets","trklets dist for TPCc ev.; ntrklets; entries",nbinTrklets,lowTrklets);
 228   fOutput->Add(fhTPCcTrklets);
 229   fhTPCncTrklets = new TH1F("fhTPCncTrklets","trklets dist for TPCnc ev.; ntrklets; entries",nbinTrklets,lowTrklets);
 230   fOutput->Add(fhTPCncTrklets);
 231
 232   Int_t nbinZreco = 16;
 233   Double_t lowZreco[17]={-15.0,-10.0,-7.0,-5,-4,-3,-2,-1,0,1,2,3,4,5,7,10,15};
 234   fhSPD3DZreco = new TH1F("fhSPD3DZreco","Zreco dist for SPD3D ev.; Zreco [cm]; entries",nbinZreco,lowZreco);
 235   fOutput->Add(fhSPD3DZreco);
 236   fhSPDZZreco = new TH1F("fhSPDZZreco","Zreco dist for SPDZ ev.; Zreco [cm]; entries",nbinZreco,lowZreco);
 237   fOutput->Add(fhSPDZZreco);
 238
 239   fhSPDContributorsPile = new TH1F("fhSPDContributorsPile","ncontributors pile up vertex; ncontributors; entries",200,-0.5,199.5);
 240   fOutput->Add(fhSPDContributorsPile);
 241
 242   fhSPDDispContributors = new TH2F("fhSPDDispContributors","ncontributors main-pile; ncontributors main; ncontributors pile",200,-0.5,199.5,200,-0.5,199.5);
 243   fOutput->Add(fhSPDDispContributors);
 244
 245   fhntrksSPDvsSPDcls = new TH2F("fhntrksSPDvsSPDcls", "ncontributors SPD3D vs number of cluster SPD; nContributors 3D; nCluster SPD1",300,0.,3000.,400,0.,8000.);
 246   fOutput->Add(fhntrksSPDvsSPDcls);
 247
 248   fhntrksZvsSPDcls = new TH2F("fhntrksZvsSPDcls", "ncontributors SPDZ vs number of cluster SPD; nContributors Z; nCluster SPD1",300,0.,3000.,400,0.,8000.);
 249   fOutput->Add(fhntrksZvsSPDcls);
 250
 251   PostData(1, fOutput);
 252
 253   return;
 254 }
 255
 256 //________________________________________________________________________
 257 void AliAnalysisTaskVertexESD::UserExec(Option_t *)
 258 {
 259   // Main loop
 260   // Called for each event
 261
 262   if (!InputEvent()) {
 263     Printf("ERROR: fESD not available");
 264     return;
 265   }
 266   AliESDEvent* esdE = (AliESDEvent*) InputEvent();
 267
 268   // Select PHYSICS events (type=7, for data) and MC events (type=0)
 269   // fCheckEventType is kFALSE if fReadMC is kTRUE, hence check is skipped
 270   if(fCheckEventType) {
 271     if(esdE->GetEventType()!=7 && esdE->GetEventType()!=0) return;
 272   }
 273
 274
 275   TArrayF mcVertex(3);
 276   mcVertex[0]=9999.; mcVertex[1]=9999.; mcVertex[2]=9999.;
 277   Float_t dNchdy=-999.;
 278   // ***********  MC info ***************
 279   if (fReadMC) {
 280     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 281     if (!eventHandler) {
 282       Printf("ERROR: Could not retrieve MC event handler");
 283       return;
 284     }
 285
 286     AliMCEvent* mcEvent = eventHandler->MCEvent();
 287     if (!mcEvent) {
 288       Printf("ERROR: Could not retrieve MC event");
 289       return;
 290     }
 291
 292     AliStack* stack = mcEvent->Stack();
 293     if (!stack) {
 294       AliDebug(AliLog::kError, "Stack not available");
 295       return;
 296     }
 297
 298     AliHeader* header = mcEvent->Header();
 299     if (!header) {
 300       AliDebug(AliLog::kError, "Header not available");
 301       return;
 302     }
 303     AliGenEventHeader* genHeader = header->GenEventHeader();
 304     genHeader->PrimaryVertex(mcVertex);
 305
 306     Int_t ngenpart = (Int_t)stack->GetNtrack();
 307     //printf("# generated particles = %d\n",ngenpart);
 308     dNchdy=0;
 309     for(Int_t ip=0; ip<ngenpart; ip++) {
 310       TParticle* part = (TParticle*)stack->Particle(ip);
 311       // keep only electorns, muons, pions, kaons and protons
 312       Int_t apdg = TMath::Abs(part->GetPdgCode());
 313       if(apdg!=11 && apdg!=13 && apdg!=211 && apdg!=321 && apdg!=2212) continue;
 314       // reject secondaries
 315       if(TMath::Sqrt((part->Vx()-mcVertex[0])*(part->Vx()-mcVertex[0])+(part->Vy()-mcVertex[1])*(part->Vy()-mcVertex[1]))>0.0010) continue;
 316       // reject incoming protons
 317       Double_t energy  = part->Energy();
 318       if(energy>900.) continue;
 319       Double_t pz = part->Pz();
 320       Double_t y = 100.;
 321       if(energy-pz+1.e-13>0.) y = 0.5*TMath::Log((energy+pz+1.e-13)/(energy-pz+1.e-13));
 322       if(TMath::Abs(y)<1.0) dNchdy += 0.5; // count 1/2 of particles in |y|<1
 323       // tracks refs
 324       TClonesArray *trefs=0;
 325       Int_t ntrefs = mcEvent->GetParticleAndTR(ip,part,trefs);
 326       if(ntrefs<0) continue;
 327       for(Int_t iref=0; iref<ntrefs; iref++) {
 328         AliTrackReference *tref = (AliTrackReference*)trefs->At(iref);
 329         if(tref->R()>10.) continue;
 330         fhTrackRefs->Fill(tref->X(),tref->Y());
 331       }
 332     }
 333     //printf("# primary particles = %7.1f\n",dNchdy);
 334   }
 335   // ***********  MC info ***************
 336
 337
 338   // Trigger
 339   //ULong64_t triggerMask;
 340   //ULong64_t spdFO = (1 << 14);
 341   //ULong64_t v0left = (1 << 10);
 342   //ULong64_t v0right = (1 << 11);
 343
 344   //triggerMask=esdE->GetTriggerMask();
 345   // MB1: SPDFO || V0L || V0R
 346   //Bool_t eventTriggered = (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)));
 347   //MB2: GFO && V0R
 348   //triggerMask & spdFO && ((triggerMask&v0left) || (triggerMask&v0right))
 349   //Bool_t eventTriggered = (triggerMask & spdFO);
 350
 351   //static AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis();
 352   Bool_t eventTriggered = 0;//triggerAnalysis->IsTriggerFired(esdE, AliTriggerAnalysis::kSPDGFO /*| AliTriggerAnalysis::kOfflineFlag*/);
 353
 354   // use response of AliPhysicsSelection
 355         eventTriggered = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerType);
 356
 357
 358   Int_t nInFile = esdE->GetEventNumberInFile();
 359
 360  const AliMultiplicity *alimult = esdE->GetMultiplicity();
 361   Int_t ntrklets=0,spd0cls=0, spd1cls=0;
 362   if(alimult) {
 363     ntrklets = alimult->GetNumberOfTracklets();
 364
 365     for(Int_t l=0;l<alimult->GetNumberOfTracklets();l++){
 366       if(alimult->GetDeltaPhi(l)<-9998.) ntrklets--;
 367     }
 368     spd0cls = alimult->GetNumberOfSingleClusters()+ntrklets;
 369         spd1cls = alimult->GetNumberOfITSClusters(1);
 370   }
 371
 372
 373   UShort_t triggered, ntrkletsS, ntrksTRKnc;
 374   UInt_t run, bx;
 375   Float_t cetTimeLHC,xTRKnc, yTRKnc, zTRKnc;
 376   fTreeBeamSpot->SetBranchAddress("run", &run);
 377   fTreeBeamSpot->SetBranchAddress("cetTimeLHC", &cetTimeLHC);
 378   fTreeBeamSpot->SetBranchAddress("bx", &bx);
 379   fTreeBeamSpot->SetBranchAddress("triggered", &triggered);
 380   fTreeBeamSpot->SetBranchAddress("ntrkletsS", &ntrkletsS);
 381   fTreeBeamSpot->SetBranchAddress("xTRKnc", &xTRKnc);
 382   fTreeBeamSpot->SetBranchAddress("yTRKnc", &yTRKnc);
 383   fTreeBeamSpot->SetBranchAddress("zTRKnc", &zTRKnc);
 384   fTreeBeamSpot->SetBranchAddress("ntrksTRKnc", &ntrksTRKnc);
 385
 386
 387   Double_t tstamp = esdE->GetTimeStamp();
 388   Float_t cetTime =(tstamp-1262304000.)+7200.;
 389
 390   Float_t cetTime1h =(tstamp-1262304000.)+3600.;
 391
 392   cetTimeLHC = (Float_t)cetTime1h;
 393
 394   Int_t ntracks = esdE->GetNumberOfTracks();
 395   Int_t nITS5or6=0,nTPCin=0,nTPCinEta09=0;
 396   //printf("Tracks # = %d\n",esdE->GetNumberOfTracks());
 397   for(Int_t itr=0; itr<ntracks; itr++) {
 398     AliESDtrack *t = esdE->GetTrack(itr);
 399     if(t->GetNcls(0)>=5) nITS5or6++;
 400     Double_t z0; t->GetZAt(0,esdE->GetMagneticField(),z0);
 401     if(t->GetNcls(1)>0 && TMath::Abs(t->GetD(0,0,esdE->GetMagneticField()))<2.8 && TMath::Abs(z0)<20) {
 402       nTPCin++;
 403       if(TMath::Abs(t->GetTgl())<1.5) nTPCinEta09++;
 404     }
 405   }
 406
 407
 408   const AliESDVertex *spdv=esdE->GetPrimaryVertexSPD();
 409   const AliESDVertex *spdvp=esdE->GetPileupVertexSPD(0);
 410   const AliESDVertex *tpcv=esdE->GetPrimaryVertexTPC();
 411   const AliESDVertex *trkv=esdE->GetPrimaryVertexTracks();
 412
 413   // fill histos
 414
 415   if(spdv) {
 416     if(spdv->GetNContributors()>0) {
 417       TString title=spdv->GetTitle();
 418       if(title.Contains("3D")) {
 419                   fhSPDVertexX->Fill(spdv->GetX());
 420                   fhSPDVertexY->Fill(spdv->GetY());
 421                   fhSPDVertexZ->Fill(spdv->GetZ());
 422
 423                   fhntrksSPDvsSPDcls->Fill(spdv->GetNContributors(), spd0cls);
 424       }
 425          if(title.Contains("Z")){
 426                  fhSPDVertexZonly->Fill(spdv->GetZ());
 427                  fhntrksZvsSPDcls->Fill(spdv->GetNContributors(), spd1cls);
 428          }
 429         }
 430   }
 431
 432   if(trkv) {
 433     if(trkv->GetNContributors()>0) {
 434       fhTRKVertexX->Fill(trkv->GetX());
 435       fhTRKVertexY->Fill(trkv->GetY());
 436       fhTRKVertexZ->Fill(trkv->GetZ());
 437     }
 438   }
 439
 440   if(tpcv) {
 441     if(tpcv->GetNContributors()>0) {
 442       fhTPCVertexX->Fill(tpcv->GetX());
 443       fhTPCVertexY->Fill(tpcv->GetY());
 444       fhTPCVertexZ->Fill(tpcv->GetZ());
 445     }
 446   }
 447
 448
 449
 450    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 451   // filling Vertex reco efficiency plots
 452
 453   if(eventTriggered ? 1. : 0.){
 454     fhTriggeredTrklets->Fill(ntrklets);
 455
 456     if(spdv){
 457       if(spdv->GetNContributors()>0.5){
 458         TString spdtitle = spdv->GetTitle();
 459         if(spdtitle.Contains("vertexer: 3D") ? 1. : 0.){
 460           fhSPD3DTrklets->Fill(ntrklets);
 461           fhSPD3DZreco->Fill(spdv->GetZ());
 462         }
 463         else{
 464           fhSPDZTrklets->Fill(ntrklets);
 465           fhSPDZZreco->Fill(spdv->GetZ());
 466         }
 467       }
 468     }
 469
 470     if(trkv){
 471       if(trkv->GetNContributors()>0.5)fhTRKTrklets->Fill(ntrklets);
 472     }
 473     if(fRecoVtxITSTPC) {
 474       AliESDVertex *trkvc = ReconstructPrimaryVertexITSTPC(kTRUE);
 475       if(trkvc->GetNContributors()>0.5)fhTRKcTrklets->Fill(ntrklets);
 476       delete trkvc; trkvc=0;
 477       AliESDVertex *trkvnc = ReconstructPrimaryVertexITSTPC(kFALSE);
 478       if(trkvnc->GetNContributors()>0.5)fhTRKncTrklets->Fill(ntrklets);
 479       delete trkvnc; trkvnc=0;
 480     }
 481
 482     if(tpcv){
 483       if(tpcv->GetNContributors()>0.5)fhTPCTrklets->Fill(ntrklets);
 484     }
 485     if(fRecoVtxTPC) {
 486         AliESDVertex *tpcvc = ReconstructPrimaryVertexTPC(kTRUE);
 487         if(tpcvc->GetNContributors()>0.5)fhTPCcTrklets->Fill(ntrklets);
 488         delete tpcvc; tpcvc=0;
 489         AliESDVertex *tpcvnc = ReconstructPrimaryVertexTPC(kFALSE);
 490         if(tpcvnc->GetNContributors()>0.5)fhTPCncTrklets->Fill(ntrklets);
 491         delete tpcvnc; tpcvnc=0;
 492       }
 493   }
 494   /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 495
 496   Float_t xpile=-999.;
 497   Float_t ypile=-999.;
 498   Float_t zpile=-999.;
 499   Float_t expile=-999.;
 500   Float_t eypile=-999.;
 501   Float_t ezpile=-999.;
 502   Int_t ntrkspile=-1;
 503   Bool_t isPileUpfromSPD=esdE->IsPileupFromSPD();
 504
 505   if(esdE->GetNumberOfPileupVerticesSPD()>0 && spdvp && spdv){
 506
 507     if(spdvp->GetNContributors()>=2) {
 508
 509       xpile=spdvp->GetX();
 510       expile=spdvp->GetXRes();
 511       ypile=spdvp->GetY();
 512       eypile=spdvp->GetYRes();
 513       zpile=spdvp->GetZ();
 514       ezpile=spdvp->GetZRes();
 515       ntrkspile=spdvp->GetNContributors();
 516
 517       if (esdE->IsPileupFromSPD(2,0.,3.,2.,5.)){
 518
 519         fhSPDVertexDiffZPileContr2->Fill(spdv->GetZ()-spdvp->GetZ());
 520         if(spdvp->GetNContributors()>=3) {fhSPDVertexDiffZPileContr3->Fill(spdv->GetZ()-spdvp->GetZ());}
 521         if(spdvp->GetNContributors()>=4) {fhSPDVertexDiffZPileContr4->Fill(spdv->GetZ()-spdvp->GetZ());}
 522         if(spdvp->GetNContributors()>=5) {fhSPDVertexDiffZPileContr5->Fill(spdv->GetZ()-spdvp->GetZ());}
 523
 524       }//end IsPileUpFromSPD
 525
 526       if (isPileUpfromSPD){
 527
 528         fhSPDVertexXPile->Fill(spdvp->GetX());
 529             fhSPDVertexYPile->Fill(spdvp->GetY());
 530             fhSPDVertexZPile->Fill(spdvp->GetZ());
 531         fhSPDContributorsPile->Fill(spdvp->GetNContributors());
 532         fhSPDDispContributors->Fill(spdv->GetNContributors(),spdvp->GetNContributors());
 533         fhSPDVertexDiffZPileDefault->Fill(spdv->GetZ()-spdvp->GetZ());
 534
 535       }
 536     }
 537   }
 538
 539   // fill ntuple
 540   Int_t isize=89;
 541   Float_t xnt[89]; for(Int_t iii=0;iii<isize;iii++) xnt[iii]=0.;
 542
 543   Int_t isizeSecNt=9;
 544   Float_t secnt[9]; for(Int_t iii=0;iii<isizeSecNt;iii++) secnt[iii]=0.;
 545
 546   Int_t index=0;
 547   Int_t indexSecNt=0;
 548
 549   xnt[index++]=(Float_t)esdE->GetRunNumber();
 550   secnt[indexSecNt++]=(Float_t)esdE->GetRunNumber();
 551   run = (Int_t)esdE->GetRunNumber();
 552
 553   xnt[index++]=cetTime; //(Float_t)esdE->GetTimeStamp();
 554   //secnt[indexSecNt++]=cetTime;
 555   secnt[indexSecNt++]=cetTime1h;
 556
 557   xnt[index++]=(Float_t)esdE->GetBunchCrossNumber();
 558   secnt[indexSecNt++]=(Float_t)esdE->GetBunchCrossNumber();
 559   bx = (Int_t)esdE->GetBunchCrossNumber();
 560
 561   xnt[index++]=(eventTriggered ? 1. : 0.);
 562   secnt[indexSecNt++]=(eventTriggered ? 1. : 0.);
 563   triggered = (UShort_t)(eventTriggered ? 1 : 0);
 564
 565   xnt[index++]=(Float_t)dNchdy;
 566
 567   xnt[index++]=(Float_t)esdE->GetDiamondX();
 568   xnt[index++]=(Float_t)esdE->GetDiamondY();
 569   xnt[index++]=(Float_t)esdE->GetDiamondZ();
 570
 571   xnt[index++]=(Float_t)(TMath::Sqrt(esdE->GetSigma2DiamondX()));
 572   xnt[index++]=(Float_t)(TMath::Sqrt(esdE->GetSigma2DiamondY()));
 573   xnt[index++]=(Float_t)(TMath::Sqrt(esdE->GetSigma2DiamondZ()));
 574
 575
 576   xnt[index++]=mcVertex[0];
 577   xnt[index++]=mcVertex[1];
 578   xnt[index++]=mcVertex[2];
 579
 580   xnt[index++]=spdv->GetX();
 581   xnt[index++]=spdv->GetXRes();
 582   xnt[index++]=spdv->GetY();
 583   xnt[index++]=spdv->GetYRes();
 584   xnt[index++]=spdv->GetZ();
 585   xnt[index++]=spdv->GetZRes();
 586   xnt[index++]=spdv->GetNContributors();
 587   TString spdtitle = spdv->GetTitle();
 588   xnt[index++]=(spdtitle.Contains("vertexer: 3D") ? 1. : 0.);
 589   xnt[index++]=spdv->GetDispersion();
 590
 591   xnt[index++]=tpcv->GetX();
 592   xnt[index++]=tpcv->GetXRes();
 593   xnt[index++]=tpcv->GetY();
 594   xnt[index++]=tpcv->GetYRes();
 595   xnt[index++]=tpcv->GetZ();
 596   xnt[index++]=tpcv->GetZRes();
 597   xnt[index++]=tpcv->GetNContributors();
 598   TString tpctitle = tpcv->GetTitle();
 599   xnt[index++]=(tpctitle.Contains("WithConstraint") ? 1. : 0.);
 600
 601   xnt[index++]=trkv->GetX();
 602   xnt[index++]=trkv->GetXRes();
 603   xnt[index++]=trkv->GetY();
 604   xnt[index++]=trkv->GetYRes();
 605   xnt[index++]=trkv->GetZ();
 606   xnt[index++]=trkv->GetZRes();
 607   xnt[index++]=trkv->GetNContributors();// tpccontrorig;
 608   TString trktitle = trkv->GetTitle();
 609   xnt[index++]=(trktitle.Contains("WithConstraint") ? 1. : 0.);
 610
 611   xnt[index++]=float(ntrklets);
 612   secnt[indexSecNt++]=float(ntrklets);
 613   ntrkletsS = (UShort_t)ntrklets;
 614
 615   xnt[index++]=float(ntracks);
 616   xnt[index++]=float(nITS5or6);
 617
 618   xnt[index++]=float(nTPCin);
 619   xnt[index++]=float(nTPCinEta09);
 620
 621   xnt[index++]=spd0cls;
 622
 623   xnt[index++]=Float_t(isPileUpfromSPD);
 624   xnt[index++]=xpile;
 625   xnt[index++]=expile;
 626   xnt[index++]=ypile;
 627   xnt[index++]=eypile;
 628   xnt[index++]=zpile;
 629   xnt[index++]=ezpile;
 630   xnt[index++]=(Float_t)ntrkspile;
 631
 632
 633
 634   // add recalculated vertices TRK and TPC
 635
 636   if(fRecoVtxTPC) {
 637     AliESDVertex *tpcvnc = ReconstructPrimaryVertexTPC(kFALSE);
 638     xnt[index++]=tpcvnc->GetX();
 639     xnt[index++]=tpcvnc->GetXRes();
 640     xnt[index++]=tpcvnc->GetY();
 641     xnt[index++]=tpcvnc->GetYRes();
 642     xnt[index++]=tpcvnc->GetZ();
 643     xnt[index++]=tpcvnc->GetZRes();
 644     xnt[index++]=tpcvnc->GetNContributors();
 645     delete tpcvnc; tpcvnc=0;
 646
 647     AliESDVertex *tpcvc = ReconstructPrimaryVertexTPC(kTRUE);
 648     xnt[index++]=tpcvc->GetX();
 649     xnt[index++]=tpcvc->GetXRes();
 650     xnt[index++]=tpcvc->GetY();
 651     xnt[index++]=tpcvc->GetYRes();
 652     xnt[index++]=tpcvc->GetZ();
 653     xnt[index++]=tpcvc->GetZRes();
 654     xnt[index++]=tpcvc->GetNContributors();
 655     delete tpcvc; tpcvc=0;
 656   } else index+=14;
 657
 658   if(fRecoVtxITSTPC) {
 659     AliESDVertex *trkvnc = ReconstructPrimaryVertexITSTPC(kFALSE);
 660     xnt[index++]=trkvnc->GetX();
 661     xnt[index++]=trkvnc->GetXRes();
 662     xnt[index++]=trkvnc->GetY();
 663     xnt[index++]=trkvnc->GetYRes();
 664     xnt[index++]=trkvnc->GetZ();
 665     xnt[index++]=trkvnc->GetZRes();
 666     xnt[index++]=trkvnc->GetNContributors();
 667
 668     secnt[indexSecNt++]=trkvnc->GetX();
 669     secnt[indexSecNt++]=trkvnc->GetY();
 670     secnt[indexSecNt++]=trkvnc->GetZ();
 671     secnt[indexSecNt++]=trkvnc->GetNContributors();
 672
 673     xTRKnc = (Float_t)trkvnc->GetX();
 674     yTRKnc = (Float_t)trkvnc->GetY();
 675     zTRKnc = (Float_t)trkvnc->GetZ();
 676     ntrksTRKnc = (trkvnc->GetNContributors()<0 ? 0 : (UShort_t)trkvnc->GetNContributors());
 677
 678     delete trkvnc; trkvnc=0;
 679
 680
 681     AliESDVertex *trkvc = ReconstructPrimaryVertexITSTPC(kTRUE);
 682     xnt[index++]=trkvc->GetX();
 683     xnt[index++]=trkvc->GetXRes();
 684     xnt[index++]=trkvc->GetY();
 685     xnt[index++]=trkvc->GetYRes();
 686     xnt[index++]=trkvc->GetZ();
 687     xnt[index++]=trkvc->GetZRes();
 688     xnt[index++]=trkvc->GetNContributors();
 689     delete trkvc; trkvc=0;
 690   } else index+=14;
 691
 692   if(fRecoVtxITSTPCHalfEvent) {
 693     AliESDVertex *trkvncodd  = ReconstructPrimaryVertexITSTPC(kFALSE,1);
 694     AliESDVertex *trkvnceven = ReconstructPrimaryVertexITSTPC(kFALSE,2);
 695     if(trkvncodd->GetNContributors()>0 &&
 696        trkvnceven->GetNContributors()>0) {
 697       xnt[index++]=trkvncodd->GetX()-trkvnceven->GetX();
 698       xnt[index++]=trkvncodd->GetY()-trkvnceven->GetY();
 699       xnt[index++]=trkvncodd->GetZ()-trkvnceven->GetZ();
 700       xnt[index++]=TMath::Sqrt(trkvncodd->GetXRes()*trkvncodd->GetXRes()+trkvnceven->GetXRes()*trkvnceven->GetXRes());
 701       xnt[index++]=TMath::Sqrt(trkvncodd->GetYRes()*trkvncodd->GetYRes()+trkvnceven->GetYRes()*trkvnceven->GetYRes());
 702       xnt[index++]=TMath::Sqrt(trkvncodd->GetZRes()*trkvncodd->GetZRes()+trkvnceven->GetZRes()*trkvnceven->GetZRes());
 703       xnt[index++]=trkvnceven->GetNContributors();
 704       xnt[index++]=trkvncodd->GetNContributors();
 705     } else {
 706       xnt[index++]=0.;
 707       xnt[index++]=0.;
 708       xnt[index++]=0.;
 709       xnt[index++]=0.;
 710       xnt[index++]=0.;
 711       xnt[index++]=0.;
 712       xnt[index++]=-1;
 713       xnt[index++]=-1;
 714     }
 715     delete trkvncodd; trkvncodd=0;
 716     delete trkvnceven; trkvnceven=0;
 717   } else index+=8;
 718
 719
 720   if(index>isize) printf("AliAnalysisTaskVertexESD: ERROR, index!=isize\n");
 721   if(fFillNtuple) fNtupleVertexESD->Fill(xnt);
 722
 723   if(indexSecNt>isizeSecNt) printf("AliAnalysisTaskVertexESD: ERROR, indexSecNt!=isizeSecNt\n");
 724
 725   // if(indexTree>isizeTree) printf("AliAnalysisTaskVertexESD: ERROR, indexTree!=isizeTree\n");
 726   // only every second event (to reduce output size)
 727   if(fFillTreeBeamSpot && (nInFile%2 == 0)) fTreeBeamSpot->Fill();
 728
 729
 730   // Post the data already here
 731   PostData(1, fOutput);
 732
 733   return;
 734 }
 735
 736 //________________________________________________________________________
 737 void AliAnalysisTaskVertexESD::Terminate(Option_t *)
 738 {
 739   // Draw result to the screen
 740   // Called once at the end of the query
 741   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 742   if (!fOutput) {
 743     Printf("ERROR: fOutput not available");
 744     return;
 745   }
 746
 747   //////////////////////////////////////////////////////
 748   /*
 749   TH1F *fhTriggeredTrklets=(TH1F*)fOutput->FindObject("fhTriggeredTrklets");
 750   TH1F *fhSPDZTrklets=(TH1F*)fOutput->FindObject("fhSPDZTrklets");
 751   TH1F *fhSPD3DTrklets=(TH1F*)fOutput->FindObject("fhSPD3DTrklets");
 752   TH1F *fhTRKTrklets=(TH1F*)fOutput->FindObject("fhTRKTrklets");
 753   TH1F *fhTRKcTrklets=(TH1F*)fOutput->FindObject("fhTRKcTrklets");
 754   TH1F *fhTRKncTrklets=(TH1F*)fOutput->FindObject("fhTRKncTrklets");
 755   TH1F *fhSPDZZreco=(TH1F*)fOutput->FindObject("fhSPDZZreco");
 756   TH1F *fhSPD3DZreco=(TH1F*)fOutput->FindObject("fhSPD3DZreco");
 757
 758   TGraphAsymmErrors *fhSPDZEffTrklets=new TGraphAsymmErrors(fhSPDZTrklets,fhTriggeredTrklets,"w");
 759   fhSPDZEffTrklets->SetName("fhSPDZEffTrklets");
 760   fhSPDZEffTrklets->SetDrawOption("AP");
 761   TGraphAsymmErrors *fhSPD3DEffTrklets=new TGraphAsymmErrors(fhSPD3DTrklets,fhTriggeredTrklets,"w");
 762   fhSPD3DEffTrklets->SetName("fhSPD3DEffTrklets");
 763   TH1F * fhSPDOverallTrklets=(TH1F*)fhSPDZTrklets->Clone("fhSPDOverallTrklets");
 764   fhSPDOverallTrklets->Add(fhSPD3DTrklets);
 765   TGraphAsymmErrors *fhSPDOverallEffTrklets=new TGraphAsymmErrors(fhSPDOverallTrklets,fhTriggeredTrklets,"w");
 766   fhSPDOverallEffTrklets->SetName("fhSPDOverallEffTrklets");
 767   TGraphAsymmErrors *fhTRKEffTrklets=new TGraphAsymmErrors(fhTRKTrklets,fhTriggeredTrklets,"w");
 768   fhTRKEffTrklets->SetName("fhTRKEffTrklets");
 769   TGraphAsymmErrors *fhTRKcEffTrklets=new TGraphAsymmErrors(fhTRKcTrklets,fhTriggeredTrklets,"w");
 770   fhTRKcEffTrklets->SetName("fhTRKcEffTrklets");
 771   TGraphAsymmErrors *fhTRKncEffTrklets=new TGraphAsymmErrors(fhTRKncTrklets,fhTriggeredTrklets,"w");
 772   fhTRKncEffTrklets->SetName("fhTRKncEffTrklets");
 773   TH1F * fhSPDOverallZreco=(TH1F*)fhSPDZZreco->Clone("fhSPDOverallZreco");
 774   fhSPDOverallZreco->Add(fhSPD3DZreco);
 775   TGraphAsymmErrors *fhSPDEffZreco=new TGraphAsymmErrors(fhSPD3DZreco,fhSPDOverallZreco,"w");
 776   fhSPDEffZreco->SetName("fhSPDEffZreco");
 777
 778   TH1F *fhEff = new TH1F("hEff","hEff",6,0.5,6.5);
 779   Int_t count=1;
 780   if(fhSPDZTrklets->GetEntries()!=0 && fhTriggeredTrklets->GetEntries()!=0){
 781     fhEff->Fill(count,fhSPDZTrklets->GetEntries()/fhTriggeredTrklets->GetEntries());
 782     fhEff->SetBinError(count,fhEff->GetBinContent(count)*TMath::Sqrt(1/fhSPDZTrklets->GetEntries()+1/fhTriggeredTrklets->GetEntries()));
 783   }
 784   fhEff->GetXaxis()->SetBinLabel(count,"SPDZ");
 785
 786   count++;
 787   if(fhSPD3DTrklets->GetEntries()!=0 && fhTriggeredTrklets->GetEntries()!=0){
 788     fhEff->Fill(count,fhSPD3DTrklets->GetEntries()/fhTriggeredTrklets->GetEntries());
 789     fhEff->SetBinError(count,fhEff->GetBinContent(count)*TMath::Sqrt(1/fhSPD3DTrklets->GetEntries()+1/fhTriggeredTrklets->GetEntries()));
 790   }
 791   fhEff->GetXaxis()->SetBinLabel(count,"SPD3D");
 792
 793   count++;
 794   if(fhSPDOverallTrklets->GetEntries()!=0 && fhTriggeredTrklets->GetEntries()!=0){
 795     fhEff->Fill(count,fhSPDOverallTrklets->GetEntries()/fhTriggeredTrklets->GetEntries());
 796     fhEff->SetBinError(count,fhEff->GetBinContent(count)*TMath::Sqrt(1/fhSPDOverallTrklets->GetEntries()+1/fhTriggeredTrklets->GetEntries()));
 797   }
 798   fhEff->GetXaxis()->SetBinLabel(count,"SPD Overall");
 799
 800   count++;
 801   if(fhTRKTrklets->GetEntries()!=0 && fhTriggeredTrklets->GetEntries()!=0){
 802     fhEff->Fill(count,fhTRKTrklets->GetEntries()/fhTriggeredTrklets->GetEntries());
 803     fhEff->SetBinError(count,fhEff->GetBinContent(count)*TMath::Sqrt(1/fhTRKTrklets->GetEntries()+1/fhTriggeredTrklets->GetEntries()));
 804   }
 805   fhEff->GetXaxis()->SetBinLabel(count,"TRK");
 806
 807   count++;
 808   if(fhTRKcTrklets->GetEntries()!=0 && fhTriggeredTrklets->GetEntries()!=0){
 809     fhEff->Fill(count,fhTRKcTrklets->GetEntries()/fhTriggeredTrklets->GetEntries());
 810     fhEff->SetBinError(count,fhEff->GetBinContent(count)*TMath::Sqrt(1/fhTRKcTrklets->GetEntries()+1/fhTriggeredTrklets->GetEntries()));
 811   }
 812   fhEff->GetXaxis()->SetBinLabel(count,"TRKc");
 813
 814   count++;
 815   if(fhTRKncTrklets->GetEntries()!=0 && fhTriggeredTrklets->GetEntries()!=0){
 816     fhEff->Fill(count,fhTRKncTrklets->GetEntries()/fhTriggeredTrklets->GetEntries());
 817     fhEff->SetBinError(count,fhEff->GetBinContent(count)*TMath::Sqrt(1/fhTRKncTrklets->GetEntries()+1/fhTriggeredTrklets->GetEntries()));
 818   }
 819   fhEff->GetXaxis()->SetBinLabel(count,"TRKnc");
 820
 821   count++;
 822   fhEff->Print("all");
 823
 824   TFile* fileEff = new TFile("VtxEff.root","recreate");
 825   fhSPDZEffTrklets->Write();
 826   fhSPD3DEffTrklets->Write();
 827   fhSPDOverallEffTrklets->Write();
 828   fhTRKEffTrklets->Write();
 829   fhTRKcEffTrklets->Write();
 830   fhTRKncEffTrklets->Write();
 831   fhSPDEffZreco->Write();
 832   fhEff->Write();
 833   fileEff->Close();
 834   delete fileEff;
 835
 836   /////////////////////////////////////////
 837   */
 838
 839
 840   if (!fNtupleVertexESD){
 841     Printf("ERROR: fNtuple not available");
 842     return;
 843   }
 844
 845   fNtupleVertexESD = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleVertexESD"));
 846
 847
 848   return;
 849 }
 850
 851 //_________________________________________________________________________
 852 AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexTPC(Bool_t constr) const {
 853   // On the fly reco of TPC vertex from ESD
 854   AliESDEvent* evt = (AliESDEvent*) fInputEvent;
 855   AliVertexerTracks vertexer(evt->GetMagneticField());
 856   if(evt->GetNumberOfTracks()<500) {
 857     vertexer.SetTPCMode(); // defaults
 858   } else {
 859     vertexer.SetTPCMode(0.1,1.0,5.0,10,1,3.,0.1,1.5,3.,30.,1,1);// PbPb
 860   }
 861   Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy);
 862   Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0};
 863   Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
 864   if (cov[0]<1e-10) cov[0] = 1e-2; // dummy error of 1mm
 865   if (cov[2]<1e-10) cov[2] = 1e-2; // dummy error of 1mm
 866   AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1);
 867   vertexer.SetVtxStart(initVertex);
 868   delete initVertex;
 869   if(!constr) vertexer.SetConstraintOff();
 870
 871   return vertexer.FindPrimaryVertex(evt);
 872 }
 873
 874 //_________________________________________________________________________
 875 AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexITSTPC(Bool_t constr,Int_t mode) const {
 876   // On the fly reco of ITS+TPC vertex from ESD
 877   // mode=0 use all tracks
 878   // mode=1 use odd-number tracks
 879   // mode=2 use even-number tracks
 880
 881   AliESDEvent* evt = (AliESDEvent*) fInputEvent;
 882   AliVertexerTracks vertexer(evt->GetMagneticField());
 883   if(evt->GetNumberOfTracks()<500) {
 884     vertexer.SetITSMode(); // defaults
 885     vertexer.SetMinClusters(3); // default is 5
 886   } else {
 887     vertexer.SetITSMode(0.1,0.1,0.5,5,1,3.,100.,1000.,3.,30.,1,1);// PbPb
 888   }
 889   //vertexer.SetITSpureSA();
 890   Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy);
 891   Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0};
 892   Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
 893   if (cov[0]<1e-10) cov[0] = 1e-2; // dummy error of 1mm
 894   if (cov[2]<1e-10) cov[2] = 1e-2; // dummy error of 1mm
 895   AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1);
 896   vertexer.SetVtxStart(initVertex);
 897   delete initVertex;
 898   if(!constr) vertexer.SetConstraintOff();
 899
 900   // use only ITS-TPC or only ITS-SA tracks
 901   if(fOnlyITSTPCTracks || fOnlyITSSATracks || mode>0) {
 902     Int_t iskip=0;
 903     Int_t ntracks = evt->GetNumberOfTracks();
 904     Int_t *skip = new Int_t[ntracks];
 905     for(Int_t i=0;i<ntracks;i++) skip[i]=-1;
 906     for(Int_t itr=0;itr<ntracks; itr++) {
 907       AliESDtrack* track = evt->GetTrack(itr);
 908       if(fOnlyITSTPCTracks && track->GetNcls(1)==0) { // skip ITSSA
 909         skip[iskip++]=itr;
 910         continue;
 911       }
 912       if(fOnlyITSSATracks && track->GetNcls(1)>0) { // skip ITSTPC
 913         skip[iskip++]=itr;
 914         continue;
 915       }
 916       if(mode==1 && itr%2==0) skip[iskip++]=itr;
 917       if(mode==2 && itr%2==1) skip[iskip++]=itr;
 918     }
 919     vertexer.SetSkipTracks(iskip,skip);
 920     delete [] skip; skip=NULL;
 921   }
 922
 923   return vertexer.FindPrimaryVertex(evt);
 924 }