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 AliAnalysisTaskSE(name),
  65 fCheckEventType(kTRUE),
  66 fReadMC(kFALSE),
  67 fRecoVtxTPC(kFALSE),
  68 fRecoVtxITSTPC(kFALSE),
  69 fRecoVtxITSTPCHalfEvent(kFALSE),
  70 fOnlyITSTPCTracks(kFALSE),
  71 fOnlyITSSATracks(kFALSE),
  72 fFillNtuple(kFALSE),
  73 fESD(0),
  74 fOutput(0),
  75 fNtupleVertexESD(0),
  76 fhSPDVertexX(0),
  77 fhSPDVertexY(0),
  78 fhSPDVertexZ(0),
  79 fhTRKVertexX(0),
  80 fhTRKVertexY(0),
  81 fhTRKVertexZ(0),
  82 fhTPCVertexX(0),
  83 fhTPCVertexY(0),
  84 fhTPCVertexZ(0),
  85 fhTrackRefs(0)
  86 {
  87   // Constructor
  88
  89   // Define input and output slots here
  90   // Output slot #0 writes into a TList container
  91   DefineOutput(1, TList::Class());  //My private output
  92 }
  93 //________________________________________________________________________
  94 AliAnalysisTaskVertexESD::~AliAnalysisTaskVertexESD()
  95 {
  96   // Destructor
  97
  98   // histograms are in the output list and deleted when the output
  99   // list is deleted by the TSelector dtor
 100
 101   if (fOutput) {
 102     delete fOutput;
 103     fOutput = 0;
 104   }
 105 }
 106
 107
 108 //________________________________________________________________________
 109 void AliAnalysisTaskVertexESD::UserCreateOutputObjects()
 110 {
 111   // Create histograms
 112   // Called once
 113
 114   // Several histograms are more conveniently managed in a TList
 115   fOutput = new TList;
 116   fOutput->SetOwner();
 117
 118   fNtupleVertexESD = new TNtuple("fNtupleVertexESD","vertices","run:tstamp:bunchcross:triggered:dndygen: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: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");
 119
 120   fOutput->Add(fNtupleVertexESD);
 121
 122   fhSPDVertexX = new TH1F("fhSPDVertexX","SPDVertex x; x vertex [cm]; events",200,-1,1);
 123   fOutput->Add(fhSPDVertexX);
 124   fhSPDVertexY = new TH1F("fhSPDVertexY","SPDVertex y; y vertex [cm]; events",200,-1,1);
 125   fOutput->Add(fhSPDVertexY);
 126   fhSPDVertexZ = new TH1F("fhSPDVertexZ","SPDVertex z; z vertex [cm]; events",200,-20,20);
 127   fOutput->Add(fhSPDVertexZ);
 128   fhTRKVertexX = new TH1F("fhTRKVertexX","TRKVertex x; x vertex [cm]; events",200,-1,1);
 129   fOutput->Add(fhTRKVertexX);
 130   fhTRKVertexY = new TH1F("fhTRKVertexY","TRKVertex y; y vertex [cm]; events",200,-1,1);
 131   fOutput->Add(fhTRKVertexY);
 132   fhTRKVertexZ = new TH1F("fhTRKVertexZ","TRKVertex z; z vertex [cm]; events",200,-20,20);
 133   fOutput->Add(fhTRKVertexZ);
 134   fhTPCVertexX = new TH1F("fhTPCVertexX","TPCVertex x; x vertex [cm]; events",200,-3,3);
 135   fOutput->Add(fhTPCVertexX);
 136   fhTPCVertexY = new TH1F("fhTPCVertexY","TPCVertex y; y vertex [cm]; events",200,-3,3);
 137   fOutput->Add(fhTPCVertexY);
 138   fhTPCVertexZ = new TH1F("fhTPCVertexZ","TPCVertex z; z vertex [cm]; events",200,-20,20);
 139   fOutput->Add(fhTPCVertexZ);
 140
 141   fhTrackRefs = new TH2F("fhTrackRefs","Track references; x; y",1000,-4,4,1000,-4,4);
 142   fOutput->Add(fhTrackRefs);
 143
 144   return;
 145 }
 146
 147 //________________________________________________________________________
 148 void AliAnalysisTaskVertexESD::UserExec(Option_t *)
 149 {
 150   // Main loop
 151   // Called for each event
 152
 153   if (!InputEvent()) {
 154     Printf("ERROR: fESD not available");
 155     return;
 156   }
 157   AliESDEvent* esdE = (AliESDEvent*) InputEvent();
 158
 159   if(fCheckEventType && (esdE->GetEventType())!=7) return;
 160
 161
 162   TArrayF mcVertex(3);
 163   mcVertex[0]=9999.; mcVertex[1]=9999.; mcVertex[2]=9999.;
 164   Float_t dNchdy=-999.;
 165   // ***********  MC info ***************
 166   if (fReadMC) {
 167     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 168     if (!eventHandler) {
 169       Printf("ERROR: Could not retrieve MC event handler");
 170       return;
 171     }
 172
 173     AliMCEvent* mcEvent = eventHandler->MCEvent();
 174     if (!mcEvent) {
 175       Printf("ERROR: Could not retrieve MC event");
 176       return;
 177     }
 178
 179     AliStack* stack = mcEvent->Stack();
 180     if (!stack) {
 181       AliDebug(AliLog::kError, "Stack not available");
 182       return;
 183     }
 184
 185     AliHeader* header = mcEvent->Header();
 186     if (!header) {
 187       AliDebug(AliLog::kError, "Header not available");
 188       return;
 189     }
 190     AliGenEventHeader* genHeader = header->GenEventHeader();
 191     genHeader->PrimaryVertex(mcVertex);
 192
 193     Int_t ngenpart = (Int_t)stack->GetNtrack();
 194     //printf("# generated particles = %d\n",ngenpart);
 195     dNchdy=0;
 196     for(Int_t ip=0; ip<ngenpart; ip++) {
 197       TParticle* part = (TParticle*)stack->Particle(ip);
 198       // keep only electorns, muons, pions, kaons and protons
 199       Int_t apdg = TMath::Abs(part->GetPdgCode());
 200       if(apdg!=11 && apdg!=13 && apdg!=211 && apdg!=321 && apdg!=2212) continue;
 201       // reject secondaries
 202       if(TMath::Sqrt((part->Vx()-mcVertex[0])*(part->Vx()-mcVertex[0])+(part->Vy()-mcVertex[1])*(part->Vy()-mcVertex[1]))>0.0010) continue;
 203       // reject incoming protons
 204       Double_t energy  = part->Energy();
 205       if(energy>900.) continue;
 206       Double_t pz = part->Pz();
 207       Double_t y = 0.5*TMath::Log((energy+pz+1.e-13)/(energy-pz+1.e-13));
 208       if(TMath::Abs(y)<1.0) dNchdy += 0.5; // count 1/2 of particles in |y|<1
 209       // tracks refs
 210       TClonesArray *trefs=0;
 211       Int_t ntrefs = mcEvent->GetParticleAndTR(ip,part,trefs);
 212       if(ntrefs<0) continue;
 213       for(Int_t iref=0; iref<ntrefs; iref++) {
 214         AliTrackReference *tref = (AliTrackReference*)trefs->At(iref);
 215         if(tref->R()>10.) continue;
 216         fhTrackRefs->Fill(tref->X(),tref->Y());
 217       }
 218     }
 219     //printf("# primary particles = %7.1f\n",dNchdy);
 220   }
 221   // ***********  MC info ***************
 222
 223
 224   // Trigger
 225   //ULong64_t triggerMask;
 226   //ULong64_t spdFO = (1 << 14);
 227   //ULong64_t v0left = (1 << 10);
 228   //ULong64_t v0right = (1 << 11);
 229
 230   //triggerMask=esdE->GetTriggerMask();
 231   // MB1: SPDFO || V0L || V0R
 232   //Bool_t eventTriggered = (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)));
 233   //MB2: GFO && V0R
 234   //triggerMask & spdFO && ((triggerMask&v0left) || (triggerMask&v0right))
 235   //Bool_t eventTriggered = (triggerMask & spdFO);
 236
 237   //static AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis();
 238   Bool_t eventTriggered = 0;//triggerAnalysis->IsTriggerFired(esdE, AliTriggerAnalysis::kSPDGFO /*| AliTriggerAnalysis::kOfflineFlag*/);
 239
 240   // use response of AliPhysicsSelection
 241   eventTriggered = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 242
 243   Int_t ntracks = esdE->GetNumberOfTracks();
 244   Int_t nITS5or6=0,nTPCin=0,nTPCinEta09=0;
 245   //printf("Tracks # = %d\n",esdE->GetNumberOfTracks());
 246   for(Int_t itr=0; itr<ntracks; itr++) {
 247     AliESDtrack *t = esdE->GetTrack(itr);
 248     if(t->GetNcls(0)>=5) nITS5or6++;
 249     Double_t z0; t->GetZAt(0,esdE->GetMagneticField(),z0);
 250     if(t->GetNcls(1)>0 && TMath::Abs(t->GetD(0,0,esdE->GetMagneticField()))<2.8 && TMath::Abs(z0)<20) {
 251       nTPCin++;
 252       if(TMath::Abs(t->GetTgl())<1.5) nTPCinEta09++;
 253     }
 254   }
 255
 256
 257   const AliESDVertex *spdv=esdE->GetPrimaryVertexSPD();
 258   const AliESDVertex *spdvp=esdE->GetPileupVertexSPD(0);
 259   const AliESDVertex *tpcv=esdE->GetPrimaryVertexTPC();
 260   const AliESDVertex *trkv=esdE->GetPrimaryVertexTracks();
 261
 262   const AliMultiplicity *alimult = esdE->GetMultiplicity();
 263   Int_t ntrklets=0,spd0cls=0;
 264   if(alimult) {
 265     ntrklets = alimult->GetNumberOfTracklets();
 266     for(Int_t l=0;l<alimult->GetNumberOfTracklets();l++){
 267       if(alimult->GetDeltaPhi(l)<-9998.) ntrklets--;
 268     }
 269     spd0cls = alimult->GetNumberOfSingleClusters()+ntrklets;
 270   }
 271
 272   // fill histos
 273
 274   if(spdv) {
 275     if(spdv->GetNContributors()>0) {
 276       TString title=spdv->GetTitle();
 277       if(title.Contains("3D")) {
 278         fhSPDVertexX->Fill(spdv->GetXv());
 279         fhSPDVertexY->Fill(spdv->GetYv());
 280       }
 281       fhSPDVertexZ->Fill(spdv->GetZv());
 282     }
 283   }
 284
 285   if(trkv) {
 286     if(trkv->GetNContributors()>0) {
 287       fhTRKVertexX->Fill(trkv->GetXv());
 288       fhTRKVertexY->Fill(trkv->GetYv());
 289       fhTRKVertexZ->Fill(trkv->GetZv());
 290     }
 291   }
 292
 293   if(tpcv) {
 294     if(tpcv->GetNContributors()>0) {
 295       fhTPCVertexX->Fill(tpcv->GetXv());
 296       fhTPCVertexY->Fill(tpcv->GetYv());
 297       fhTPCVertexZ->Fill(tpcv->GetZv());
 298     }
 299   }
 300
 301   Float_t xpile=-999.;
 302   Float_t ypile=-999.;
 303   Float_t zpile=-999.;
 304   Float_t expile=-999.;
 305   Float_t eypile=-999.;
 306   Float_t ezpile=-999.;
 307   Int_t ntrkspile=-1;
 308
 309   if(esdE->GetNumberOfPileupVerticesSPD()>0 && spdvp){
 310     xpile=spdvp->GetXv();
 311     expile=spdvp->GetXRes();
 312     ypile=spdvp->GetYv();
 313     eypile=spdvp->GetYRes();
 314     zpile=spdvp->GetZv();
 315     ezpile=spdvp->GetZRes();
 316     ntrkspile=spdvp->GetNContributors();
 317   }
 318
 319   // fill ntuple
 320   Int_t isize=82;
 321   Float_t xnt[82]; for(Int_t iii=0;iii<isize;iii++) xnt[iii]=0.;
 322
 323   Int_t index=0;
 324
 325   xnt[index++]=(Float_t)esdE->GetRunNumber();
 326   xnt[index++]=(Float_t)esdE->GetTimeStamp();
 327   xnt[index++]=(Float_t)esdE->GetBunchCrossNumber();
 328   xnt[index++]=(eventTriggered ? 1. : 0.);
 329
 330   xnt[index++]=(Float_t)dNchdy;
 331
 332
 333   xnt[index++]=mcVertex[0];
 334   xnt[index++]=mcVertex[1];
 335   xnt[index++]=mcVertex[2];
 336
 337   xnt[index++]=spdv->GetXv();
 338   xnt[index++]=spdv->GetXRes();
 339   xnt[index++]=spdv->GetYv();
 340   xnt[index++]=spdv->GetYRes();
 341   xnt[index++]=spdv->GetZv();
 342   xnt[index++]=spdv->GetZRes();
 343   xnt[index++]=spdv->GetNContributors();
 344   TString spdtitle = spdv->GetTitle();
 345   xnt[index++]=(spdtitle.Contains("vertexer: 3D") ? 1. : 0.);
 346   xnt[index++]=spdv->GetDispersion();
 347
 348   xnt[index++]=tpcv->GetXv();
 349   xnt[index++]=tpcv->GetXRes();
 350   xnt[index++]=tpcv->GetYv();
 351   xnt[index++]=tpcv->GetYRes();
 352   xnt[index++]=tpcv->GetZv();
 353   xnt[index++]=tpcv->GetZRes();
 354   xnt[index++]=tpcv->GetNContributors();
 355   TString tpctitle = tpcv->GetTitle();
 356   xnt[index++]=(tpctitle.Contains("WithConstraint") ? 1. : 0.);
 357
 358   xnt[index++]=trkv->GetXv();
 359   xnt[index++]=trkv->GetXRes();
 360   xnt[index++]=trkv->GetYv();
 361   xnt[index++]=trkv->GetYRes();
 362   xnt[index++]=trkv->GetZv();
 363   xnt[index++]=trkv->GetZRes();
 364   xnt[index++]=trkv->GetNContributors();// tpccontrorig;
 365   TString trktitle = trkv->GetTitle();
 366   xnt[index++]=(trktitle.Contains("WithConstraint") ? 1. : 0.);
 367
 368   xnt[index++]=float(ntrklets);
 369   xnt[index++]=float(ntracks);
 370   xnt[index++]=float(nITS5or6);
 371   xnt[index++]=float(nTPCin);
 372   xnt[index++]=float(nTPCinEta09);
 373
 374   xnt[index++]=spd0cls;
 375
 376   xnt[index++]=xpile;
 377   xnt[index++]=expile;
 378   xnt[index++]=ypile;
 379   xnt[index++]=eypile;
 380   xnt[index++]=zpile;
 381   xnt[index++]=ezpile;
 382   xnt[index++]=(Float_t)ntrkspile;
 383
 384
 385
 386   // add recalculated vertices TRK and TPC
 387
 388   if(fRecoVtxTPC) {
 389     AliESDVertex *tpcvnc = ReconstructPrimaryVertexTPC(kFALSE);
 390     xnt[index++]=tpcvnc->GetXv();
 391     xnt[index++]=tpcvnc->GetXRes();
 392     xnt[index++]=tpcvnc->GetYv();
 393     xnt[index++]=tpcvnc->GetYRes();
 394     xnt[index++]=tpcvnc->GetZv();
 395     xnt[index++]=tpcvnc->GetZRes();
 396     xnt[index++]=tpcvnc->GetNContributors();
 397     delete tpcvnc; tpcvnc=0;
 398
 399     AliESDVertex *tpcvc = ReconstructPrimaryVertexTPC(kTRUE);
 400     xnt[index++]=tpcvc->GetXv();
 401     xnt[index++]=tpcvc->GetXRes();
 402     xnt[index++]=tpcvc->GetYv();
 403     xnt[index++]=tpcvc->GetYRes();
 404     xnt[index++]=tpcvc->GetZv();
 405     xnt[index++]=tpcvc->GetZRes();
 406     xnt[index++]=tpcvc->GetNContributors();
 407     delete tpcvc; tpcvc=0;
 408   } else index+=14;
 409
 410   if(fRecoVtxITSTPC) {
 411     AliESDVertex *trkvnc = ReconstructPrimaryVertexITSTPC(kFALSE);
 412     xnt[index++]=trkvnc->GetXv();
 413     xnt[index++]=trkvnc->GetXRes();
 414     xnt[index++]=trkvnc->GetYv();
 415     xnt[index++]=trkvnc->GetYRes();
 416     xnt[index++]=trkvnc->GetZv();
 417     xnt[index++]=trkvnc->GetZRes();
 418     xnt[index++]=trkvnc->GetNContributors();
 419     delete trkvnc; trkvnc=0;
 420
 421     AliESDVertex *trkvc = ReconstructPrimaryVertexITSTPC(kTRUE);
 422     xnt[index++]=trkvc->GetXv();
 423     xnt[index++]=trkvc->GetXRes();
 424     xnt[index++]=trkvc->GetYv();
 425     xnt[index++]=trkvc->GetYRes();
 426     xnt[index++]=trkvc->GetZv();
 427     xnt[index++]=trkvc->GetZRes();
 428     xnt[index++]=trkvc->GetNContributors();
 429     delete trkvc; trkvc=0;
 430   } else index+=14;
 431
 432   if(fRecoVtxITSTPCHalfEvent) {
 433     AliESDVertex *trkvncodd  = ReconstructPrimaryVertexITSTPC(kFALSE,1);
 434     AliESDVertex *trkvnceven = ReconstructPrimaryVertexITSTPC(kFALSE,2);
 435     if(trkvncodd->GetNContributors()>0 &&
 436        trkvnceven->GetNContributors()>0) {
 437       xnt[index++]=trkvncodd->GetXv()-trkvnceven->GetXv();
 438       xnt[index++]=trkvncodd->GetYv()-trkvnceven->GetYv();
 439       xnt[index++]=trkvncodd->GetZv()-trkvnceven->GetZv();
 440       xnt[index++]=TMath::Sqrt(trkvncodd->GetXRes()*trkvncodd->GetXRes()+trkvnceven->GetXRes()*trkvnceven->GetXRes());
 441       xnt[index++]=TMath::Sqrt(trkvncodd->GetYRes()*trkvncodd->GetYRes()+trkvnceven->GetYRes()*trkvnceven->GetYRes());
 442       xnt[index++]=TMath::Sqrt(trkvncodd->GetZRes()*trkvncodd->GetZRes()+trkvnceven->GetZRes()*trkvnceven->GetZRes());
 443       xnt[index++]=trkvnceven->GetNContributors();
 444       xnt[index++]=trkvncodd->GetNContributors();
 445     } else {
 446       xnt[index++]=0.;
 447       xnt[index++]=0.;
 448       xnt[index++]=0.;
 449       xnt[index++]=0.;
 450       xnt[index++]=0.;
 451       xnt[index++]=0.;
 452       xnt[index++]=-1;
 453       xnt[index++]=-1;
 454     }
 455     delete trkvncodd; trkvncodd=0;
 456     delete trkvnceven; trkvnceven=0;
 457   } else index+=8;
 458
 459
 460   if(index>isize) printf("AliAnalysisTaskVertexESD: ERROR, index!=isize\n");
 461   if(fFillNtuple) fNtupleVertexESD->Fill(xnt);
 462
 463   // Post the data already here
 464   PostData(1, fOutput);
 465
 466   return;
 467 }
 468
 469 //________________________________________________________________________
 470 void AliAnalysisTaskVertexESD::Terminate(Option_t *)
 471 {
 472   // Draw result to the screen
 473   // Called once at the end of the query
 474   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 475   if (!fOutput) {
 476     Printf("ERROR: fOutput not available");
 477     return;
 478   }
 479
 480   fNtupleVertexESD = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleVertexESD"));
 481
 482   return;
 483 }
 484
 485 //_________________________________________________________________________
 486 AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexTPC(Bool_t constr) const {
 487   // On the fly reco of TPC vertex from ESD
 488   AliESDEvent* evt = (AliESDEvent*) fInputEvent;
 489   AliVertexerTracks vertexer(evt->GetMagneticField());
 490   vertexer.SetTPCMode(); // defaults
 491   Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy);
 492   Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0};
 493   Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
 494   AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1);
 495   vertexer.SetVtxStart(initVertex);
 496   delete initVertex;
 497   if(!constr) vertexer.SetConstraintOff();
 498
 499   return vertexer.FindPrimaryVertex(evt);
 500 }
 501
 502 //_________________________________________________________________________
 503 AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexITSTPC(Bool_t constr,Int_t mode) const {
 504   // On the fly reco of ITS+TPC vertex from ESD
 505   // mode=0 use all tracks
 506   // mode=1 use odd-number tracks
 507   // mode=2 use even-number tracks
 508
 509   AliESDEvent* evt = (AliESDEvent*) fInputEvent;
 510   AliVertexerTracks vertexer(evt->GetMagneticField());
 511   vertexer.SetITSMode(); // defaults
 512   vertexer.SetMinClusters(4); // default is 5
 513   Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy);
 514   Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0};
 515   Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
 516   AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1);
 517   vertexer.SetVtxStart(initVertex);
 518   delete initVertex;
 519   if(!constr) vertexer.SetConstraintOff();
 520
 521   // use only ITS-TPC or only ITS-SA tracks
 522   if(fOnlyITSTPCTracks || fOnlyITSSATracks || mode>0) {
 523     Int_t iskip=0;
 524     Int_t *skip = new Int_t[evt->GetNumberOfTracks()];
 525     for(Int_t itr=0;itr<evt->GetNumberOfTracks(); itr++) {
 526       AliESDtrack* track = evt->GetTrack(itr);
 527       if(fOnlyITSTPCTracks && track->GetNcls(1)==0) { // skip ITSSA
 528         skip[iskip++]=itr;
 529         continue;
 530       }
 531       if(fOnlyITSSATracks && track->GetNcls(1)>0) { // skip ITSTPC
 532         skip[iskip++]=itr;
 533         continue;
 534       }
 535       if(mode==1 && itr%2==0) skip[iskip++]=itr;
 536       if(mode==2 && itr%2==1) skip[iskip++]=itr;
 537     }
 538     vertexer.SetSkipTracks(iskip,skip);
 539     delete [] skip; skip=NULL;
 540   }
 541
 542   return vertexer.FindPrimaryVertex(evt);
 543 }