Adding abstract class for comparison components
[u/mrichter/AliRoot.git] / PWG1 / AliComparisonDCA.cxx
1 //------------------------------------------------------------------------------
2 // Implementation of AliComparisonDCA class. It keeps information from 
3 // comparison of reconstructed and MC particle tracks. In addtion, 
4 // it keeps selection cuts used during comparison. The comparison 
5 // information is stored in the ROOT histograms. Analysis of these 
6 // histograms can be done by using Analyse() class function. The result of 
7 // the analysis (histograms/graphs) are stored in the folder
8 // which is a data member of AliComparisonDCA.
9 //  
10 // Author: J.Otwinowski 04/02/2008 
11 //------------------------------------------------------------------------------
12
13 /*
14  
15   // after running comparison task, read the file, and get component
16   gSystem->Load("libPWG1.so");
17   TFile f("Output.root");
18   AliComparisonDCA * compObj = (AliComparisonDCA*)f.Get("AliComparisonDCA");
19
20   // analyse comparison data
21   compObj->Analyse();
22
23   // the output histograms/graphs will be stored in the folder "folderDCA" 
24   compObj->GetAnalysisFolder()->ls("*");
25  
26   // user can save whole comparison object (or only folder with anlysed histograms) 
27   // in the seperate output file (e.g.)
28   TFile fout("Analysed_DCA.root","recreate");
29   compObj->Write(); // compObj->GetAnalysisFolder()->Write();
30   fout.Close();
31
32 */
33
34 #include <iostream>
35
36 #include "TFile.h"
37 #include "TCint.h"
38 #include "TH3F.h"
39 #include "TH2F.h"
40 #include "TF1.h"
41 #include "TProfile.h"
42 #include "TProfile2D.h"
43 #include "TGraph2D.h"
44 #include "TCanvas.h"
45 #include "TGraph.h"
46 // 
47 #include "AliTracker.h"   
48 #include "AliESDEvent.h"   
49 #include "AliESD.h"
50 #include "AliESDfriend.h"
51 #include "AliESDfriendTrack.h"
52 #include "AliRecInfoCuts.h" 
53 #include "AliMCInfoCuts.h" 
54 #include "AliLog.h" 
55 #include "AliESDVertex.h" 
56 //
57 #include "AliMathBase.h"
58 #include "AliTreeDraw.h" 
59
60 #include "AliMCInfo.h" 
61 #include "AliESDRecInfo.h" 
62 #include "AliComparisonDCA.h" 
63
64 using namespace std;
65
66 ClassImp(AliComparisonDCA)
67
68 //_____________________________________________________________________________
69 AliComparisonDCA::AliComparisonDCA():
70 //  TNamed("AliComparisonDCA","AliComparisonDCA"),
71   AliComparisonObject("AliComparisonDCA"),
72
73   // DCA histograms
74   fD0TanSPtB1(0),
75   fD1TanSPtB1(0),
76   fD0TanSPtL1(0),
77   fD1TanSPtL1(0),
78   fD0TanSPtInTPC(0),
79   fD1TanSPtInTPC(0),
80   fVertex(0),
81
82   // Cuts 
83   fCutsRC(0), 
84   fCutsMC(0),  
85
86   // histogram folder 
87   fAnalysisFolder(0)
88 {
89   Init();
90
91   // vertex (0,0,0)
92   fVertex = new AliESDVertex();
93   fVertex->SetXv(0.0);
94   fVertex->SetYv(0.0);
95   fVertex->SetZv(0.0);
96 }
97
98 //_____________________________________________________________________________
99 AliComparisonDCA::~AliComparisonDCA()
100 {
101   //
102   if(fD0TanSPtB1) delete fD0TanSPtB1; fD0TanSPtB1=0;
103   if(fD1TanSPtB1) delete fD1TanSPtB1; fD1TanSPtB1=0;
104   if(fD0TanSPtL1) delete fD0TanSPtL1; fD0TanSPtL1=0;
105   if(fD1TanSPtL1) delete fD1TanSPtL1; fD1TanSPtL1=0;
106   if(fD0TanSPtInTPC) delete fD0TanSPtInTPC; fD0TanSPtInTPC=0;
107   if(fD1TanSPtInTPC) delete fD1TanSPtInTPC; fD1TanSPtInTPC=0;
108   if(fVertex) delete fVertex; fVertex=0;
109   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
110
111 }
112
113 //_____________________________________________________________________________
114 void AliComparisonDCA::Init()
115 {
116   // DCA histograms
117   fD0TanSPtB1 = new TH3F("DCAyTanSPtB1","DCAyTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);
118   fD0TanSPtB1->SetXTitle("tan(#theta)");
119   fD0TanSPtB1->SetYTitle("#sqrt{p_{t}(GeV/c)}");
120   fD0TanSPtB1->SetZTitle("DCA_{xy}");
121
122   fD1TanSPtB1 = new TH3F("DCAzTanSPtB1","DCAzTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);
123   fD1TanSPtB1->SetXTitle("tan(#theta)");
124   fD1TanSPtB1->SetYTitle("#sqrt(p_{t}(GeV/c))");
125   fD1TanSPtB1->SetZTitle("DCA_{z}");
126
127   fD0TanSPtL1 = new TH3F("DCAyTanSPtL1","DCAyTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);
128   fD0TanSPtL1->SetXTitle("tan(#theta)");
129   fD0TanSPtL1->SetYTitle("#sqrt{p_{t}(GeV/c)}");
130   fD0TanSPtL1->SetZTitle("DCA_{xy}");
131
132   fD1TanSPtL1 = new TH3F("DCAzTanSPtL1","DCAzTanSPt",40,-2,2, 10,0.3,3, 100, -1,1);
133   fD1TanSPtL1->SetXTitle("tan(#theta)");
134   fD1TanSPtL1->SetYTitle("#sqrt{p_{t}(GeV/c)}");
135   fD1TanSPtL1->SetZTitle("DCA_{z}");
136
137   fD0TanSPtInTPC = new TH3F("DCAyTanSPtInTPC","DCAyTanSPt",40,-2,2, 10,0.3,3, 100,-1,1);
138   fD0TanSPtInTPC->SetXTitle("tan(#theta)");
139   fD0TanSPtInTPC->SetYTitle("#sqrt{p_{t}(GeV/c)}");
140   fD0TanSPtInTPC->SetZTitle("DCA_{xy}");
141
142   fD1TanSPtInTPC = new TH3F("DCAzTanSPtInTPC","DCAzTanSPt",40,-2,2, 10,0.3,3, 100, -1,1);
143   fD1TanSPtInTPC->SetXTitle("tan(#theta)");
144   fD1TanSPtInTPC->SetYTitle("#sqrt{p_{t}(GeV/c)}");
145   fD1TanSPtInTPC->SetZTitle("DCA_{z}");
146
147   // init cuts
148   if(!fCutsMC) 
149     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
150   if(!fCutsRC) 
151     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
152  
153   // init folder
154   fAnalysisFolder = CreateFolder("folderDCA","Analysis DCA Folder");
155 }
156
157 //_____________________________________________________________________________
158 void AliComparisonDCA::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
159 {
160   // Fill DCA comparison information
161   AliExternalTrackParam *track = 0;
162   Double_t field      = AliTracker::GetBz(); // nominal Bz field [kG]
163   Double_t kMaxD      = 123456.0; // max distance
164
165   Int_t clusterITS[200];
166   Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
167
168   Float_t mcpt = infoMC->GetParticle().Pt();
169   Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
170   Float_t spt = TMath::Sqrt(mcpt);
171
172   // distance to Prim. vertex 
173   const Double_t* dv = infoMC->GetVDist(); 
174
175   Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
176
177   // Check selection cuts
178   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; 
179   if (!isPrim) return;
180   if (infoRC->GetStatus(1)!=3) return;
181   if (!infoRC->GetESDtrack()) return;  
182   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
183   if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
184
185   // calculate and set prim. vertex
186   fVertex->SetXv( infoMC->GetParticle().Vx() - dv[0] );
187   fVertex->SetYv( infoMC->GetParticle().Vy() - dv[1] );
188   fVertex->SetZv( infoMC->GetParticle().Vz() - dv[2] );
189
190   // calculate track parameters at vertex
191   if (infoRC->GetESDtrack()->GetTPCInnerParam())
192   {
193     if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )
194     {
195       Bool_t bDCAStatus = track->PropagateToDCA(fVertex,field,kMaxD,dca,cov);
196
197       if(bDCAStatus) {
198         fD0TanSPtInTPC->Fill(tantheta,spt,dca[0]);
199         fD1TanSPtInTPC->Fill(tantheta,spt,dca[1]);
200           }
201
202           delete track;
203     }
204   }
205
206  if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)==0){
207     fD0TanSPtB1->Fill(tantheta,spt,dca[0]);
208     fD1TanSPtB1->Fill(tantheta,spt,dca[1]);
209   }
210     fD0TanSPtL1->Fill(tantheta,spt,dca[0]);
211     fD1TanSPtL1->Fill(tantheta,spt,dca[1]);
212 }
213
214 //_____________________________________________________________________________
215 Long64_t AliComparisonDCA::Merge(TCollection* list) 
216 {
217   // Merge list of objects (needed by PROOF)
218
219   if (!list)
220   return 0;
221
222   if (list->IsEmpty())
223   return 1;
224
225   TIterator* iter = list->MakeIterator();
226   TObject* obj = 0;
227
228   // collection of generated histograms
229   Int_t count=0;
230   while((obj = iter->Next()) != 0) 
231   {
232     AliComparisonDCA* entry = dynamic_cast<AliComparisonDCA*>(obj);
233     if (entry == 0) continue; 
234     
235
236     fD0TanSPtB1->Add(entry->fD0TanSPtB1);
237     fD1TanSPtB1->Add(entry->fD1TanSPtB1);
238     fD0TanSPtL1->Add(entry->fD0TanSPtL1);
239     fD1TanSPtL1->Add(entry->fD1TanSPtL1);
240     fD0TanSPtInTPC->Add(entry->fD0TanSPtInTPC);
241     fD1TanSPtInTPC->Add(entry->fD1TanSPtInTPC);
242
243     count++;
244   }
245
246 return count;
247 }
248
249 //_____________________________________________________________________________
250 void AliComparisonDCA::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
251   // Process comparison information
252   Process(infoMC,infoRC);
253 }
254
255 //_____________________________________________________________________________
256 void AliComparisonDCA::Analyse()
257 {
258   //
259   // Analyse comparison information and store output histograms
260   // in the analysis folder "folderDCA" 
261   //
262   
263   TH1::AddDirectory(kFALSE);
264
265   TGraph2D *gr=0;
266   TGraph * gr0=0;
267   AliComparisonDCA * comp=this;
268   TFolder *folder = comp->GetAnalysisFolder();
269
270   // recreate folder every time
271   if(folder) delete folder;
272   folder = CreateFolder("folderDCA","Analysis DCA Folder");
273   folder->SetOwner();
274
275   // write results in the folder 
276   // Canvas to draw analysed histograms
277   TCanvas * c = new TCanvas("canDCA","DCA resolution");
278   c->Divide(1,2);
279   c->cd(1);
280   //
281   // DCA resolution
282   //
283   gr0 = AliMathBase::MakeStat1D(comp->fD0TanSPtB1,2,5);
284   gr0->GetXaxis()->SetTitle("Tan(#theta)");
285   gr0->GetYaxis()->SetTitle("#sigmaDCA (cm)");
286   gr0->SetName("DCAResolTan");
287   gr0->Draw("Al*");
288
289   //if(folder) folder->Add(gr0->GetHistogram());
290   if(folder) folder->Add(gr0);
291   //
292   c->cd(2);
293   gr = AliMathBase::MakeStat2D(comp->fD0TanSPtB1,4,2,5); 
294   gr->GetXaxis()->SetTitle("Tan(#theta)");
295   gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV/c)}");
296   gr->GetYaxis()->SetTitle("#sigmaDCA (cm)");
297   gr->SetName("DCAResolSPTTan");
298   gr->GetHistogram()->Draw("colz");
299
300   if(folder) folder->Add(gr->GetHistogram());
301
302   // set pointer to fAnalysisFolder
303   fAnalysisFolder = folder;
304 }
305
306 //_____________________________________________________________________________
307 TFolder* AliComparisonDCA::CreateFolder(TString name,TString title) { 
308 // create folder for analysed histograms
309 TFolder *folder = 0;
310   folder = new TFolder(name.Data(),title.Data());
311
312   return folder;
313 }