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