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