Added bunch crossing number to ntuple
[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 fOnlyITSTPCTracks(kFALSE),
70 fOnlyITSSATracks(kFALSE),
71 fFillNtuple(kFALSE),
72 fESD(0), 
73 fOutput(0), 
74 fNtupleVertexESD(0),
75 fhSPDVertexX(0),
76 fhSPDVertexY(0),
77 fhSPDVertexZ(0),
78 fhTRKVertexX(0),
79 fhTRKVertexY(0),
80 fhTRKVertexZ(0),
81 fhTPCVertexX(0),
82 fhTPCVertexY(0),
83 fhTPCVertexZ(0),
84 fhTrackRefs(0)
85 {
86   // Constructor
87
88   // Define input and output slots here
89   // Output slot #0 writes into a TList container
90   DefineOutput(1, TList::Class());  //My private output
91 }
92 //________________________________________________________________________
93 AliAnalysisTaskVertexESD::~AliAnalysisTaskVertexESD()
94 {
95   // Destructor
96
97   // histograms are in the output list and deleted when the output
98   // list is deleted by the TSelector dtor
99
100   if (fOutput) {
101     delete fOutput;
102     fOutput = 0;
103   }
104 }
105
106
107 //________________________________________________________________________
108 void AliAnalysisTaskVertexESD::UserCreateOutputObjects()
109 {
110   // Create histograms
111   // Called once
112
113   // Several histograms are more conveniently managed in a TList
114   fOutput = new TList;
115   fOutput->SetOwner();
116
117   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:xTPCnc:xerrTPCnc:yTPCnc:yerrTPCnc:zTPCnc:zerrTPCnc:ntrksTPCnc:xTRKnc:xerrTRKnc:yTRKnc:yerrTRKnc:zTRKnc:zerrTRKnc:ntrksTRKnc:xTPCc:xerrTPCc:yTPCc:yerrTPCc:zTPCc:zerrTPCc:ntrksTPCc:xTRKc:xerrTRKc:yTRKc:yerrTRKc:zTRKc:zerrTRKc:ntrksTRKc");
118
119   fOutput->Add(fNtupleVertexESD);
120
121   fhSPDVertexX = new TH1F("fhSPDVertexX","SPDVertex x; x vertex [cm]; events",200,-1,1);
122   fOutput->Add(fhSPDVertexX);
123   fhSPDVertexY = new TH1F("fhSPDVertexY","SPDVertex y; y vertex [cm]; events",200,-1,1);
124   fOutput->Add(fhSPDVertexY);
125   fhSPDVertexZ = new TH1F("fhSPDVertexZ","SPDVertex z; z vertex [cm]; events",200,-20,20);
126   fOutput->Add(fhSPDVertexZ);
127   fhTRKVertexX = new TH1F("fhTRKVertexX","TRKVertex x; x vertex [cm]; events",200,-1,1);
128   fOutput->Add(fhTRKVertexX);
129   fhTRKVertexY = new TH1F("fhTRKVertexY","TRKVertex y; y vertex [cm]; events",200,-1,1);
130   fOutput->Add(fhTRKVertexY);
131   fhTRKVertexZ = new TH1F("fhTRKVertexZ","TRKVertex z; z vertex [cm]; events",200,-20,20);
132   fOutput->Add(fhTRKVertexZ);
133   fhTPCVertexX = new TH1F("fhTPCVertexX","TPCVertex x; x vertex [cm]; events",200,-3,3);
134   fOutput->Add(fhTPCVertexX);
135   fhTPCVertexY = new TH1F("fhTPCVertexY","TPCVertex y; y vertex [cm]; events",200,-3,3);
136   fOutput->Add(fhTPCVertexY);
137   fhTPCVertexZ = new TH1F("fhTPCVertexZ","TPCVertex z; z vertex [cm]; events",200,-20,20);
138   fOutput->Add(fhTPCVertexZ);
139
140   fhTrackRefs = new TH2F("fhTrackRefs","Track references; x; y",1000,-4,4,1000,-4,4);
141   fOutput->Add(fhTrackRefs);
142
143   return;
144 }
145
146 //________________________________________________________________________
147 void AliAnalysisTaskVertexESD::UserExec(Option_t *) 
148 {
149   // Main loop
150   // Called for each event
151   
152   if (!InputEvent()) {
153     Printf("ERROR: fESD not available");
154     return;
155   }
156
157   AliESDEvent* esdE = (AliESDEvent*) InputEvent();
158   
159   if(fCheckEventType && (esdE->GetEventType())!=7) return; 
160
161
162   TArrayF mcVertex(3);
163   mcVertex[0]=9999.; mcVertex[1]=9999.; mcVertex[2]=9999.;
164   Float_t dNchdy=-999.;
165
166   // ***********  MC info ***************
167   if (fReadMC) {
168     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
169     if (!eventHandler) {
170       Printf("ERROR: Could not retrieve MC event handler");
171       return;
172     }
173     
174     AliMCEvent* mcEvent = eventHandler->MCEvent();
175     if (!mcEvent) {
176       Printf("ERROR: Could not retrieve MC event");
177       return;
178     }
179     
180     AliStack* stack = mcEvent->Stack();
181     if (!stack) {
182       AliDebug(AliLog::kError, "Stack not available");
183       return;
184     }
185     
186     AliHeader* header = mcEvent->Header();
187     if (!header) {
188       AliDebug(AliLog::kError, "Header not available");
189       return;
190     }
191     AliGenEventHeader* genHeader = header->GenEventHeader();
192     genHeader->PrimaryVertex(mcVertex);
193
194     Int_t ngenpart = (Int_t)stack->GetNtrack();
195     //printf("# generated particles = %d\n",ngenpart);
196     dNchdy=0;
197     for(Int_t ip=0; ip<ngenpart; ip++) {
198       TParticle* part = (TParticle*)stack->Particle(ip);
199       // keep only electorns, muons, pions, kaons and protons
200       Int_t apdg = TMath::Abs(part->GetPdgCode());
201       if(apdg!=11 && apdg!=13 && apdg!=211 && apdg!=321 && apdg!=2212) continue;      
202       // reject secondaries
203       if(TMath::Sqrt((part->Vx()-mcVertex[0])*(part->Vx()-mcVertex[0])+(part->Vy()-mcVertex[1])*(part->Vy()-mcVertex[1]))>0.0010) continue;
204       // reject incoming protons
205       Double_t energy  = part->Energy();
206       if(energy>900.) continue;
207       Double_t pz = part->Pz();
208       Double_t y = 0.5*TMath::Log((energy+pz+1.e-13)/(energy-pz+1.e-13));
209       if(TMath::Abs(y)<1.0) dNchdy += 0.5; // count 1/2 of particles in |y|<1
210       // tracks refs
211       TClonesArray *trefs=0;
212       Int_t ntrefs = mcEvent->GetParticleAndTR(ip,part,trefs);
213       if(ntrefs<0) continue;
214       for(Int_t iref=0; iref<ntrefs; iref++) {
215         AliTrackReference *tref = (AliTrackReference*)trefs->At(iref);
216         if(tref->R()>10.) continue;
217         fhTrackRefs->Fill(tref->X(),tref->Y());
218       }
219     }
220     //printf("# primary particles = %7.1f\n",dNchdy);
221   } 
222   // ***********  MC info ***************
223
224     
225   // Trigger
226   //ULong64_t triggerMask;
227   //ULong64_t spdFO = (1 << 14);
228   //ULong64_t v0left = (1 << 10);
229   //ULong64_t v0right = (1 << 11);
230   
231   //triggerMask=esdE->GetTriggerMask();
232   // MB1: SPDFO || V0L || V0R
233   //Bool_t eventTriggered = (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right))); 
234   //MB2: GFO && V0R
235   //triggerMask & spdFO && ((triggerMask&v0left) || (triggerMask&v0right))
236   //Bool_t eventTriggered = (triggerMask & spdFO); 
237  
238   //static AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis();
239   Bool_t eventTriggered = 0;//triggerAnalysis->IsTriggerFired(esdE, AliTriggerAnalysis::kSPDGFO /*| AliTriggerAnalysis::kOfflineFlag*/); 
240
241   // use response of AliPhysicsSelection
242   eventTriggered = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
243
244   Int_t ntracks = esdE->GetNumberOfTracks();
245   Int_t nITS5or6=0,nTPCin=0,nTPCinEta09=0;
246   //printf("Tracks # = %d\n",esdE->GetNumberOfTracks());
247   for(Int_t itr=0; itr<ntracks; itr++) {
248     AliESDtrack *t = esdE->GetTrack(itr);
249     if(t->GetNcls(0)>=5) nITS5or6++;
250     Double_t z0; t->GetZAt(0,esdE->GetMagneticField(),z0);
251     if(t->GetNcls(1)>0 && TMath::Abs(t->GetD(0,0,esdE->GetMagneticField()))<2.8 && TMath::Abs(z0)<20) {
252       nTPCin++;
253       if(TMath::Abs(t->GetTgl())<1.5) nTPCinEta09++;
254     }
255   }
256
257     
258   const AliESDVertex *spdv=esdE->GetPrimaryVertexSPD();
259   const AliESDVertex *tpcv=esdE->GetPrimaryVertexTPC();
260   const AliESDVertex *trkv=esdE->GetPrimaryVertexTracks();
261
262
263   const AliMultiplicity *alimult = esdE->GetMultiplicity();
264   Int_t ntrklets=0,spd0cls=0;
265   if(alimult) {
266     ntrklets = alimult->GetNumberOfTracklets();
267     for(Int_t l=0;l<alimult->GetNumberOfTracklets();l++){
268       if(alimult->GetDeltaPhi(l)<-9998.) ntrklets--;
269     }
270     spd0cls = alimult->GetNumberOfSingleClusters()+ntrklets;
271   }
272   
273   // fill histos
274   
275   if(spdv) {
276     if(spdv->GetNContributors()>0) {
277       TString title=spdv->GetTitle();
278       if(title.Contains("3D")) {
279         fhSPDVertexX->Fill(spdv->GetXv());
280         fhSPDVertexY->Fill(spdv->GetYv());
281       }
282       fhSPDVertexZ->Fill(spdv->GetZv());
283     }
284   }
285   
286   if(trkv) {
287     if(trkv->GetNContributors()>0) {
288       fhTRKVertexX->Fill(trkv->GetXv());
289       fhTRKVertexY->Fill(trkv->GetYv());
290       fhTRKVertexZ->Fill(trkv->GetZv());
291     }
292   }
293   
294   if(tpcv) {
295     if(tpcv->GetNContributors()>0) {
296       fhTPCVertexX->Fill(tpcv->GetXv());
297       fhTPCVertexY->Fill(tpcv->GetYv());
298       fhTPCVertexZ->Fill(tpcv->GetZv());
299     }
300   } 
301   
302
303   // fill ntuple
304   Int_t isize=67;
305   Float_t xnt[67]; for(Int_t iii=0;iii<isize;iii++) xnt[iii]=0.;
306   
307   Int_t index=0;
308
309   xnt[index++]=(Float_t)esdE->GetRunNumber();
310   xnt[index++]=(Float_t)esdE->GetTimeStamp();
311   xnt[index++]=(Float_t)esdE->GetBunchCrossNumber();
312   xnt[index++]=(eventTriggered ? 1. : 0.);
313
314   xnt[index++]=(Float_t)dNchdy;
315
316
317   xnt[index++]=mcVertex[0];
318   xnt[index++]=mcVertex[1];
319   xnt[index++]=mcVertex[2];
320   
321   xnt[index++]=spdv->GetXv();
322   xnt[index++]=spdv->GetXRes();
323   xnt[index++]=spdv->GetYv();
324   xnt[index++]=spdv->GetYRes();
325   xnt[index++]=spdv->GetZv();
326   xnt[index++]=spdv->GetZRes();
327   xnt[index++]=spdv->GetNContributors();
328   TString spdtitle = spdv->GetTitle();
329   xnt[index++]=(spdtitle.Contains("vertexer: 3D") ? 1. : 0.);
330   xnt[index++]=spdv->GetDispersion();
331   
332   xnt[index++]=tpcv->GetXv();
333   xnt[index++]=tpcv->GetXRes();
334   xnt[index++]=tpcv->GetYv();
335   xnt[index++]=tpcv->GetYRes();
336   xnt[index++]=tpcv->GetZv();
337   xnt[index++]=tpcv->GetZRes();
338   xnt[index++]=tpcv->GetNContributors();
339   TString tpctitle = tpcv->GetTitle();
340   xnt[index++]=(tpctitle.Contains("WithConstraint") ? 1. : 0.);
341   
342   xnt[index++]=trkv->GetXv();
343   xnt[index++]=trkv->GetXRes();
344   xnt[index++]=trkv->GetYv();
345   xnt[index++]=trkv->GetYRes();
346   xnt[index++]=trkv->GetZv();
347   xnt[index++]=trkv->GetZRes();
348   xnt[index++]=trkv->GetNContributors();// tpccontrorig;
349   TString trktitle = trkv->GetTitle();
350   xnt[index++]=(trktitle.Contains("WithConstraint") ? 1. : 0.);  
351
352   xnt[index++]=float(ntrklets);
353   xnt[index++]=float(ntracks);
354   xnt[index++]=float(nITS5or6);
355   xnt[index++]=float(nTPCin);
356   xnt[index++]=float(nTPCinEta09);
357
358   xnt[index++]=spd0cls;
359
360   // add recalculated vertices TRK and TPC
361
362   if(fRecoVtxTPC) {
363     AliESDVertex *tpcvnc = ReconstructPrimaryVertexTPC(kFALSE);
364     xnt[index++]=tpcvnc->GetXv();
365     xnt[index++]=tpcvnc->GetXRes();
366     xnt[index++]=tpcvnc->GetYv();
367     xnt[index++]=tpcvnc->GetYRes();
368     xnt[index++]=tpcvnc->GetZv();
369     xnt[index++]=tpcvnc->GetZRes();
370     xnt[index++]=tpcvnc->GetNContributors();
371     delete tpcvnc; tpcvnc=0;
372     
373     AliESDVertex *tpcvc = ReconstructPrimaryVertexTPC(kTRUE);
374     xnt[index++]=tpcvc->GetXv();
375     xnt[index++]=tpcvc->GetXRes();
376     xnt[index++]=tpcvc->GetYv();
377     xnt[index++]=tpcvc->GetYRes();
378     xnt[index++]=tpcvc->GetZv();
379     xnt[index++]=tpcvc->GetZRes();
380     xnt[index++]=tpcvc->GetNContributors();
381     delete tpcvc; tpcvc=0;
382   }
383
384   if(fRecoVtxITSTPC) {
385     AliESDVertex *trkvnc = ReconstructPrimaryVertexITSTPC(kFALSE);
386     xnt[index++]=trkvnc->GetXv();
387     xnt[index++]=trkvnc->GetXRes();
388     xnt[index++]=trkvnc->GetYv();
389     xnt[index++]=trkvnc->GetYRes();
390     xnt[index++]=trkvnc->GetZv();
391     xnt[index++]=trkvnc->GetZRes();
392     xnt[index++]=trkvnc->GetNContributors();
393     delete trkvnc; trkvnc=0;
394     
395     AliESDVertex *trkvc = ReconstructPrimaryVertexITSTPC(kTRUE);
396     xnt[index++]=trkvc->GetXv();
397     xnt[index++]=trkvc->GetXRes();
398     xnt[index++]=trkvc->GetYv();
399     xnt[index++]=trkvc->GetYRes();
400     xnt[index++]=trkvc->GetZv();
401     xnt[index++]=trkvc->GetZRes();
402     xnt[index++]=trkvc->GetNContributors();
403     delete trkvc; trkvc=0;
404   }
405
406
407   if(index>isize) printf("AliAnalysisTaskVertexESD: ERROR, index!=isize\n");
408
409   if(fFillNtuple) fNtupleVertexESD->Fill(xnt);
410   
411   // Post the data already here
412   PostData(1, fOutput);
413
414   return;
415 }      
416
417 //________________________________________________________________________
418 void AliAnalysisTaskVertexESD::Terminate(Option_t *) 
419 {
420   // Draw result to the screen
421   // Called once at the end of the query
422   fOutput = dynamic_cast<TList*> (GetOutputData(1));
423   if (!fOutput) {     
424     Printf("ERROR: fOutput not available");
425     return;
426   }
427
428   fNtupleVertexESD = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleVertexESD"));
429
430   return;
431 }
432
433 //_________________________________________________________________________
434 AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexTPC(Bool_t constr) const {
435   // On the fly reco of TPC vertex from ESD
436   AliESDEvent* evt = (AliESDEvent*) fInputEvent;
437   AliVertexerTracks vertexer(evt->GetMagneticField());
438   vertexer.SetTPCMode(); // defaults
439   Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy);
440   Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0}; 
441   Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
442   AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1);
443   vertexer.SetVtxStart(initVertex);
444   delete initVertex;
445   if(!constr) vertexer.SetConstraintOff();
446
447   return vertexer.FindPrimaryVertex(evt);
448 }
449
450 //_________________________________________________________________________
451 AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexITSTPC(Bool_t constr) const {
452   // On the fly reco of ITS+TPC vertex from ESD
453   AliESDEvent* evt = (AliESDEvent*) fInputEvent;
454   AliVertexerTracks vertexer(evt->GetMagneticField());
455   vertexer.SetITSMode(); // defaults
456   vertexer.SetMinClusters(4); // default is 5
457   Float_t diamondcovxy[3]; evt->GetDiamondCovXY(diamondcovxy);
458   Double_t pos[3]={evt->GetDiamondX(),evt->GetDiamondY(),0}; 
459   Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
460   AliESDVertex *initVertex = new AliESDVertex(pos,cov,1.,1);
461   vertexer.SetVtxStart(initVertex);
462   delete initVertex;
463   if(!constr) vertexer.SetConstraintOff();
464
465   // use only ITS-TPC or only ITS-SA tracks
466   if(fOnlyITSTPCTracks || fOnlyITSSATracks) {
467     Int_t iskip=0;
468     Int_t *skip = new Int_t[evt->GetNumberOfTracks()];
469     for(Int_t itr=0;itr<evt->GetNumberOfTracks(); itr++) {
470       AliESDtrack* track = evt->GetTrack(itr);
471       if(fOnlyITSTPCTracks && track->GetNcls(1)==0) { // skip ITSSA
472         skip[iskip++]=itr;
473       }
474       if(fOnlyITSSATracks && track->GetNcls(1)>0) { // skip ITSTPC
475         skip[iskip++]=itr;
476       }
477     }
478     vertexer.SetSkipTracks(iskip,skip);
479     delete [] skip; skip=NULL;
480   }
481
482   return vertexer.FindPrimaryVertex(evt);
483 }