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