64be43d34b21639faa85af7714f2384c7ce48e6c
[u/mrichter/AliRoot.git] / PWG1 / 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
38 #include "AliAnalysisTask.h"
39 #include "AliAnalysisManager.h"
40
41 #include "AliMultiplicity.h"
42 #include "AliVertexerTracks.h"
43 #include "AliESDtrack.h"
44 #include "AliExternalTrackParam.h"
45 #include "AliESDVertex.h"
46 #include "AliESDEvent.h"
47 #include "AliESDInputHandler.h"
48 #include "AliTrackReference.h"
49 //#include "AliTriggerAnalysis.h"
50
51 #include "AliMCEventHandler.h"
52 #include "AliMCEvent.h"
53 #include "AliStack.h"
54 #include "AliLog.h"
55
56 #include "AliGenEventHeader.h" 
57 #include "AliAnalysisTaskVertexESD.h"
58
59
60 ClassImp(AliAnalysisTaskVertexESD)
61
62 //________________________________________________________________________
63 AliAnalysisTaskVertexESD::AliAnalysisTaskVertexESD(const char *name) : 
64 AliAnalysisTaskSE(name), 
65 fCheckEventType(kTRUE),
66 fReadMC(kFALSE),
67 fRecoVtxTPC(kFALSE),
68 fRecoVtxITSTPC(kFALSE),
69 fRecoVtxITSTPCHalfEvent(kFALSE),
70 fOnlyITSTPCTracks(kFALSE),
71 fOnlyITSSATracks(kFALSE),
72 fFillNtuple(kFALSE),
73 fFillTreeBeamSpot(kFALSE),
74 fESD(0), 
75 fOutput(0), 
76 fNtupleVertexESD(0),
77 fhSPDVertexX(0),
78 fhSPDVertexY(0),
79 fhSPDVertexZ(0),
80 fhTRKVertexX(0),
81 fhTRKVertexY(0),
82 fhTRKVertexZ(0),
83 fhTPCVertexX(0),
84 fhTPCVertexY(0),
85 fhTPCVertexZ(0),
86 fhTrackRefs(0),
87 fTreeBeamSpot(0)
88 {
89   // Constructor
90
91   // Define input and output slots here
92   // Output slot #0 writes into a TList container
93   DefineOutput(1, TList::Class());  //My private output
94 }
95 //________________________________________________________________________
96 AliAnalysisTaskVertexESD::~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
109
110 //________________________________________________________________________
111 void AliAnalysisTaskVertexESD::UserCreateOutputObjects()
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
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");
121
122   fOutput->Add(fNtupleVertexESD);
123
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
143   fhTrackRefs = new TH2F("fhTrackRefs","Track references; x; y",1000,-4,4,1000,-4,4);
144   fOutput->Add(fhTrackRefs);
145
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);
160
161   PostData(1, fOutput);
162
163   return;
164 }
165
166 //________________________________________________________________________
167 void AliAnalysisTaskVertexESD::UserExec(Option_t *) 
168 {
169   // Main loop
170   // Called for each event
171   
172   if (!InputEvent()) {
173     Printf("ERROR: fESD not available");
174     return;
175   }
176   AliESDEvent* esdE = (AliESDEvent*) InputEvent();
177   
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   }
183
184
185   TArrayF mcVertex(3);
186   mcVertex[0]=9999.; mcVertex[1]=9999.; mcVertex[2]=9999.;
187   Float_t dNchdy=-999.;
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
246     
247   // Trigger
248   //ULong64_t triggerMask;
249   //ULong64_t spdFO = (1 << 14);
250   //ULong64_t v0left = (1 << 10);
251   //ULong64_t v0right = (1 << 11);
252   
253   //triggerMask=esdE->GetTriggerMask();
254   // MB1: SPDFO || V0L || V0R
255   //Bool_t eventTriggered = (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right))); 
256   //MB2: GFO && V0R
257   //triggerMask & spdFO && ((triggerMask&v0left) || (triggerMask&v0right))
258   //Bool_t eventTriggered = (triggerMask & spdFO); 
259  
260   //static AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis();
261   Bool_t eventTriggered = 0;//triggerAnalysis->IsTriggerFired(esdE, AliTriggerAnalysis::kSPDGFO /*| AliTriggerAnalysis::kOfflineFlag*/); 
262
263   // use response of AliPhysicsSelection
264   eventTriggered = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
265
266
267   Int_t nInFile = esdE->GetEventNumberInFile();
268
269  const AliMultiplicity *alimult = esdE->GetMultiplicity();
270   Int_t ntrklets=0,spd0cls=0;
271   if(alimult) {
272     ntrklets = alimult->GetNumberOfTracklets();
273
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
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
295   Double_t tstamp = esdE->GetTimeStamp();
296   Float_t cetTime =(tstamp-1262304000.)+7200.;
297
298   Float_t cetTime1h =(tstamp-1262304000.)+3600.;
299
300   cetTimeLHC = (Float_t)cetTime1h;
301
302   Int_t ntracks = esdE->GetNumberOfTracks();
303   Int_t nITS5or6=0,nTPCin=0,nTPCinEta09=0;
304   //printf("Tracks # = %d\n",esdE->GetNumberOfTracks());
305   for(Int_t itr=0; itr<ntracks; itr++) {
306     AliESDtrack *t = esdE->GetTrack(itr);
307     if(t->GetNcls(0)>=5) nITS5or6++;
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) {
310       nTPCin++;
311       if(TMath::Abs(t->GetTgl())<1.5) nTPCinEta09++;
312     }
313   }
314
315
316   const AliESDVertex *spdv=esdE->GetPrimaryVertexSPD();
317   const AliESDVertex *spdvp=esdE->GetPileupVertexSPD(0);
318   const AliESDVertex *tpcv=esdE->GetPrimaryVertexTPC();
319   const AliESDVertex *trkv=esdE->GetPrimaryVertexTracks();
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   } 
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   }
367
368   // fill ntuple
369   Int_t isize=82;
370   Float_t xnt[82]; for(Int_t iii=0;iii<isize;iii++) xnt[iii]=0.;
371
372   Int_t isizeSecNt=9;
373   Float_t secnt[9]; for(Int_t iii=0;iii<isizeSecNt;iii++) secnt[iii]=0.;
374
375   Int_t index=0;
376   Int_t indexSecNt=0;
377
378   xnt[index++]=(Float_t)esdE->GetRunNumber();
379   secnt[indexSecNt++]=(Float_t)esdE->GetRunNumber();
380   run = (Int_t)esdE->GetRunNumber();
381
382   xnt[index++]=cetTime; //(Float_t)esdE->GetTimeStamp();
383   //secnt[indexSecNt++]=cetTime;
384   secnt[indexSecNt++]=cetTime1h;
385
386   xnt[index++]=(Float_t)esdE->GetBunchCrossNumber();
387   secnt[indexSecNt++]=(Float_t)esdE->GetBunchCrossNumber();
388   bx = (Int_t)esdE->GetBunchCrossNumber();
389
390   xnt[index++]=(eventTriggered ? 1. : 0.);
391   secnt[indexSecNt++]=(eventTriggered ? 1. : 0.);
392   triggered = (UShort_t)(eventTriggered ? 1 : 0);
393
394   xnt[index++]=(Float_t)dNchdy;
395   
396   xnt[index++]=mcVertex[0];
397   xnt[index++]=mcVertex[1];
398   xnt[index++]=mcVertex[2];
399   
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();
407   TString spdtitle = spdv->GetTitle();
408   xnt[index++]=(spdtitle.Contains("vertexer: 3D") ? 1. : 0.);
409   xnt[index++]=spdv->GetDispersion();
410   
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();
418   TString tpctitle = tpcv->GetTitle();
419   xnt[index++]=(tpctitle.Contains("WithConstraint") ? 1. : 0.);
420   
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;
428   TString trktitle = trkv->GetTitle();
429   xnt[index++]=(trktitle.Contains("WithConstraint") ? 1. : 0.);  
430
431   xnt[index++]=float(ntrklets);
432   secnt[indexSecNt++]=float(ntrklets);
433   ntrkletsS = (UShort_t)ntrklets;
434
435   xnt[index++]=float(ntracks);
436   xnt[index++]=float(nITS5or6);
437
438   xnt[index++]=float(nTPCin);
439   xnt[index++]=float(nTPCinEta09);
440
441   xnt[index++]=spd0cls;
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  
453   // add recalculated vertices TRK and TPC
454
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;
475   } else index+=14;
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();
486   
487     secnt[indexSecNt++]=trkvnc->GetXv();
488     secnt[indexSecNt++]=trkvnc->GetYv();
489     secnt[indexSecNt++]=trkvnc->GetZv();
490     secnt[indexSecNt++]=trkvnc->GetNContributors();
491
492     xTRKnc = (Float_t)trkvnc->GetXv();
493     yTRKnc = (Float_t)trkvnc->GetYv();
494     zTRKnc = (Float_t)trkvnc->GetZv();
495     ntrksTRKnc = (trkvnc->GetNContributors()<0 ? 0 : (UShort_t)trkvnc->GetNContributors());
496
497     delete trkvnc; trkvnc=0;
498
499
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;
509   } else index+=14;
510   
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();
517       xnt[index++]=trkvncodd->GetYv()-trkvnceven->GetYv();
518       xnt[index++]=trkvncodd->GetZv()-trkvnceven->GetZv();
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());
521       xnt[index++]=TMath::Sqrt(trkvncodd->GetZRes()*trkvncodd->GetZRes()+trkvnceven->GetZRes()*trkvnceven->GetZRes());
522       xnt[index++]=trkvnceven->GetNContributors();
523       xnt[index++]=trkvncodd->GetNContributors();
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;
532       xnt[index++]=-1;
533     }
534     delete trkvncodd; trkvncodd=0;
535     delete trkvnceven; trkvnceven=0;    
536   } else index+=8;
537   
538
539   if(index>isize) printf("AliAnalysisTaskVertexESD: ERROR, index!=isize\n");
540   if(fFillNtuple) fNtupleVertexESD->Fill(xnt);
541
542   if(indexSecNt>isizeSecNt) printf("AliAnalysisTaskVertexESD: ERROR, indexSecNt!=isizeSecNt\n");
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
548   
549   // Post the data already here
550   PostData(1, fOutput);
551
552   return;
553 }      
554
555 //________________________________________________________________________
556 void AliAnalysisTaskVertexESD::Terminate(Option_t *) 
557 {
558   // Draw result to the screen
559   // Called once at the end of the query
560   fOutput = dynamic_cast<TList*> (GetOutputData(1));
561   if (!fOutput) {     
562     Printf("ERROR: fOutput not available");
563     return;
564   }
565   
566   if (!fNtupleVertexESD){
567     Printf("ERROR: fNtuple not available");
568     return;
569   }
570   
571   fNtupleVertexESD = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleVertexESD"));
572
573
574   return;
575 }
576
577 //_________________________________________________________________________
578 AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexTPC(Bool_t constr) const {
579   // On the fly reco of TPC vertex from ESD
580   AliESDEvent* evt = (AliESDEvent*) fInputEvent;
581   AliVertexerTracks vertexer(evt->GetMagneticField());
582   vertexer.SetTPCMode(); // defaults
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);
587   vertexer.SetVtxStart(initVertex);
588   delete initVertex;
589   if(!constr) vertexer.SetConstraintOff();
590
591   return vertexer.FindPrimaryVertex(evt);
592 }
593
594 //_________________________________________________________________________
595 AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexITSTPC(Bool_t constr,Int_t mode) const {
596   // On the fly reco of ITS+TPC vertex from ESD
597   // mode=0 use all tracks
598   // mode=1 use odd-number tracks
599   // mode=2 use even-number tracks
600
601   AliESDEvent* evt = (AliESDEvent*) fInputEvent;
602   AliVertexerTracks vertexer(evt->GetMagneticField());
603   vertexer.SetITSMode(); // defaults
604   vertexer.SetMinClusters(4); // default is 5
605   //vertexer.SetITSpureSA();
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);
610   vertexer.SetVtxStart(initVertex);
611   delete initVertex;
612   if(!constr) vertexer.SetConstraintOff();
613
614   // use only ITS-TPC or only ITS-SA tracks
615   if(fOnlyITSTPCTracks || fOnlyITSSATracks || mode>0) {
616     Int_t iskip=0;
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);
620       if(fOnlyITSTPCTracks && track->GetNcls(1)==0) { // skip ITSSA
621         skip[iskip++]=itr;
622         continue;
623       }
624       if(fOnlyITSSATracks && track->GetNcls(1)>0) { // skip ITSTPC
625         skip[iskip++]=itr;
626         continue;
627       }
628       if(mode==1 && itr%2==0) skip[iskip++]=itr;
629       if(mode==2 && itr%2==1) skip[iskip++]=itr;
630     }
631     vertexer.SetSkipTracks(iskip,skip);
632     delete [] skip; skip=NULL;
633   }
634
635   return vertexer.FindPrimaryVertex(evt);
636 }