Added correlation histos
[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 fRecoVtxTPC(kFALSE),
66 fRecoVtxITSTPC(kFALSE),
67 fESD(0), 
68 fOutput(0), 
69 fNtupleVertexESD(0),
70 fhTrackRefs(0),
71 fhVtxTPCx(0),
72 fhVtxTPCy(0),
73 fhVtxTPCz(0),
74 fhVtxTRKx(0),
75 fhVtxTRKy(0),
76 fhVtxTRKz(0),
77 fhVtxSPDx(0),
78 fhVtxSPDy(0),
79 fhVtxSPDz3D(0),
80 fhVtxSPDzZ(0),
81 fhVtxTPCvsSPDx(0),
82 fhVtxTPCvsSPDy(0),
83 fhVtxTPCvsSPDz(0),
84 fhVtxTRKvsSPDx(0),
85 fhVtxTRKvsSPDy(0),
86 fhVtxTRKvsSPDz(0),
87 fhVtxSPDContrvsMult(0),
88 fhVtxTRKContrvsTrks56(0),
89 fhVtxTPCContrvsTrks(0)
90 {
91   // Constructor
92
93   // Define input and output slots here
94   // Input slot #0 works with a TChain
95   DefineInput(0, TChain::Class());
96   // Output slot #0 writes into a TList container
97   DefineOutput(0, TList::Class());  //My private output
98 }
99 //________________________________________________________________________
100 AliAnalysisTaskVertexESD::~AliAnalysisTaskVertexESD()
101 {
102   // Destructor
103
104   // histograms are in the output list and deleted when the output
105   // list is deleted by the TSelector dtor
106
107   if (fOutput) {
108     delete fOutput;
109     fOutput = 0;
110   }
111 }
112
113 //________________________________________________________________________
114 void AliAnalysisTaskVertexESD::ConnectInputData(Option_t *) 
115 {
116   // Connect ESD or AOD here
117   // Called once
118
119   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
120   if (!tree) {
121     Printf("ERROR: Could not read chain from input slot 0");
122   } else {
123
124     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
125     
126     if (!esdH) {
127       Printf("ERROR: Could not get ESDInputHandler");
128     } else {
129       fESD = esdH->GetEvent();
130     }
131
132   }
133   
134   return;
135 }
136
137 //________________________________________________________________________
138 void AliAnalysisTaskVertexESD::CreateOutputObjects()
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:xtrue:ytrue:ztrue:xSPD:xerrSPD:ySPD:yerrSPD:zSPD:zerrSPD:ntrksSPD:xTPC:xerrTPC:yTPC:yerrTPC:zTPC:zerrTPC:ntrksTPC:xTRK:xerrTRK:yTRK:yerrTRK:zTRK:zerrTRK:ntrksTRK:ntrklets:nESDtracks:nITSrefit5or6:nTPCin:nTPCinEta09:dndygen:triggered:SPD3D:SPD0cls:constrTRK:constrTPC");
148
149   fOutput->Add(fNtupleVertexESD);
150
151   fhTrackRefs = new TH2F("fhTrackRefs","Track references; x; y",1000,-4,4,1000,-4,4);
152   fOutput->Add(fhTrackRefs);
153
154   fhVtxTPCx = new TH1F("fhVtxTPCx","VtxTPC; x [cm]; events",1000,-1,1);
155   fOutput->Add(fhVtxTPCx);
156   fhVtxTPCy = new TH1F("fhVtxTPCy","VtxTPC; y [cm]; events",1000,-1,1);
157   fOutput->Add(fhVtxTPCx);
158   fhVtxTPCz = new TH1F("fhVtxTPCz","VtxTPC; z [cm]; events",1000,-20,20);
159   fOutput->Add(fhVtxTPCz);
160   fhVtxTRKx = new TH1F("fhVtxTRKx","VtxTRK; x [cm]; events",1000,-1,1);
161   fOutput->Add(fhVtxTRKx);
162   fhVtxTRKy = new TH1F("fhVtxTRKy","VtxTRK; y [cm]; events",1000,-1,1);
163   fOutput->Add(fhVtxTRKx);
164   fhVtxTRKz = new TH1F("fhVtxTRKz","VtxTRK; z [cm]; events",1000,-20,20);
165   fOutput->Add(fhVtxTRKz);
166   fhVtxSPDx = new TH1F("fhVtxSPDx","VtxSPD3D; x [cm]; events",1000,-1,1);
167   fOutput->Add(fhVtxSPDx);
168   fhVtxSPDy = new TH1F("fhVtxSPDy","VtxSPD3D; y [cm]; events",1000,-1,1);
169   fOutput->Add(fhVtxSPDx);
170   fhVtxSPDz3D = new TH1F("fhVtxSPDz3D","VtxSPD3D; z [cm]; events",1000,-20,20);
171   fOutput->Add(fhVtxSPDz3D);
172   fhVtxSPDzZ = new TH1F("fhVtxSPDzZ","VtxSPDZ; z [cm]; events",1000,-20,20);
173   fOutput->Add(fhVtxSPDzZ);
174
175   fhVtxTPCvsSPDx = new TH2F("fhVtxTPCvsSPDx","TPC vs SPD ; x SPD [cm]; x TPC [cm]",100,-.1,.1,100,-1,1);
176   fOutput->Add(fhVtxTPCvsSPDx);
177   fhVtxTPCvsSPDy = new TH2F("fhVtxTPCvsSPDy","TPC vs SPD ; y SPD [cm]; y TPC [cm]",100,-.1,.1,100,-1,1);
178   fOutput->Add(fhVtxTPCvsSPDy);
179   fhVtxTPCvsSPDz = new TH2F("fhVtxTPCvsSPDz","TPC vs SPD ; z SPD [cm]; z TPC [cm]",100,-20,20,100,20,20);
180   fOutput->Add(fhVtxTPCvsSPDz);
181
182   fhVtxTRKvsSPDx = new TH2F("fhVtxTRKvsSPDx","TRK vs SPD ; x SPD [cm]; x TRK [cm]",100,-.1,.1,100,-.1,.1);
183   fOutput->Add(fhVtxTRKvsSPDx);
184   fhVtxTRKvsSPDy = new TH2F("fhVtxTRKvsSPDy","TRK vs SPD ; y SPD [cm]; y TRK [cm]",100,-.1,.1,100,-.1,.1);
185   fOutput->Add(fhVtxTRKvsSPDy);
186   fhVtxTRKvsSPDz = new TH2F("fhVtxTRKvsSPDz","TRK vs SPD ; z SPD [cm]; z TRK [cm]",100,-20,20,100,20,20);
187   fOutput->Add(fhVtxTRKvsSPDz);
188
189   fhVtxSPDContrvsMult = new TH2F("fhVtxSPDContrvsMult","SPD vertex: contributors VS SPD tracklets; contributors; SPD tracklets (AliMult)",100,-0.5,99.5,100,-0.5,99.5);
190   fOutput->Add(fhVtxSPDContrvsMult);
191
192   fhVtxTRKContrvsTrks56 = new TH2F("fhVtxTRKContrvsTrks56","TRK vertex: contributors VS trks with #ge 5 ITS cls; contributors; tracks",100,-0.5,99.5,100,-0.5,99.5);
193   fOutput->Add(fhVtxTRKContrvsTrks56);
194
195   fhVtxTPCContrvsTrks = new TH2F("fhVtxTPCContrvsTrks","TPC vertex: contributors VS TPC trks (all); contributors; tracks",100,-0.5,99.5,100,-0.5,99.5);
196   fOutput->Add(fhVtxTPCContrvsTrks);
197
198   return;
199 }
200
201 //________________________________________________________________________
202 void AliAnalysisTaskVertexESD::Exec(Option_t *) 
203 {
204   // Main loop
205   // Called for each event
206   
207   if (!fESD) {
208     Printf("ERROR: fESD not available");
209     return;
210   }
211
212   
213   TArrayF mcVertex(3);
214   mcVertex[0]=9999.; mcVertex[1]=9999.; mcVertex[2]=9999.;
215   Float_t dNchdy=-999.;
216
217   // ***********  MC info ***************
218   if (fReadMC) {
219     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
220     if (!eventHandler) {
221       Printf("ERROR: Could not retrieve MC event handler");
222       return;
223     }
224     
225     AliMCEvent* mcEvent = eventHandler->MCEvent();
226     if (!mcEvent) {
227       Printf("ERROR: Could not retrieve MC event");
228       return;
229     }
230     
231     AliStack* stack = mcEvent->Stack();
232     if (!stack) {
233       AliDebug(AliLog::kError, "Stack not available");
234       return;
235     }
236     
237     AliHeader* header = mcEvent->Header();
238     if (!header) {
239       AliDebug(AliLog::kError, "Header not available");
240       return;
241     }
242     AliGenEventHeader* genHeader = header->GenEventHeader();
243     genHeader->PrimaryVertex(mcVertex);
244
245     Int_t ngenpart = (Int_t)stack->GetNtrack();
246     //printf("# generated particles = %d\n",ngenpart);
247     dNchdy=0;
248     for(Int_t ip=0; ip<ngenpart; ip++) {
249       TParticle* part = (TParticle*)stack->Particle(ip);
250       // keep only electorns, muons, pions, kaons and protons
251       Int_t apdg = TMath::Abs(part->GetPdgCode());
252       if(apdg!=11 && apdg!=13 && apdg!=211 && apdg!=321 && apdg!=2212) continue;      
253       // reject secondaries
254       if(TMath::Sqrt((part->Vx()-mcVertex[0])*(part->Vx()-mcVertex[0])+(part->Vy()-mcVertex[1])*(part->Vy()-mcVertex[1]))>0.0010) continue;
255       // reject incoming protons
256       Double_t energy  = part->Energy();
257       if(energy>900.) continue;
258       Double_t pz = part->Pz();
259       Double_t y = 0.5*TMath::Log((energy+pz+1.e-13)/(energy-pz+1.e-13));
260       if(TMath::Abs(y)<1.0) dNchdy += 0.5; // count 1/2 of particles in |y|<1
261       // tracks refs
262       TClonesArray *trefs=0;
263       Int_t ntrefs = mcEvent->GetParticleAndTR(ip,part,trefs);
264       if(ntrefs<0) continue;
265       for(Int_t iref=0; iref<ntrefs; iref++) {
266         AliTrackReference *tref = (AliTrackReference*)trefs->At(iref);
267         if(tref->R()>10.) continue;
268         fhTrackRefs->Fill(tref->X(),tref->Y());
269       }
270     }
271     //printf("# primary particles = %7.1f\n",dNchdy);
272   } 
273   // ***********  MC info ***************
274
275     
276   // Trigger
277   ULong64_t triggerMask;
278   ULong64_t spdFO = (1 << 14);
279   ULong64_t v0left = (1 << 11);
280   ULong64_t v0right = (1 << 12);
281   
282   triggerMask=fESD->GetTriggerMask();
283   // MB1: SPDFO || V0L || V0R
284   Bool_t eventTriggered = (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right))); 
285   //MB2: GFO && V0R
286   //triggerMask & spdFO && ((triggerMask&v0left) || (triggerMask&v0right))
287  
288
289   Int_t ntracks = fESD->GetNumberOfTracks();
290   Int_t nITS5or6=0,nTPCin=0,nTPCinEta09=0;
291   //printf("Tracks # = %d\n",fESD->GetNumberOfTracks());
292   for(Int_t itr=0; itr<ntracks; itr++) {
293     AliESDtrack *t = fESD->GetTrack(itr);
294     if(t->GetNcls(0)>=5) nITS5or6++;
295     Double_t z0; t->GetZAt(0,fESD->GetMagneticField(),z0);
296     if(t->GetNcls(1)>0 && TMath::Abs(t->GetD(0,0,fESD->GetMagneticField()))<2.8 && TMath::Abs(z0)<20) {
297       nTPCin++;
298       if(TMath::Abs(t->GetTgl())<1.5) nTPCinEta09++;
299     }
300   }
301
302     
303   const AliESDVertex *spdv=fESD->GetPrimaryVertexSPD();
304   const AliESDVertex *tpcv=fESD->GetPrimaryVertexTPC();
305   const AliESDVertex *trkv=fESD->GetPrimaryVertex();
306
307   //Float_t tpccontrorig=tpcv->GetNContributors();
308
309   if(fRecoVtxTPC) {
310     tpcv = 0;
311     tpcv = ReconstructPrimaryVertexTPC();
312   }
313   if(fRecoVtxITSTPC) {
314     trkv = 0;
315     trkv = ReconstructPrimaryVertexITSTPC();
316   }
317
318   const AliMultiplicity *alimult = fESD->GetMultiplicity();
319   Int_t ntrklets=0,spd0cls=0;
320   if(alimult) {
321     ntrklets = alimult->GetNumberOfTracklets();
322     for(Int_t l=0;l<alimult->GetNumberOfTracklets();l++){
323       if(alimult->GetDeltaPhi(l)<-9998.) ntrklets--;
324     }
325     spd0cls = alimult->GetNumberOfSingleClusters()+ntrklets;
326   }
327
328
329   Int_t isize=37;
330   Float_t xnt[37];
331   
332   Int_t index=0;
333
334   xnt[index++]=(Float_t)fESD->GetRunNumber();
335   xnt[index++]=(Float_t)fESD->GetTimeStamp();
336
337   xnt[index++]=mcVertex[0];
338   xnt[index++]=mcVertex[1];
339   xnt[index++]=mcVertex[2];
340   
341   xnt[index++]=spdv->GetXv();
342   xnt[index++]=spdv->GetXRes();
343   xnt[index++]=spdv->GetYv();
344   xnt[index++]=spdv->GetYRes();
345   xnt[index++]=spdv->GetZv();
346   xnt[index++]=spdv->GetZRes();
347   xnt[index++]=spdv->GetNContributors();
348   
349   xnt[index++]=tpcv->GetXv();
350   xnt[index++]=tpcv->GetXRes();
351   xnt[index++]=tpcv->GetYv();
352   xnt[index++]=tpcv->GetYRes();
353   xnt[index++]=tpcv->GetZv();
354   xnt[index++]=tpcv->GetZRes();
355   xnt[index++]=tpcv->GetNContributors();
356   
357   xnt[index++]=trkv->GetXv();
358   xnt[index++]=trkv->GetXRes();
359   xnt[index++]=trkv->GetYv();
360   xnt[index++]=trkv->GetYRes();
361   xnt[index++]=trkv->GetZv();
362   xnt[index++]=trkv->GetZRes();
363   xnt[index++]=trkv->GetNContributors();// tpccontrorig;
364   
365
366   xnt[index++]=float(ntrklets);
367   xnt[index++]=float(ntracks);
368   xnt[index++]=float(nITS5or6);
369   xnt[index++]=float(nTPCin);
370   xnt[index++]=float(nTPCinEta09);
371
372   xnt[index++]=float(dNchdy);
373
374   xnt[index++]=(eventTriggered ? 1. : 0.);
375
376   TString spdtitle = spdv->GetTitle();
377   xnt[index++]=(spdtitle.Contains("vertexer: 3D") ? 1. : 0.);
378
379   xnt[index++]=spd0cls;
380
381   TString trktitle = trkv->GetTitle();
382   xnt[index++]=(trktitle.Contains("WithConstraint") ? 1. : 0.);
383
384   TString tpctitle = tpcv->GetTitle();
385   xnt[index++]=(tpctitle.Contains("WithConstraint") ? 1. : 0.);
386
387
388   if(index!=isize) printf("AliAnalysisTaskVertexESD: ERROR, index!=isize\n");
389
390   fNtupleVertexESD->Fill(xnt);
391   
392   fhVtxTPCx->Fill(tpcv->GetXv());
393   fhVtxTPCy->Fill(tpcv->GetYv());
394   fhVtxTPCz->Fill(tpcv->GetZv());
395   fhVtxTRKx->Fill(trkv->GetXv());
396   fhVtxTRKy->Fill(trkv->GetYv());
397   fhVtxTRKz->Fill(trkv->GetZv());
398   fhVtxSPDx->Fill(spdv->GetXv());
399   fhVtxSPDy->Fill(spdv->GetYv());
400   if(spdtitle.Contains("vertexer: 3D")) {
401     fhVtxSPDz3D->Fill(spdv->GetZv());
402   } else {
403     fhVtxSPDzZ->Fill(spdv->GetZv());
404   }
405   fhVtxTPCvsSPDx->Fill(spdv->GetXv(),tpcv->GetXv());
406   fhVtxTPCvsSPDy->Fill(spdv->GetYv(),tpcv->GetYv());
407   fhVtxTPCvsSPDz->Fill(spdv->GetZv(),tpcv->GetZv());
408   fhVtxTRKvsSPDx->Fill(spdv->GetXv(),trkv->GetXv());
409   fhVtxTRKvsSPDy->Fill(spdv->GetYv(),trkv->GetYv());
410   fhVtxTRKvsSPDz->Fill(spdv->GetZv(),trkv->GetZv());
411   fhVtxSPDContrvsMult->Fill(ntrklets,spdv->GetNContributors());
412   fhVtxTRKContrvsTrks56->Fill(nITS5or6,trkv->GetNContributors());   
413   fhVtxTPCContrvsTrks->Fill(nTPCin,tpcv->GetNContributors());
414   
415   // Post the data already here
416   PostData(0, fOutput);
417
418   return;
419 }      
420
421 //________________________________________________________________________
422 void AliAnalysisTaskVertexESD::Terminate(Option_t *) 
423 {
424   // Draw result to the screen
425   // Called once at the end of the query
426   fOutput = dynamic_cast<TList*> (GetOutputData(0));
427   if (!fOutput) {     
428     Printf("ERROR: fOutput not available");
429     return;
430   }
431
432   fNtupleVertexESD = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleVertexESD"));
433
434   return;
435 }
436
437 //_________________________________________________________________________
438 AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexTPC() const {
439   // On the fly reco of TPC vertex from ESD
440
441   AliVertexerTracks vertexer(fESD->GetMagneticField());
442   vertexer.SetTPCMode(); // defaults
443   //vertexer.SetTPCMode(0.1,1.0,5.,0,1,3.,0.1,1.5);
444   Double_t pos[3]={+0.0220,-0.0340,+0.270}; 
445   Double_t err[3]={0.0200,0.0200,7.5};
446   AliESDVertex *initVertex = new AliESDVertex(pos,err);
447   vertexer.SetVtxStart(initVertex);
448   delete initVertex;
449   vertexer.SetConstraintOff();
450
451   return vertexer.FindPrimaryVertex(fESD);
452 }
453
454 //_________________________________________________________________________
455 AliESDVertex* AliAnalysisTaskVertexESD::ReconstructPrimaryVertexITSTPC() const {
456   // On the fly reco of ITS+TPC vertex from ESD
457
458   AliVertexerTracks vertexer(fESD->GetMagneticField());
459   vertexer.SetITSMode(); // defaults
460   //vertexer.SetTPCMode(0.1,1.0,5.,0,1,3.,0.1,1.5);
461   Double_t pos[3]={+0.0220,-0.0340,+0.270}; 
462   Double_t err[3]={0.0200,0.0200,7.5};
463   AliESDVertex *initVertex = new AliESDVertex(pos,err);
464   vertexer.SetVtxStart(initVertex);
465   delete initVertex;
466   vertexer.SetConstraintOff();
467
468   return vertexer.FindPrimaryVertex(fESD);
469 }