773355b1e2544230ed009eec61197db2137b4485
[u/mrichter/AliRoot.git] / PWGPP / global / AliAnalysisTaskVertexESD.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //*************************************************************************
17 // Class AliAnalysisTaskVertexESD
18 // AliAnalysisTask to extract from ESD the information for the analysis
19 // of the primary vertex reconstruction efficiency and resolution
20 // (for MC events) and distributions (for real data). Three vertices:
21 // - SPD tracklets
22 // - ITS+TPC tracks
23 // - TPC-only tracks
24 //
25 // Author: A.Dainese, andrea.dainese@pd.infn.it
26 //*************************************************************************
27
28 #include <TChain.h>
29 #include <TTree.h>
30 #include <TNtuple.h>
31 #include <TBranch.h>
32 #include <TClonesArray.h>
33 #include <TObjArray.h>
34 #include <TH1F.h>
35 #include <TH2F.h>  
36 #include <TCanvas.h>
37 #include <TGraphAsymmErrors.h>
38 #include <TFile.h>
39
40 #include "AliAnalysisTask.h"
41 #include "AliAnalysisManager.h"
42
43 #include "AliMultiplicity.h"
44 #include "AliVertexerTracks.h"
45 #include "AliESDtrack.h"
46 #include "AliExternalTrackParam.h"
47 #include "AliESDVertex.h"
48 #include "AliVEvent.h"
49 #include "AliESDInputHandler.h"
50 #include "AliTrackReference.h"
51 //#include "AliTriggerAnalysis.h"
52
53 #include "AliMCEventHandler.h"
54 #include "AliMCEvent.h"
55 #include "AliStack.h"
56 #include "AliLog.h"
57
58 #include "AliGenEventHeader.h" 
59 #include "AliAnalysisTaskVertexESD.h"
60
61
62 ClassImp(AliAnalysisTaskVertexESD)
63
64 //________________________________________________________________________
65 AliAnalysisTaskVertexESD::AliAnalysisTaskVertexESD(const char *name) : 
66 AliAnalysisTaskSE(name), 
67   fCheckEventType(kTRUE),
68   fReadMC(kFALSE),
69   fRecoVtxTPC(kTRUE),
70   fRecoVtxITSTPC(kTRUE),
71   fRecoVtxITSTPCHalfEvent(kFALSE),
72   fOnlyITSTPCTracks(kFALSE),
73   fOnlyITSSATracks(kFALSE),
74   fFillNtuple(kFALSE),
75   fFillTreeBeamSpot(kFALSE),
76   fESD(0), 
77   fOutput(0), 
78   fNtupleVertexESD(0),
79   fhSPDVertexX(0),
80   fhSPDVertexY(0),
81   fhSPDVertexZ(0),
82   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 //________________________________________________________________________
123 AliAnalysisTaskVertexESD::~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 //________________________________________________________________________
138 void 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 //________________________________________________________________________
257 void 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 //________________________________________________________________________
737 void 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 //_________________________________________________________________________
852 AliESDVertex* 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 //_________________________________________________________________________
875 AliESDVertex* 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 }