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