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