[u/mrichter/AliRoot.git] / PWG1 / global / AliAnalysisTaskVertexESD.cxx

/**************************************************************************
 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//*************************************************************************
// Class AliAnalysisTaskVertexESD
// AliAnalysisTask to extract from ESD the information for the analysis
// of the primary vertex reconstruction efficiency and resolution
// (for MC events) and distributions (for real data). Three vertices:
// - SPD tracklets
// - ITS+TPC tracks
// - TPC-only tracks
//
// Author: A.Dainese, andrea.dainese@pd.infn.it
//*************************************************************************

#include <TChain.h>
#include <TTree.h>
#include <TNtuple.h>
#include <TBranch.h>
#include <TClonesArray.h>
#include <TObjArray.h>
#include <TH1F.h>
#include <TH2F.h>  
#include <TCanvas.h>

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"

#include "AliMultiplicity.h"
#include "AliVertexerTracks.h"
#include "AliESDtrack.h"
#include "AliExternalTrackParam.h"
#include "AliESDVertex.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliTrackReference.h"
//#include "AliTriggerAnalysis.h"

#include "AliMCEventHandler.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliLog.h"

#include "AliGenEventHeader.h" 
#include "AliAnalysisTaskVertexESD.h"


ClassImp(AliAnalysisTaskVertexESD)

//________________________________________________________________________
AliAnalysisTaskVertexESD::AliAnalysisTaskVertexESD(const char *name) : 
AliAnalysisTaskSE(name), 
fCheckEventType(kTRUE),
fReadMC(kFALSE),
fRecoVtxTPC(kFALSE),
fRecoVtxITSTPC(kFALSE),
fRecoVtxITSTPCHalfEvent(kFALSE),
fOnlyITSTPCTracks(kFALSE),
fOnlyITSSATracks(kFALSE),
fFillNtuple(kFALSE),
fESD(0), 
fOutput(0), 
fNtupleVertexESD(0),
fhSPDVertexX(0),
fhSPDVertexY(0),
fhSPDVertexZ(0),
fhTRKVertexX(0),
fhTRKVertexY(0),
fhTRKVertexZ(0),
fhTPCVertexX(0),
fhTPCVertexY(0),
fhTPCVertexZ(0),
fhTrackRefs(0)
{
  // Constructor

  // Define input and output slots here
  // Output slot #0 writes into a TList container
  DefineOutput(1, TList::Class());  //My private output
}
//________________________________________________________________________
AliAnalysisTaskVertexESD::~AliAnalysisTaskVertexESD()
{
  // Destructor

  // histograms are in the output list and deleted when the output
  // list is deleted by the TSelector dtor

  if (fOutput) {
    delete fOutput;
    fOutput = 0;
  }
}


//________________________________________________________________________
void AliAnalysisTaskVertexESD::UserCreateOutputObjects()
{
  // Create histograms
  // Called once

  // Several histograms are more conveniently managed in a TList
  fOutput = new TList;
  fOutput->SetOwner();

  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:avntrksTRKnc");

  fOutput->Add(fNtupleVertexESD);

  fhSPDVertexX = new TH1F("fhSPDVertexX","SPDVertex x; x vertex [cm]; events",200,-1,1);
  fOutput->Add(fhSPDVertexX);
  fhSPDVertexY = new TH1F("fhSPDVertexY","SPDVertex y; y vertex [cm]; events",200,-1,1);
  fOutput->Add(fhSPDVertexY);
  fhSPDVertexZ = new TH1F("fhSPDVertexZ","SPDVertex z; z vertex [cm]; events",200,-20,20);
  fOutput->Add(fhSPDVertexZ);
  fhTRKVertexX = new TH1F("fhTRKVertexX","TRKVertex x; x vertex [cm]; events",200,-1,1);
  fOutput->Add(fhTRKVertexX);
  fhTRKVertexY = new TH1F("fhTRKVertexY","TRKVertex y; y vertex [cm]; events",200,-1,1);
  fOutput->Add(fhTRKVertexY);
  fhTRKVertexZ = new TH1F("fhTRKVertexZ","TRKVertex z; z vertex [cm]; events",200,-20,20);
  fOutput->Add(fhTRKVertexZ);
  fhTPCVertexX = new TH1F("fhTPCVertexX","TPCVertex x; x vertex [cm]; events",200,-3,3);
  fOutput->Add(fhTPCVertexX);
  fhTPCVertexY = new TH1F("fhTPCVertexY","TPCVertex y; y vertex [cm]; events",200,-3,3);
  fOutput->Add(fhTPCVertexY);
  fhTPCVertexZ = new TH1F("fhTPCVertexZ","TPCVertex z; z vertex [cm]; events",200,-20,20);
  fOutput->Add(fhTPCVertexZ);

  fhTrackRefs = new TH2F("fhTrackRefs","Track references; x; y",1000,-4,4,1000,-4,4);
  fOutput->Add(fhTrackRefs);

  return;
}

//________________________________________________________________________
void AliAnalysisTaskVertexESD::UserExec(Option_t *) 
{
  // Main loop
  // Called for each event
  
  if (!InputEvent()) {
    Printf("ERROR: fESD not available");
    return;
  }
  AliESDEvent* esdE = (AliESDEvent*) InputEvent();
  
  if(fCheckEventType && (esdE->GetEventType())!=7) return; 


  TArrayF mcVertex(3);
  mcVertex[0]=9999.; mcVertex[1]=9999.; mcVertex[2]=9999.;
  Float_t dNchdy=-999.;
  // ***********  MC info ***************
  if (fReadMC) {
    AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
    if (!eventHandler) {
      Printf("ERROR: Could not retrieve MC event handler");
      return;
    }
    
    AliMCEvent* mcEvent = eventHandler->MCEvent();
    if (!mcEvent) {
      Printf("ERROR: Could not retrieve MC event");
      return;
    }
    
    AliStack* stack = mcEvent->Stack();
    if (!stack) {
      AliDebug(AliLog::kError, "Stack not available");
      return;
    }
    
    AliHeader* header = mcEvent->Header();
    if (!header) {
      AliDebug(AliLog::kError, "Header not available");
      return;
    }
    AliGenEventHeader* genHeader = header->GenEventHeader();
    genHeader->PrimaryVertex(mcVertex);

    Int_t ngenpart = (Int_t)stack->GetNtrack();
    //printf("# generated particles = %d\n",ngenpart);
    dNchdy=0;
    for(Int_t ip=0; ip<ngenpart; ip++) {
      TParticle* part = (TParticle*)stack->Particle(ip);
      // keep only electorns, muons, pions, kaons and protons
      Int_t apdg = TMath::Abs(part->GetPdgCode());
      if(apdg!=11 && apdg!=13 && apdg!=211 && apdg!=321 && apdg!=2212) continue;      
      // reject secondaries
      if(TMath::Sqrt((part->Vx()-mcVertex[0])*(part->Vx()-mcVertex[0])+(part->Vy()-mcVertex[1])*(part->Vy()-mcVertex[1]))>0.0010) continue;
      // reject incoming protons
      Double_t energy  = part->Energy();
      if(energy>900.) continue;
      Double_t pz = part->Pz();
      Double_t y = 0.5*TMath::Log((energy+pz+1.e-13)/(energy-pz+1.e-13));
      if(TMath::Abs(y)<1.0) dNchdy += 0.5; // count 1/2 of particles in |y|<1
      // tracks refs
      TClonesArray *trefs=0;
      Int_t ntrefs = mcEvent->GetParticleAndTR(ip,part,trefs);
      if(ntrefs<0) continue;
      for(Int_t iref=0; iref<ntrefs; iref++) {
	AliTrackReference *tref = (AliTrackReference*)trefs->At(iref);
	if(tref->R()>10.) continue;
	fhTrackRefs->Fill(tref->X(),tref->Y());
      }
    }
    //printf("# primary particles = %7.1f\n",dNchdy);
  } 
  // ***********  MC info ***************

    
  // Trigger
  //ULong64_t triggerMask;
  //ULong64_t spdFO = (1 << 14);
  //ULong64_t v0left = (1 << 10);
  //ULong64_t v0right = (1 << 11);
  
  //triggerMask=esdE->GetTriggerMask();
  // MB1: SPDFO || V0L || V0R
  //Bool_t eventTriggered = (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right))); 
  //MB2: GFO && V0R
  //triggerMask & spdFO && ((triggerMask&v0left) || (triggerMask&v0right))
  //Bool_t eventTriggered = (triggerMask & spdFO); 
 
  //static AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis();
  Bool_t eventTriggered = 0;//triggerAnalysis->IsTriggerFired(esdE, AliTriggerAnalysis::kSPDGFO /*| AliTriggerAnalysis::kOfflineFlag*/); 

  // use response of AliPhysicsSelection
  eventTriggered = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();

  Int_t ntracks = esdE->GetNumberOfTracks();
  Int_t nITS5or6=0,nTPCin=0,nTPCinEta09=0;
  //printf("Tracks # = %d\n",esdE->GetNumberOfTracks());
  for(Int_t itr=0; itr<ntracks; itr++) {
    AliESDtrack *t = esdE->GetTrack(itr);
    if(t->GetNcls(0)>=5) nITS5or6++;
    Double_t z0; t->GetZAt(0,esdE->GetMagneticField(),z0);
    if(t->GetNcls(1)>0 && TMath::Abs(t->GetD(0,0,esdE->GetMagneticField()))<2.8 && TMath::Abs(z0)<20) {
      nTPCin++;
      if(TMath::Abs(t->GetTgl())<1.5) nTPCinEta09++;
    }
  }

    
  const AliESDVertex *spdv=esdE->GetPrimaryVertexSPD();
  const AliESDVertex *spdvp=esdE->GetPileupVertexSPD(0);
  const AliESDVertex *tpcv=esdE->GetPrimaryVertexTPC();
  const AliESDVertex *trkv=esdE->GetPrimaryVertexTracks();

  const AliMultiplicity *alimult = esdE->GetMultiplicity();
  Int_t ntrklets=0,spd0cls=0;
  if(alimult) {
    ntrklets = alimult->GetNumberOfTracklets();
    for(Int_t l=0;l<alimult->GetNumberOfTracklets();l++){
      if(alimult->GetDeltaPhi(l)<-9998.) ntrklets--;
    }
    spd0cls = alimult->GetNumberOfSingleClusters()+ntrklets;
  }
  
  // fill histos
  
  if(spdv) {
    if(spdv->GetNContributors()>0) {
      TString title=spdv->GetTitle();
      if(title.Contains("3D")) {
	fhSPDVertexX->Fill(spdv->GetXv());
	fhSPDVertexY->Fill(spdv->GetYv());
      }
      fhSPDVertexZ->Fill(spdv->GetZv());
    }
  }
  
  if(trkv) {
    if(trkv->GetNContributors()>0) {
      fhTRKVertexX->Fill(trkv->GetXv());
      fhTRKVertexY->Fill(trkv->GetYv());
      fhTRKVertexZ->Fill(trkv->GetZv());
    }
  }
  
  if(tpcv) {
    if(tpcv->GetNContributors()>0) {
      fhTPCVertexX->Fill(tpcv->GetXv());
      fhTPCVertexY->Fill(tpcv->GetYv());
      fhTPCVertexZ->Fill(tpcv->GetZv());
    }
  } 

  Float_t xpile=-999.;
  Float_t ypile=-999.;
  Float_t zpile=-999.;
  Float_t expile=-999.;
  Float_t eypile=-999.;
  Float_t ezpile=-999.;
  Int_t ntrkspile=-1;

  if(esdE->GetNumberOfPileupVerticesSPD()>0 && spdvp){
    xpile=spdvp->GetXv();
    expile=spdvp->GetXRes();
    ypile=spdvp->GetYv();
    eypile=spdvp->GetYRes();
    zpile=spdvp->GetZv();
    ezpile=spdvp->GetZRes();
    ntrkspile=spdvp->GetNContributors();  
  }

  // fill ntuple
  Int_t isize=81;
  Float_t xnt[81]; for(Int_t iii=0;iii<isize;iii++) xnt[iii]=0.;
  
  Int_t index=0;

  xnt[index++]=(Float_t)esdE->GetRunNumber();
  xnt[index++]=(Float_t)esdE->GetTimeStamp();
  xnt[index++]=(Float_t)esdE->GetBunchCrossNumber();
  xnt[index++]=(eventTriggered ? 1. : 0.);

  xnt[index++]=(Float_t)dNchdy;


  xnt[index++]=mcVertex[0];
  xnt[index++]=mcVertex[1];
  xnt[index++]=mcVertex[2];
  
  xnt[index++]=spdv->GetXv();
  xnt[index++]=spdv->GetXRes();
  xnt[index++]=spdv->GetYv();
  xnt[index++]=spdv->GetYRes();
  xnt[index++]=spdv->GetZv();
  xnt[index++]=spdv->GetZRes();
  xnt[index++]=spdv->GetNContributors();
  TString spdtitle = spdv->GetTitle();
  xnt[index++]=(spdtitle.Contains("vertexer: 3D") ? 1. : 0.);
  xnt[index++]=spdv->GetDispersion();
  
  xnt[index++]=tpcv->GetXv();
  xnt[index++]=tpcv->GetXRes();
  xnt[index++]=tpcv->GetYv();
  xnt[index++]=tpcv->GetYRes();
  xnt[index++]=tpcv->GetZv();
  xnt[index++]=tpcv->GetZRes();
  xnt[index++]=tpcv->GetNContributors();
  TString tpctitle = tpcv->GetTitle();
  xnt[index++]=(tpctitle.Contains("WithConstraint") ? 1. : 0.);
  
  xnt[index++]=trkv->GetXv();
  xnt[index++]=trkv->GetXRes();
  xnt[index++]=trkv->GetYv();
  xnt[index++]=trkv->GetYRes();
  xnt[index++]=trkv->GetZv();
  xnt[index++]=trkv->GetZRes();
  xnt[index++]=trkv->GetNContributors();// tpccontrorig;
  TString trktitle = trkv->GetTitle();
  xnt[index++]=(trktitle.Contains("WithConstraint") ? 1. : 0.);  

  xnt[index++]=float(ntrklets);
  xnt[index++]=float(ntracks);
  xnt[index++]=float(nITS5or6);
  xnt[index++]=float(nTPCin);
  xnt[index++]=float(nTPCinEta09);

  xnt[index++]=spd0cls;
    
  xnt[index++]=xpile;
  xnt[index++]=expile;
  xnt[index++]=ypile;
  xnt[index++]=eypile;
  xnt[index++]=zpile;
  xnt[index++]=ezpile;
  xnt[index++]=(Float_t)ntrkspile;  
    
    
  // add recalculated vertices TRK and TPC

  if(fRecoVtxTPC) {
    AliESDVertex *tpcvnc = ReconstructPrimaryVertexTPC(kFALSE);
    xnt[index++]=tpcvnc->GetXv();
    xnt[index++]=tpcvnc->GetXRes();
    xnt[index++]=tpcvnc->GetYv();
    xnt[index++]=tpcvnc->GetYRes();
    xnt[index++]=tpcvnc->GetZv();
    xnt[index++]=tpcvnc->GetZRes();
    xnt[index++]=tpcvnc->GetNContributors();
    delete tpcvnc; tpcvnc=0;
    
    AliESDVertex *tpcvc = ReconstructPrimaryVertexTPC(kTRUE);
    xnt[index++]=tpcvc->GetXv();
    xnt[index++]=tpcvc->GetXRes();
    xnt[index++]=tpcvc->GetYv();
    xnt[index++]=tpcvc->GetYRes();
    xnt[index++]=tpcvc->GetZv();
    xnt[index++]=tpcvc->GetZRes();
    xnt[index++]=tpcvc->GetNContributors();
    delete tpcvc; tpcvc=0;
  } else index+=14;

  if(fRecoVtxITSTPC) {
    AliESDVertex *trkvnc = ReconstructPrimaryVertexITSTPC(kFALSE);
    xnt[index++]=trkvnc->GetXv();
    xnt[index++]=trkvnc->GetXRes();
    xnt[index++]=trkvnc->GetYv();
    xnt[index++]=trkvnc->GetYRes();
    xnt[index++]=trkvnc->GetZv();
    xnt[index++]=trkvnc->GetZRes();
    xnt[index++]=trkvnc->GetNContributors();
    delete trkvnc; trkvnc=0;
    
    AliESDVertex *trkvc = ReconstructPrimaryVertexITSTPC(kTRUE);
    xnt[index++]=trkvc->GetXv();
    xnt[index++]=trkvc->GetXRes();
    xnt[index++]=trkvc->GetYv();
    xnt[index++]=trkvc->GetYRes();
    xnt[index++]=trkvc->GetZv();
    xnt[index++]=trkvc->GetZRes();
    xnt[index++]=trkvc->GetNContributors();
    delete trkvc; trkvc=0;
  } else index+=14;
  
  if(fRecoVtxITSTPCHalfEvent) {
    AliESDVertex *trkvncodd  = ReconstructPrimaryVertexITSTPC(kFALSE,1);
    AliESDVertex *trkvnceven = ReconstructPrimaryVertexITSTPC(kFALSE,2);
    if(trkvncodd->GetNContributors()>0 && 
       trkvnceven->GetNContributors()>0) {
      xnt[index++]=trkvncodd->GetXv()-trkvnceven->GetXv();
      xnt[index++]=trkvncodd->GetYv()-trkvnceven->GetYv();
      xnt[index++]=trkvncodd->GetZv()-trkvnceven->GetZv();
      xnt[index++]=TMath::Sqrt(trkvncodd->GetXRes()*trkvncodd->GetXRes()+trkvnceven->GetXRes()*trkvnceven->GetXRes());
      xnt[index++]=TMath::Sqrt(trkvncodd->GetYRes()*trkvncodd->GetYRes()+trkvnceven->GetYRes()*trkvnceven->GetYRes());
      xnt[index++]=TMath::Sqrt(trkvncodd->GetZRes()*trkvncodd->GetZRes()+trkvnceven->GetZRes()*trkvnceven->GetZRes());
      xnt[index++]=0.5*(trkvncodd->GetNContributors()+trkvnceven->GetNContributors());
    } else {
      xnt[index++]=0.;
      xnt[index++]=0.;
      xnt[index++]=0.;
      xnt[index++]=0.;
      xnt[index++]=0.;
      xnt[index++]=0.;
      xnt[index++]=-1;
    }
    delete trkvncodd; trkvncodd=0;
    delete trkvnceven; trkvnceven=0;    
  } else index+=7;
  

  if(index>isize) printf("AliAnalysisTaskVertexESD: ERROR, index!=isize\n");
  if(fFillNtuple) fNtupleVertexESD->Fill(xnt);
  
  // Post the data already here
  PostData(1, fOutput);

  return;
}      

//________________________________________________________________________
void AliAnalysisTaskVertexESD::Terminate(Option_t *) 
{
  // Draw result to the screen
  // Called once at the end of the query
  fOutput = dynamic_cast<TList*> (GetOutputData(1));
  if (!fOutput) {     
    Printf("ERROR: fOutput not available");
    return;
  }

  fNtupleVertexESD = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleVertexESD"));

  return;
}

//_________________________________________________________________________
AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexTPC(Bool_t constr) const {
  // On the fly reco of TPC vertex from ESD
  AliESDEvent* evt = (AliESDEvent*) fInputEvent;
  AliVertexerTracks vertexer(evt->GetMagneticField());
  vertexer.SetTPCMode(); // defaults
  Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy);
  Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0}; 
  Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
  AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1);
  vertexer.SetVtxStart(initVertex);
  delete initVertex;
  if(!constr) vertexer.SetConstraintOff();

  return vertexer.FindPrimaryVertex(evt);
}

//_________________________________________________________________________
AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexITSTPC(Bool_t constr,Int_t mode) const {
  // On the fly reco of ITS+TPC vertex from ESD
  // mode=0 use all tracks
  // mode=1 use odd-number tracks
  // mode=2 use even-number tracks

  AliESDEvent* evt = (AliESDEvent*) fInputEvent;
  AliVertexerTracks vertexer(evt->GetMagneticField());
  vertexer.SetITSMode(); // defaults
  vertexer.SetMinClusters(4); // default is 5
  Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy);
  Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0}; 
  Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
  AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1);
  vertexer.SetVtxStart(initVertex);
  delete initVertex;
  if(!constr) vertexer.SetConstraintOff();

  // use only ITS-TPC or only ITS-SA tracks
  if(fOnlyITSTPCTracks || fOnlyITSSATracks || mode>0) {
    Int_t iskip=0;
    Int_t *skip = new Int_t[evt->GetNumberOfTracks()];
    for(Int_t itr=0;itr<evt->GetNumberOfTracks(); itr++) {
      AliESDtrack* track = evt->GetTrack(itr);
      if(fOnlyITSTPCTracks && track->GetNcls(1)==0) { // skip ITSSA
	skip[iskip++]=itr;
	continue;
      }
      if(fOnlyITSSATracks && track->GetNcls(1)>0) { // skip ITSTPC
	skip[iskip++]=itr;
	continue;
      }
      if(mode==1 && itr%2==0) skip[iskip++]=itr;
      if(mode==2 && itr%2==1) skip[iskip++]=itr;
    }
    vertexer.SetSkipTracks(iskip,skip);
    delete [] skip; skip=NULL;
  }

  return vertexer.FindPrimaryVertex(evt);
}
Commit	Line	Data
388ca814	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"
388ca814	47	#include "AliESDInputHandler.h"
388ca814	48	#include "AliTrackReference.h"
38d12819	49	//#include "AliTriggerAnalysis.h"
388ca814	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) :
74917c1f	64	AliAnalysisTaskSE(name),
4148ee5b	65	fCheckEventType(kTRUE),
388ca814	66	fReadMC(kFALSE),
94f15f9b	67	fRecoVtxTPC(kFALSE),
94f15f9b	68	fRecoVtxITSTPC(kFALSE),
4d95fd6c	69	fRecoVtxITSTPCHalfEvent(kFALSE),
cb393b39	70	fOnlyITSTPCTracks(kFALSE),
cb393b39	71	fOnlyITSSATracks(kFALSE),
b588553c	72	fFillNtuple(kFALSE),
388ca814	73	fESD(0),
388ca814	74	fOutput(0),
388ca814	75	fNtupleVertexESD(0),
b588553c	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),
38d12819	85	fhTrackRefs(0)
388ca814	86	{
	87	// Constructor
	88
	89	// Define input and output slots here
388ca814	90	// Output slot #0 writes into a TList container
41deb8e3	91	DefineOutput(1, TList::Class()); //My private output
388ca814	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
388ca814	107
388ca814	108	//________________________________________________________________________
74917c1f	109	void AliAnalysisTaskVertexESD::UserCreateOutputObjects()
388ca814	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
4d95fd6c	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:avntrksTRKnc");
388ca814	119
	120	fOutput->Add(fNtupleVertexESD);
	121
b588553c	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
388ca814	141	fhTrackRefs = new TH2F("fhTrackRefs","Track references; x; y",1000,-4,4,1000,-4,4);
	142	fOutput->Add(fhTrackRefs);
	143
f6232e6e	144	return;
388ca814	145	}
	146
	147	//________________________________________________________________________
74917c1f	148	void AliAnalysisTaskVertexESD::UserExec(Option_t *)
388ca814	149	{
	150	// Main loop
	151	// Called for each event
	152
74917c1f	153	if (!InputEvent()) {
388ca814	154	Printf("ERROR: fESD not available");
	155	return;
	156	}
74917c1f	157	AliESDEvent* esdE = (AliESDEvent*) InputEvent();
	158
	159	if(fCheckEventType && (esdE->GetEventType())!=7) return;
38d12819	160
38d12819	161
388ca814	162	TArrayF mcVertex(3);
	163	mcVertex[0]=9999.; mcVertex[1]=9999.; mcVertex[2]=9999.;
	164	Float_t dNchdy=-999.;
388ca814	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
388ca814	223
388ca814	224	// Trigger
38d12819	225	//ULong64_t triggerMask;
	226	//ULong64_t spdFO = (1 << 14);
	227	//ULong64_t v0left = (1 << 10);
	228	//ULong64_t v0right = (1 << 11);
388ca814	229
74917c1f	230	//triggerMask=esdE->GetTriggerMask();
388ca814	231	// MB1: SPDFO \|\| V0L \|\| V0R
38d12819	232	//Bool_t eventTriggered = (triggerMask & spdFO \|\| ((triggerMask & v0left) \|\| (triggerMask & v0right)));
388ca814	233	//MB2: GFO && V0R
388ca814	234	//triggerMask & spdFO && ((triggerMask&v0left) \|\| (triggerMask&v0right))
38d12819	235	//Bool_t eventTriggered = (triggerMask & spdFO);
388ca814	236
38d12819	237	//static AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis();
74917c1f	238	Bool_t eventTriggered = 0;//triggerAnalysis->IsTriggerFired(esdE, AliTriggerAnalysis::kSPDGFO /\| AliTriggerAnalysis::kOfflineFlag/);
388ca814	239
8c7764e0	240	// use response of AliPhysicsSelection
	241	eventTriggered = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
	242
74917c1f	243	Int_t ntracks = esdE->GetNumberOfTracks();
388ca814	244	Int_t nITS5or6=0,nTPCin=0,nTPCinEta09=0;
74917c1f	245	//printf("Tracks # = %d\n",esdE->GetNumberOfTracks());
388ca814	246	for(Int_t itr=0; itr<ntracks; itr++) {
74917c1f	247	AliESDtrack *t = esdE->GetTrack(itr);
388ca814	248	if(t->GetNcls(0)>=5) nITS5or6++;
74917c1f	249	Double_t z0; t->GetZAt(0,esdE->GetMagneticField(),z0);
74917c1f	250	if(t->GetNcls(1)>0 && TMath::Abs(t->GetD(0,0,esdE->GetMagneticField()))<2.8 && TMath::Abs(z0)<20) {
388ca814	251	nTPCin++;
	252	if(TMath::Abs(t->GetTgl())<1.5) nTPCinEta09++;
	253	}
388ca814	254	}
	255
	256
74917c1f	257	const AliESDVertex *spdv=esdE->GetPrimaryVertexSPD();
39258194	258	const AliESDVertex *spdvp=esdE->GetPileupVertexSPD(0);
74917c1f	259	const AliESDVertex *tpcv=esdE->GetPrimaryVertexTPC();
74917c1f	260	const AliESDVertex *trkv=esdE->GetPrimaryVertexTracks();
388ca814	261
74917c1f	262	const AliMultiplicity *alimult = esdE->GetMultiplicity();
388ca814	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	}
b588553c	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	}
39258194	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	}
388ca814	318
b588553c	319	// fill ntuple
4d95fd6c	320	Int_t isize=81;
4d95fd6c	321	Float_t xnt[81]; for(Int_t iii=0;iii<isize;iii++) xnt[iii]=0.;
388ca814	322
97805f7d	323	Int_t index=0;
97805f7d	324
74917c1f	325	xnt[index++]=(Float_t)esdE->GetRunNumber();
74917c1f	326	xnt[index++]=(Float_t)esdE->GetTimeStamp();
08c95603	327	xnt[index++]=(Float_t)esdE->GetBunchCrossNumber();
8c7764e0	328	xnt[index++]=(eventTriggered ? 1. : 0.);
	329
	330	xnt[index++]=(Float_t)dNchdy;
	331
97805f7d	332
	333	xnt[index++]=mcVertex[0];
	334	xnt[index++]=mcVertex[1];
	335	xnt[index++]=mcVertex[2];
388ca814	336
97805f7d	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();
8c7764e0	344	TString spdtitle = spdv->GetTitle();
	345	xnt[index++]=(spdtitle.Contains("vertexer: 3D") ? 1. : 0.);
	346	xnt[index++]=spdv->GetDispersion();
388ca814	347
97805f7d	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();
8c7764e0	355	TString tpctitle = tpcv->GetTitle();
8c7764e0	356	xnt[index++]=(tpctitle.Contains("WithConstraint") ? 1. : 0.);
388ca814	357
97805f7d	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;
8c7764e0	365	TString trktitle = trkv->GetTitle();
8c7764e0	366	xnt[index++]=(trktitle.Contains("WithConstraint") ? 1. : 0.);
388ca814	367
97805f7d	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);
388ca814	373
97805f7d	374	xnt[index++]=spd0cls;
39258194	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
8c7764e0	386	// add recalculated vertices TRK and TPC
97805f7d	387
8c7764e0	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;
39258194	408	} else index+=14;
8c7764e0	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;
39258194	430	} else index+=14;
39258194	431
4d95fd6c	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();
4d95fd6c	438	xnt[index++]=trkvncodd->GetYv()-trkvnceven->GetYv();
4d95fd6c	439	xnt[index++]=trkvncodd->GetZv()-trkvnceven->GetZv();
007ec307	440	xnt[index++]=TMath::Sqrt(trkvncodd->GetXRes()trkvncodd->GetXRes()+trkvnceven->GetXRes()trkvnceven->GetXRes());
007ec307	441	xnt[index++]=TMath::Sqrt(trkvncodd->GetYRes()trkvncodd->GetYRes()+trkvnceven->GetYRes()trkvnceven->GetYRes());
4d95fd6c	442	xnt[index++]=TMath::Sqrt(trkvncodd->GetZRes()trkvncodd->GetZRes()+trkvnceven->GetZRes()trkvnceven->GetZRes());
	443	xnt[index++]=0.5*(trkvncodd->GetNContributors()+trkvnceven->GetNContributors());
	444	} else {
	445	xnt[index++]=0.;
	446	xnt[index++]=0.;
	447	xnt[index++]=0.;
	448	xnt[index++]=0.;
	449	xnt[index++]=0.;
	450	xnt[index++]=0.;
	451	xnt[index++]=-1;
	452	}
	453	delete trkvncodd; trkvncodd=0;
	454	delete trkvnceven; trkvnceven=0;
	455	} else index+=7;
	456
	457
8c7764e0	458	if(index>isize) printf("AliAnalysisTaskVertexESD: ERROR, index!=isize\n");
b588553c	459	if(fFillNtuple) fNtupleVertexESD->Fill(xnt);
f6232e6e	460
388ca814	461	// Post the data already here
41deb8e3	462	PostData(1, fOutput);
388ca814	463
	464	return;
	465	}
	466
	467	//________________________________________________________________________
	468	void AliAnalysisTaskVertexESD::Terminate(Option_t *)
	469	{
	470	// Draw result to the screen
	471	// Called once at the end of the query
2ac70243	472	fOutput = dynamic_cast<TList*> (GetOutputData(1));
388ca814	473	if (!fOutput) {
	474	Printf("ERROR: fOutput not available");
	475	return;
	476	}
	477
	478	fNtupleVertexESD = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleVertexESD"));
	479
	480	return;
	481	}
	482
	483	//_________________________________________________________________________
8c7764e0	484	AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexTPC(Bool_t constr) const {
388ca814	485	// On the fly reco of TPC vertex from ESD
74917c1f	486	AliESDEvent* evt = (AliESDEvent*) fInputEvent;
74917c1f	487	AliVertexerTracks vertexer(evt->GetMagneticField());
97805f7d	488	vertexer.SetTPCMode(); // defaults
8c7764e0	489	Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy);
	490	Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0};
	491	Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
	492	AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1);
388ca814	493	vertexer.SetVtxStart(initVertex);
388ca814	494	delete initVertex;
8c7764e0	495	if(!constr) vertexer.SetConstraintOff();
388ca814	496
74917c1f	497	return vertexer.FindPrimaryVertex(evt);
388ca814	498	}
94f15f9b	499
94f15f9b	500	//_________________________________________________________________________
4d95fd6c	501	AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexITSTPC(Bool_t constr,Int_t mode) const {
94f15f9b	502	// On the fly reco of ITS+TPC vertex from ESD
4d95fd6c	503	// mode=0 use all tracks
	504	// mode=1 use odd-number tracks
	505	// mode=2 use even-number tracks
	506
74917c1f	507	AliESDEvent* evt = (AliESDEvent*) fInputEvent;
74917c1f	508	AliVertexerTracks vertexer(evt->GetMagneticField());
94f15f9b	509	vertexer.SetITSMode(); // defaults
8c7764e0	510	vertexer.SetMinClusters(4); // default is 5
	511	Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy);
	512	Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0};
	513	Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
	514	AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1);
94f15f9b	515	vertexer.SetVtxStart(initVertex);
94f15f9b	516	delete initVertex;
8c7764e0	517	if(!constr) vertexer.SetConstraintOff();
94f15f9b	518
cb393b39	519	// use only ITS-TPC or only ITS-SA tracks
5457d525	520	if(fOnlyITSTPCTracks \|\| fOnlyITSSATracks \|\| mode>0) {
cb393b39	521	Int_t iskip=0;
74917c1f	522	Int_t *skip = new Int_t[evt->GetNumberOfTracks()];
	523	for(Int_t itr=0;itr<evt->GetNumberOfTracks(); itr++) {
	524	AliESDtrack* track = evt->GetTrack(itr);
cb393b39	525	if(fOnlyITSTPCTracks && track->GetNcls(1)==0) { // skip ITSSA
cb393b39	526	skip[iskip++]=itr;
4d95fd6c	527	continue;
cb393b39	528	}
	529	if(fOnlyITSSATracks && track->GetNcls(1)>0) { // skip ITSTPC
	530	skip[iskip++]=itr;
4d95fd6c	531	continue;
cb393b39	532	}
4d95fd6c	533	if(mode==1 && itr%2==0) skip[iskip++]=itr;
4d95fd6c	534	if(mode==2 && itr%2==1) skip[iskip++]=itr;
cb393b39	535	}
	536	vertexer.SetSkipTracks(iskip,skip);
	537	delete [] skip; skip=NULL;
	538	}
	539
74917c1f	540	return vertexer.FindPrimaryVertex(evt);
94f15f9b	541	}