[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 "AliESDfriend.h"
#include "AliESDInputHandler.h"
#include "AliTrackPointArray.h"
#include "AliTrackReference.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) : 
AliAnalysisTask(name, "VertexESDTask"), 
fReadMC(kFALSE),
fESD(0), 
fESDfriend(0), 
fOutput(0), 
fNtupleVertexESD(0),
fhTrackPoints(0),
fhTrackRefs(0)
{
  // Constructor

  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #0 writes into a TList container
  DefineOutput(0, 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::ConnectInputData(Option_t *) 
{
  // Connect ESD or AOD here
  // Called once

  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  if (!tree) {
    Printf("ERROR: Could not read chain from input slot 0");
  } else {
    // Disable all branches and enable only the needed ones
    // The next two lines are different when data produced as AliESDEvent is read
    //    tree->SetBranchStatus("*", kFALSE);
    tree->SetBranchStatus("fTriggerMask", 1);
    tree->SetBranchStatus("fSPDVertex*", 1);
    tree->SetBranchStatus("fSPDMult*", 1);
    
    tree->SetBranchStatus("ESDfriend*", 1);
    tree->SetBranchAddress("ESDfriend.",&fESDfriend);

    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
    
    if (!esdH) {
      Printf("ERROR: Could not get ESDInputHandler");
    } else fESD = esdH->GetEvent();
  }
  
  return;
}

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

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

  fNtupleVertexESD = new TNtuple("fNtupleVertexESD","vertices","xtrue:ytrue:ztrue:xSPD:xdiffSPD:xerrSPD:ySPD:ydiffSPD:yerrSPD:zSPD:zdiffSPD:zerrSPD:ntrksSPD:xTPC:xdiffTPC:xerrTPC:yTPC:ydiffTPC:yerrTPC:zTPC:zdiffTPC:zerrTPC:ntrksTPC:xTRK:xdiffTRK:xerrTRK:yTRK:ydiffTRK:yerrTRK:zTRK:zdiffTRK:zerrTRK:ntrksTRK:ntrklets:nESDtracks:nITSrefit5or6:nTPCin:nTPCinEta09:dndygen:triggered:SPD3D:SPD0cls:constrTRK");

  fOutput->Add(fNtupleVertexESD);

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

 return;
}

//________________________________________________________________________
void AliAnalysisTaskVertexESD::Exec(Option_t *) 
{
  // Main loop
  // Called for each event
  
  if (!fESD) {
    Printf("ERROR: fESD not available");
    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 ***************

  // ***********  ESD friends ***********
  fESD->SetESDfriend(fESDfriend); //Attach the friend to the ESD
  // ***********  ESD friends ***********

  //
    
  // Trigger
  ULong64_t triggerMask;
  ULong64_t spdFO = (1 << 14);
  ULong64_t v0left = (1 << 11);
  ULong64_t v0right = (1 << 12);
  
  triggerMask=fESD->GetTriggerMask();
  // MB1: SPDFO || V0L || V0R
  Bool_t eventTriggered = (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right))); 
  //MB2: GFO && V0R
  //triggerMask & spdFO && ((triggerMask&v0left) || (triggerMask&v0right))
 

  Int_t ntracks = fESD->GetNumberOfTracks();
  Int_t nITS5or6=0,nTPCin=0,nTPCinEta09=0;
  //printf("Tracks # = %d\n",fESD->GetNumberOfTracks());
  for(Int_t itr=0; itr<ntracks; itr++) {
    AliESDtrack *t = fESD->GetTrack(itr);
    if(t->GetNcls(0)>=5) nITS5or6++;
    Double_t z0; t->GetZAt(0,fESD->GetMagneticField(),z0);
    if(t->GetNcls(1)>0 && TMath::Abs(t->GetD(0,0,fESD->GetMagneticField()))<2.8 && TMath::Abs(z0)<20) {
      nTPCin++;
      if(TMath::Abs(t->GetTgl())<1.5) nTPCinEta09++;
    }
    // tracks with two SPD points
    if(!TESTBIT(t->GetITSClusterMap(),0) ||
       !TESTBIT(t->GetITSClusterMap(),1)) continue; 
    // AliTrackPoints
    const AliTrackPointArray *array = t->GetTrackPointArray();
    AliTrackPoint point;
    for(Int_t ipt=0; ipt<array->GetNPoints(); ipt++) {
      array->GetPoint(point,ipt);
      fhTrackPoints->Fill(point.GetX(),point.GetY());
    }
  }

    
  const AliESDVertex *spdv=fESD->GetPrimaryVertexSPD();
  const AliESDVertex *tpcv=fESD->GetPrimaryVertexTPC();
  const AliESDVertex *trkv=fESD->GetPrimaryVertex();

  //Float_t tpccontrorig=tpcv->GetNContributors();
  //tpcv = 0;
  //tpcv = ReconstructPrimaryVertexTPC();

  const AliMultiplicity *alimult = fESD->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;
  }


  Float_t xnt[43];
  
  xnt[0]=mcVertex[0];
  xnt[1]=mcVertex[1];
  xnt[2]=mcVertex[2];
  
  xnt[3]=spdv->GetXv();
  xnt[4]=(spdv->GetXv()-mcVertex[0])*10000.;
  xnt[5]=spdv->GetXRes();
  xnt[6]=spdv->GetYv();
  xnt[7]=(spdv->GetYv()-mcVertex[1])*10000.;
  xnt[8]=spdv->GetYRes();
  xnt[9]=spdv->GetZv();
  xnt[10]=(spdv->GetZv()-mcVertex[2])*10000.;
  xnt[11]=spdv->GetZRes();
  xnt[12]=spdv->GetNContributors();
  
  xnt[13]=tpcv->GetXv();
  xnt[14]=(tpcv->GetXv()-mcVertex[0])*10000.;
  xnt[15]=tpcv->GetXRes();
  xnt[16]=tpcv->GetYv();
  xnt[17]=(tpcv->GetYv()-mcVertex[1])*10000.;
  xnt[18]=tpcv->GetYRes();
  xnt[19]=tpcv->GetZv();
  xnt[20]=(tpcv->GetZv()-mcVertex[2])*10000.;
  xnt[21]=tpcv->GetZRes();
  xnt[22]=tpcv->GetNContributors();
  
  xnt[23]=trkv->GetXv();
  xnt[24]=(trkv->GetXv()-mcVertex[0])*10000.;
  xnt[25]=trkv->GetXRes();
  xnt[26]=trkv->GetYv();
  xnt[27]=(trkv->GetYv()-mcVertex[1])*10000.;
  xnt[28]=trkv->GetYRes();
  xnt[29]=trkv->GetZv();
  xnt[30]=(trkv->GetZv()-mcVertex[2])*10000.;
  xnt[31]=trkv->GetZRes();
  xnt[32]=trkv->GetNContributors();// tpccontrorig;
  

  xnt[33]=float(ntrklets);
  xnt[34]=float(ntracks);
  xnt[35]=float(nITS5or6);
  xnt[36]=float(nTPCin);
  xnt[37]=float(nTPCinEta09);

  xnt[38]=float(dNchdy);

  xnt[39]=0.; 
  if(eventTriggered) xnt[39]=1.;

  xnt[40]=0.; 
  TString spdtitle = spdv->GetTitle();
  if(spdtitle.Contains("vertexer: 3D")) xnt[40]=1.;

  xnt[42]=0.; 
  TString trktitle = trkv->GetTitle();
  if(trktitle.Contains("WithConstraint")) xnt[42]=1.;

  xnt[41]=spd0cls;

  fNtupleVertexESD->Fill(xnt);
    
  // Post the data already here
  PostData(0, 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(0));
  if (!fOutput) {     
    Printf("ERROR: fOutput not available");
    return;
  }

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

  return;
}

//_________________________________________________________________________
AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexTPC() {
  // On the fly reco of TPC vertex from ESD

  AliVertexerTracks vertexer(fESD->GetMagneticField());
  vertexer.SetTPCMode(0.1,1.0,5.,0,1,3.,0.1,1.5); // defaults
  //vertexer.SetTPCMode(0.1,1.0,5.,0,1,3.,0.1,1.5); // defaults
  Double_t pos[3]={+0.0220,-0.0340,+0.270}; 
  Double_t err[3]={0.0200,0.0200,7.5};
  AliESDVertex *initVertex = new AliESDVertex(pos,err);
  vertexer.SetVtxStart(initVertex);
  delete initVertex;
  //vertexer.SetConstraintOff();

  return vertexer.FindPrimaryVertex(fESD);
}
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"
	47	#include "AliESDfriend.h"
	48	#include "AliESDInputHandler.h"
	49	#include "AliTrackPointArray.h"
	50	#include "AliTrackReference.h"
	51
	52	#include "AliMCEventHandler.h"
	53	#include "AliMCEvent.h"
	54	#include "AliStack.h"
	55	#include "AliLog.h"
	56
	57	#include "AliGenEventHeader.h"
	58	#include "AliAnalysisTaskVertexESD.h"
	59
	60
	61	ClassImp(AliAnalysisTaskVertexESD)
	62
	63	//________________________________________________________________________
	64	AliAnalysisTaskVertexESD::AliAnalysisTaskVertexESD(const char *name) :
65	AliAnalysisTask(name, "VertexESDTask"),
66	fReadMC(kFALSE),
67	fESD(0),
68	fESDfriend(0),
69	fOutput(0),
70	fNtupleVertexESD(0),
71	fhTrackPoints(0),
72	fhTrackRefs(0)
73	{
74	// Constructor
75
76	// Define input and output slots here
77	// Input slot #0 works with a TChain
78	DefineInput(0, TChain::Class());
79	// Output slot #0 writes into a TList container
80	DefineOutput(0, TList::Class()); //My private output
81	}
82	//________________________________________________________________________
83	AliAnalysisTaskVertexESD::~AliAnalysisTaskVertexESD()
84	{
85	// Destructor
86
87	// histograms are in the output list and deleted when the output
88	// list is deleted by the TSelector dtor
89
90	if (fOutput) {
91	delete fOutput;
92	fOutput = 0;
93	}
94	}
95
96	//________________________________________________________________________
97	void AliAnalysisTaskVertexESD::ConnectInputData(Option_t *)
98	{
99	// Connect ESD or AOD here
100	// Called once
101
102	TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
103	if (!tree) {
104	Printf("ERROR: Could not read chain from input slot 0");
105	} else {
106	// Disable all branches and enable only the needed ones
107	// The next two lines are different when data produced as AliESDEvent is read
108	// tree->SetBranchStatus("*", kFALSE);
109	tree->SetBranchStatus("fTriggerMask", 1);
110	tree->SetBranchStatus("fSPDVertex*", 1);
111	tree->SetBranchStatus("fSPDMult*", 1);
112
113	tree->SetBranchStatus("ESDfriend*", 1);
114	tree->SetBranchAddress("ESDfriend.",&fESDfriend);
115
116	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
117
118	if (!esdH) {
119	Printf("ERROR: Could not get ESDInputHandler");
120	} else fESD = esdH->GetEvent();
121	}
122
123	return;
124	}
125
126	//________________________________________________________________________
127	void AliAnalysisTaskVertexESD::CreateOutputObjects()
128	{
129	// Create histograms
130	// Called once
131
132	// Several histograms are more conveniently managed in a TList
133	fOutput = new TList;
134	fOutput->SetOwner();
135
136	fNtupleVertexESD = new TNtuple("fNtupleVertexESD","vertices","xtrue:ytrue:ztrue:xSPD:xdiffSPD:xerrSPD:ySPD:ydiffSPD:yerrSPD:zSPD:zdiffSPD:zerrSPD:ntrksSPD:xTPC:xdiffTPC:xerrTPC:yTPC:ydiffTPC:yerrTPC:zTPC:zdiffTPC:zerrTPC:ntrksTPC:xTRK:xdiffTRK:xerrTRK:yTRK:ydiffTRK:yerrTRK:zTRK:zdiffTRK:zerrTRK:ntrksTRK:ntrklets:nESDtracks:nITSrefit5or6:nTPCin:nTPCinEta09:dndygen:triggered:SPD3D:SPD0cls:constrTRK");
137
138	fOutput->Add(fNtupleVertexESD);
139
140	fhTrackPoints = new TH2F("fhTrackPoints","Track points; x; y",1000,-4,4,1000,-4,4);
141	fOutput->Add(fhTrackPoints);
142	fhTrackRefs = new TH2F("fhTrackRefs","Track references; x; y",1000,-4,4,1000,-4,4);
143	fOutput->Add(fhTrackRefs);
144
145	return;
146	}
147
148	//________________________________________________________________________
149	void AliAnalysisTaskVertexESD::Exec(Option_t *)
150	{
151	// Main loop
152	// Called for each event
153
154	if (!fESD) {
155	Printf("ERROR: fESD not available");
156	return;
157	}
158
159
160	TArrayF mcVertex(3);
161	mcVertex[0]=9999.; mcVertex[1]=9999.; mcVertex[2]=9999.;
162	Float_t dNchdy=-999.;
163
164	// ********* MC info *************
165	if (fReadMC) {
166	AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
167	if (!eventHandler) {
168	Printf("ERROR: Could not retrieve MC event handler");
169	return;
170	}
171
172	AliMCEvent* mcEvent = eventHandler->MCEvent();
173	if (!mcEvent) {
174	Printf("ERROR: Could not retrieve MC event");
175	return;
176	}
177
178	AliStack* stack = mcEvent->Stack();
179	if (!stack) {
180	AliDebug(AliLog::kError, "Stack not available");
181	return;
182	}
183
184	AliHeader* header = mcEvent->Header();
185	if (!header) {
186	AliDebug(AliLog::kError, "Header not available");
187	return;
188	}
189	AliGenEventHeader* genHeader = header->GenEventHeader();
190	genHeader->PrimaryVertex(mcVertex);
191
192	Int_t ngenpart = (Int_t)stack->GetNtrack();
193	//printf("# generated particles = %d\n",ngenpart);
194	dNchdy=0;
195	for(Int_t ip=0; ip<ngenpart; ip++) {
196	TParticle* part = (TParticle*)stack->Particle(ip);
197	// keep only electorns, muons, pions, kaons and protons
198	Int_t apdg = TMath::Abs(part->GetPdgCode());
199	if(apdg!=11 && apdg!=13 && apdg!=211 && apdg!=321 && apdg!=2212) continue;
200	// reject secondaries
201	if(TMath::Sqrt((part->Vx()-mcVertex[0])(part->Vx()-mcVertex[0])+(part->Vy()-mcVertex[1])(part->Vy()-mcVertex[1]))>0.0010) continue;
202	// reject incoming protons
203	Double_t energy = part->Energy();
204	if(energy>900.) continue;
205	Double_t pz = part->Pz();
206	Double_t y = 0.5*TMath::Log((energy+pz+1.e-13)/(energy-pz+1.e-13));
207	if(TMath::Abs(y)<1.0) dNchdy += 0.5; // count 1/2 of particles in \|y\|<1
208	// tracks refs
209	TClonesArray *trefs=0;
210	Int_t ntrefs = mcEvent->GetParticleAndTR(ip,part,trefs);
211	if(ntrefs<0) continue;
212	for(Int_t iref=0; iref<ntrefs; iref++) {
213	AliTrackReference tref = (AliTrackReference)trefs->At(iref);
214	if(tref->R()>10.) continue;
215	fhTrackRefs->Fill(tref->X(),tref->Y());
216	}
217	}
218	//printf("# primary particles = %7.1f\n",dNchdy);
219	}
220	// ********* MC info *************
221
222	// ********* ESD friends *********
223	fESD->SetESDfriend(fESDfriend); //Attach the friend to the ESD
224	// ********* ESD friends *********
225
226	//
227
228	// Trigger
229	ULong64_t triggerMask;
230	ULong64_t spdFO = (1 << 14);
231	ULong64_t v0left = (1 << 11);
232	ULong64_t v0right = (1 << 12);
233
234	triggerMask=fESD->GetTriggerMask();
235	// MB1: SPDFO \|\| V0L \|\| V0R
236	Bool_t eventTriggered = (triggerMask & spdFO \|\| ((triggerMask & v0left) \|\| (triggerMask & v0right)));
237	//MB2: GFO && V0R
238	//triggerMask & spdFO && ((triggerMask&v0left) \|\| (triggerMask&v0right))
239
240
241	Int_t ntracks = fESD->GetNumberOfTracks();
242	Int_t nITS5or6=0,nTPCin=0,nTPCinEta09=0;
243	//printf("Tracks # = %d\n",fESD->GetNumberOfTracks());
244	for(Int_t itr=0; itr<ntracks; itr++) {
245	AliESDtrack *t = fESD->GetTrack(itr);
246	if(t->GetNcls(0)>=5) nITS5or6++;
247	Double_t z0; t->GetZAt(0,fESD->GetMagneticField(),z0);
248	if(t->GetNcls(1)>0 && TMath::Abs(t->GetD(0,0,fESD->GetMagneticField()))<2.8 && TMath::Abs(z0)<20) {
249	nTPCin++;
250	if(TMath::Abs(t->GetTgl())<1.5) nTPCinEta09++;
251	}
252	// tracks with two SPD points
253	if(!TESTBIT(t->GetITSClusterMap(),0) \|\|
254	!TESTBIT(t->GetITSClusterMap(),1)) continue;
255	// AliTrackPoints
256	const AliTrackPointArray *array = t->GetTrackPointArray();
257	AliTrackPoint point;
258	for(Int_t ipt=0; ipt<array->GetNPoints(); ipt++) {
259	array->GetPoint(point,ipt);
260	fhTrackPoints->Fill(point.GetX(),point.GetY());
261	}
262	}
263
264
265	const AliESDVertex *spdv=fESD->GetPrimaryVertexSPD();
266	const AliESDVertex *tpcv=fESD->GetPrimaryVertexTPC();
267	const AliESDVertex *trkv=fESD->GetPrimaryVertex();
268
269	//Float_t tpccontrorig=tpcv->GetNContributors();
270	//tpcv = 0;
271	//tpcv = ReconstructPrimaryVertexTPC();
272
273	const AliMultiplicity *alimult = fESD->GetMultiplicity();
274	Int_t ntrklets=0,spd0cls=0;
275	if(alimult) {
276	ntrklets = alimult->GetNumberOfTracklets();
277	for(Int_t l=0;l<alimult->GetNumberOfTracklets();l++){
278	if(alimult->GetDeltaPhi(l)<-9998.) ntrklets--;
279	}
280	spd0cls = alimult->GetNumberOfSingleClusters()+ntrklets;
281	}
282
283
284	Float_t xnt[43];
285
286	xnt[0]=mcVertex[0];
287	xnt[1]=mcVertex[1];
288	xnt[2]=mcVertex[2];
289
290	xnt[3]=spdv->GetXv();
291	xnt[4]=(spdv->GetXv()-mcVertex[0])*10000.;
292	xnt[5]=spdv->GetXRes();
293	xnt[6]=spdv->GetYv();
294	xnt[7]=(spdv->GetYv()-mcVertex[1])*10000.;
295	xnt[8]=spdv->GetYRes();
296	xnt[9]=spdv->GetZv();
297	xnt[10]=(spdv->GetZv()-mcVertex[2])*10000.;
298	xnt[11]=spdv->GetZRes();
299	xnt[12]=spdv->GetNContributors();
300
301	xnt[13]=tpcv->GetXv();
302	xnt[14]=(tpcv->GetXv()-mcVertex[0])*10000.;
303	xnt[15]=tpcv->GetXRes();
304	xnt[16]=tpcv->GetYv();
305	xnt[17]=(tpcv->GetYv()-mcVertex[1])*10000.;
306	xnt[18]=tpcv->GetYRes();
307	xnt[19]=tpcv->GetZv();
308	xnt[20]=(tpcv->GetZv()-mcVertex[2])*10000.;
309	xnt[21]=tpcv->GetZRes();
310	xnt[22]=tpcv->GetNContributors();
311
312	xnt[23]=trkv->GetXv();
313	xnt[24]=(trkv->GetXv()-mcVertex[0])*10000.;
314	xnt[25]=trkv->GetXRes();
315	xnt[26]=trkv->GetYv();
316	xnt[27]=(trkv->GetYv()-mcVertex[1])*10000.;
317	xnt[28]=trkv->GetYRes();
318	xnt[29]=trkv->GetZv();
319	xnt[30]=(trkv->GetZv()-mcVertex[2])*10000.;
320	xnt[31]=trkv->GetZRes();
321	xnt[32]=trkv->GetNContributors();// tpccontrorig;
322
323
324	xnt[33]=float(ntrklets);
325	xnt[34]=float(ntracks);
326	xnt[35]=float(nITS5or6);
327	xnt[36]=float(nTPCin);
328	xnt[37]=float(nTPCinEta09);
329
330	xnt[38]=float(dNchdy);
331
332	xnt[39]=0.;
333	if(eventTriggered) xnt[39]=1.;
334
335	xnt[40]=0.;
336	TString spdtitle = spdv->GetTitle();
337	if(spdtitle.Contains("vertexer: 3D")) xnt[40]=1.;
338
339	xnt[42]=0.;
340	TString trktitle = trkv->GetTitle();
341	if(trktitle.Contains("WithConstraint")) xnt[42]=1.;
342
343	xnt[41]=spd0cls;
344
345	fNtupleVertexESD->Fill(xnt);
346
347	// Post the data already here
348	PostData(0, fOutput);
349
350	return;
351	}
352
353	//________________________________________________________________________
354	void AliAnalysisTaskVertexESD::Terminate(Option_t *)
355	{
356	// Draw result to the screen
357	// Called once at the end of the query
358	fOutput = dynamic_cast<TList*> (GetOutputData(0));
359	if (!fOutput) {
360	Printf("ERROR: fOutput not available");
361	return;
362	}
363
364	fNtupleVertexESD = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleVertexESD"));
365
366	return;
367	}
368
369	//_________________________________________________________________________
370	AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexTPC() {
371	// On the fly reco of TPC vertex from ESD
372
373	AliVertexerTracks vertexer(fESD->GetMagneticField());
374	vertexer.SetTPCMode(0.1,1.0,5.,0,1,3.,0.1,1.5); // defaults
375	//vertexer.SetTPCMode(0.1,1.0,5.,0,1,3.,0.1,1.5); // defaults
376	Double_t pos[3]={+0.0220,-0.0340,+0.270};
377	Double_t err[3]={0.0200,0.0200,7.5};
378	AliESDVertex *initVertex = new AliESDVertex(pos,err);
379	vertexer.SetVtxStart(initVertex);
380	delete initVertex;
381	//vertexer.SetConstraintOff();
382
383	return vertexer.FindPrimaryVertex(fESD);
384	}