Remove access to ESDfriends
[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
50 #include "AliMCEventHandler.h"
51 #include "AliMCEvent.h"
52 #include "AliStack.h"
53 #include "AliLog.h"
54
55 #include "AliGenEventHeader.h" 
56 #include "AliAnalysisTaskVertexESD.h"
57
58
59 ClassImp(AliAnalysisTaskVertexESD)
60
61 //________________________________________________________________________
62 AliAnalysisTaskVertexESD::AliAnalysisTaskVertexESD(const char *name) : 
63 AliAnalysisTask(name, "VertexESDTask"), 
64 fReadMC(kFALSE),
65 fESD(0), 
66 fOutput(0), 
67 fNtupleVertexESD(0),
68 fhTrackRefs(0)
69 {
70   // Constructor
71
72   // Define input and output slots here
73   // Input slot #0 works with a TChain
74   DefineInput(0, TChain::Class());
75   // Output slot #0 writes into a TList container
76   DefineOutput(0, TList::Class());  //My private output
77 }
78 //________________________________________________________________________
79 AliAnalysisTaskVertexESD::~AliAnalysisTaskVertexESD()
80 {
81   // Destructor
82
83   // histograms are in the output list and deleted when the output
84   // list is deleted by the TSelector dtor
85
86   if (fOutput) {
87     delete fOutput;
88     fOutput = 0;
89   }
90 }
91
92 //________________________________________________________________________
93 void AliAnalysisTaskVertexESD::ConnectInputData(Option_t *) 
94 {
95   // Connect ESD or AOD here
96   // Called once
97
98   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
99   if (!tree) {
100     Printf("ERROR: Could not read chain from input slot 0");
101   } else {
102
103     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
104     
105     if (!esdH) {
106       Printf("ERROR: Could not get ESDInputHandler");
107     } else {
108       fESD = esdH->GetEvent();
109     }
110
111   }
112   
113   return;
114 }
115
116 //________________________________________________________________________
117 void AliAnalysisTaskVertexESD::CreateOutputObjects()
118 {
119   // Create histograms
120   // Called once
121
122   // Several histograms are more conveniently managed in a TList
123   fOutput = new TList;
124   fOutput->SetOwner();
125
126   fNtupleVertexESD = new TNtuple("fNtupleVertexESD","vertices","xtrue:ytrue:ztrue:xSPD:xdiffSPD:xerrSPD:ySPD:ydiffSPD:yerrSPD:zSPD:zdiffSPD:zerrSPD:ntrksSPD:xTPC:xdiffTPC:xerrTPC:yTPC:ydiffTPC:yerrTPC:zTPC:zdiffTPC:zerrTPC:ntrksTPC:xTRK:xdiffTRK:xerrTRK:yTRK:ydiffTRK:yerrTRK:zTRK:zdiffTRK:zerrTRK:ntrksTRK:ntrklets:nESDtracks:nITSrefit5or6:nTPCin:nTPCinEta09:dndygen:triggered:SPD3D:SPD0cls:constrTRK");
127
128   fOutput->Add(fNtupleVertexESD);
129
130   fhTrackRefs = new TH2F("fhTrackRefs","Track references; x; y",1000,-4,4,1000,-4,4);
131   fOutput->Add(fhTrackRefs);
132
133  return;
134 }
135
136 //________________________________________________________________________
137 void AliAnalysisTaskVertexESD::Exec(Option_t *) 
138 {
139   // Main loop
140   // Called for each event
141   
142   if (!fESD) {
143     Printf("ERROR: fESD not available");
144     return;
145   }
146
147   
148   TArrayF mcVertex(3);
149   mcVertex[0]=9999.; mcVertex[1]=9999.; mcVertex[2]=9999.;
150   Float_t dNchdy=-999.;
151
152   // ***********  MC info ***************
153   if (fReadMC) {
154     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
155     if (!eventHandler) {
156       Printf("ERROR: Could not retrieve MC event handler");
157       return;
158     }
159     
160     AliMCEvent* mcEvent = eventHandler->MCEvent();
161     if (!mcEvent) {
162       Printf("ERROR: Could not retrieve MC event");
163       return;
164     }
165     
166     AliStack* stack = mcEvent->Stack();
167     if (!stack) {
168       AliDebug(AliLog::kError, "Stack not available");
169       return;
170     }
171     
172     AliHeader* header = mcEvent->Header();
173     if (!header) {
174       AliDebug(AliLog::kError, "Header not available");
175       return;
176     }
177     AliGenEventHeader* genHeader = header->GenEventHeader();
178     genHeader->PrimaryVertex(mcVertex);
179
180     Int_t ngenpart = (Int_t)stack->GetNtrack();
181     //printf("# generated particles = %d\n",ngenpart);
182     dNchdy=0;
183     for(Int_t ip=0; ip<ngenpart; ip++) {
184       TParticle* part = (TParticle*)stack->Particle(ip);
185       // keep only electorns, muons, pions, kaons and protons
186       Int_t apdg = TMath::Abs(part->GetPdgCode());
187       if(apdg!=11 && apdg!=13 && apdg!=211 && apdg!=321 && apdg!=2212) continue;      
188       // reject secondaries
189       if(TMath::Sqrt((part->Vx()-mcVertex[0])*(part->Vx()-mcVertex[0])+(part->Vy()-mcVertex[1])*(part->Vy()-mcVertex[1]))>0.0010) continue;
190       // reject incoming protons
191       Double_t energy  = part->Energy();
192       if(energy>900.) continue;
193       Double_t pz = part->Pz();
194       Double_t y = 0.5*TMath::Log((energy+pz+1.e-13)/(energy-pz+1.e-13));
195       if(TMath::Abs(y)<1.0) dNchdy += 0.5; // count 1/2 of particles in |y|<1
196       // tracks refs
197       TClonesArray *trefs=0;
198       Int_t ntrefs = mcEvent->GetParticleAndTR(ip,part,trefs);
199       if(ntrefs<0) continue;
200       for(Int_t iref=0; iref<ntrefs; iref++) {
201         AliTrackReference *tref = (AliTrackReference*)trefs->At(iref);
202         if(tref->R()>10.) continue;
203         fhTrackRefs->Fill(tref->X(),tref->Y());
204       }
205     }
206     //printf("# primary particles = %7.1f\n",dNchdy);
207   } 
208   // ***********  MC info ***************
209
210     
211   // Trigger
212   ULong64_t triggerMask;
213   ULong64_t spdFO = (1 << 14);
214   ULong64_t v0left = (1 << 11);
215   ULong64_t v0right = (1 << 12);
216   
217   triggerMask=fESD->GetTriggerMask();
218   // MB1: SPDFO || V0L || V0R
219   Bool_t eventTriggered = (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right))); 
220   //MB2: GFO && V0R
221   //triggerMask & spdFO && ((triggerMask&v0left) || (triggerMask&v0right))
222  
223
224   Int_t ntracks = fESD->GetNumberOfTracks();
225   Int_t nITS5or6=0,nTPCin=0,nTPCinEta09=0;
226   //printf("Tracks # = %d\n",fESD->GetNumberOfTracks());
227   for(Int_t itr=0; itr<ntracks; itr++) {
228     AliESDtrack *t = fESD->GetTrack(itr);
229     if(t->GetNcls(0)>=5) nITS5or6++;
230     Double_t z0; t->GetZAt(0,fESD->GetMagneticField(),z0);
231     if(t->GetNcls(1)>0 && TMath::Abs(t->GetD(0,0,fESD->GetMagneticField()))<2.8 && TMath::Abs(z0)<20) {
232       nTPCin++;
233       if(TMath::Abs(t->GetTgl())<1.5) nTPCinEta09++;
234     }
235   }
236
237     
238   const AliESDVertex *spdv=fESD->GetPrimaryVertexSPD();
239   const AliESDVertex *tpcv=fESD->GetPrimaryVertexTPC();
240   const AliESDVertex *trkv=fESD->GetPrimaryVertex();
241
242   //Float_t tpccontrorig=tpcv->GetNContributors();
243   //tpcv = 0;
244   //tpcv = ReconstructPrimaryVertexTPC();
245
246   const AliMultiplicity *alimult = fESD->GetMultiplicity();
247   Int_t ntrklets=0,spd0cls=0;
248   if(alimult) {
249     ntrklets = alimult->GetNumberOfTracklets();
250     for(Int_t l=0;l<alimult->GetNumberOfTracklets();l++){
251       if(alimult->GetDeltaPhi(l)<-9998.) ntrklets--;
252     }
253     spd0cls = alimult->GetNumberOfSingleClusters()+ntrklets;
254   }
255
256
257   Float_t xnt[43];
258   
259   xnt[0]=mcVertex[0];
260   xnt[1]=mcVertex[1];
261   xnt[2]=mcVertex[2];
262   
263   xnt[3]=spdv->GetXv();
264   xnt[4]=(spdv->GetXv()-mcVertex[0])*10000.;
265   xnt[5]=spdv->GetXRes();
266   xnt[6]=spdv->GetYv();
267   xnt[7]=(spdv->GetYv()-mcVertex[1])*10000.;
268   xnt[8]=spdv->GetYRes();
269   xnt[9]=spdv->GetZv();
270   xnt[10]=(spdv->GetZv()-mcVertex[2])*10000.;
271   xnt[11]=spdv->GetZRes();
272   xnt[12]=spdv->GetNContributors();
273   
274   xnt[13]=tpcv->GetXv();
275   xnt[14]=(tpcv->GetXv()-mcVertex[0])*10000.;
276   xnt[15]=tpcv->GetXRes();
277   xnt[16]=tpcv->GetYv();
278   xnt[17]=(tpcv->GetYv()-mcVertex[1])*10000.;
279   xnt[18]=tpcv->GetYRes();
280   xnt[19]=tpcv->GetZv();
281   xnt[20]=(tpcv->GetZv()-mcVertex[2])*10000.;
282   xnt[21]=tpcv->GetZRes();
283   xnt[22]=tpcv->GetNContributors();
284   
285   xnt[23]=trkv->GetXv();
286   xnt[24]=(trkv->GetXv()-mcVertex[0])*10000.;
287   xnt[25]=trkv->GetXRes();
288   xnt[26]=trkv->GetYv();
289   xnt[27]=(trkv->GetYv()-mcVertex[1])*10000.;
290   xnt[28]=trkv->GetYRes();
291   xnt[29]=trkv->GetZv();
292   xnt[30]=(trkv->GetZv()-mcVertex[2])*10000.;
293   xnt[31]=trkv->GetZRes();
294   xnt[32]=trkv->GetNContributors();// tpccontrorig;
295   
296
297   xnt[33]=float(ntrklets);
298   xnt[34]=float(ntracks);
299   xnt[35]=float(nITS5or6);
300   xnt[36]=float(nTPCin);
301   xnt[37]=float(nTPCinEta09);
302
303   xnt[38]=float(dNchdy);
304
305   xnt[39]=0.; 
306   if(eventTriggered) xnt[39]=1.;
307
308   xnt[40]=0.; 
309   TString spdtitle = spdv->GetTitle();
310   if(spdtitle.Contains("vertexer: 3D")) xnt[40]=1.;
311
312   xnt[42]=0.; 
313   TString trktitle = trkv->GetTitle();
314   if(trktitle.Contains("WithConstraint")) xnt[42]=1.;
315
316   xnt[41]=spd0cls;
317
318   fNtupleVertexESD->Fill(xnt);
319     
320   // Post the data already here
321   PostData(0, fOutput);
322
323   return;
324 }      
325
326 //________________________________________________________________________
327 void AliAnalysisTaskVertexESD::Terminate(Option_t *) 
328 {
329   // Draw result to the screen
330   // Called once at the end of the query
331   fOutput = dynamic_cast<TList*> (GetOutputData(0));
332   if (!fOutput) {     
333     Printf("ERROR: fOutput not available");
334     return;
335   }
336
337   fNtupleVertexESD = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleVertexESD"));
338
339   return;
340 }
341
342 //_________________________________________________________________________
343 AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexTPC() {
344   // On the fly reco of TPC vertex from ESD
345
346   AliVertexerTracks vertexer(fESD->GetMagneticField());
347   vertexer.SetTPCMode(0.1,1.0,5.,0,1,3.,0.1,1.5); // defaults
348   //vertexer.SetTPCMode(0.1,1.0,5.,0,1,3.,0.1,1.5); // defaults
349   Double_t pos[3]={+0.0220,-0.0340,+0.270}; 
350   Double_t err[3]={0.0200,0.0200,7.5};
351   AliESDVertex *initVertex = new AliESDVertex(pos,err);
352   vertexer.SetVtxStart(initVertex);
353   delete initVertex;
354   //vertexer.SetConstraintOff();
355
356   return vertexer.FindPrimaryVertex(fESD);
357 }