Adding abstract class for comparison components
[u/mrichter/AliRoot.git] / PWG1 / AliComparisonDEdx.cxx
1 //------------------------------------------------------------------------------
2 // Implementation of AliComparisonDEdx 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 which is 
8 // a data of AliComparisonDEdx.
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   AliComparisonDEdx * compObj = (AliComparisonDEdx*)f.Get("AliComparisonDEdx");
19
20   // analyse comparison data
21   compObj->Analyse();
22
23   // the output histograms/graphs will be stored in the folder "folderDEdx" 
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_DEdx.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 "AliESDEvent.h"
48 #include "AliESD.h"
49 #include "AliESDfriend.h"
50 #include "AliESDfriendTrack.h"
51 #include "AliRecInfoCuts.h" 
52 #include "AliMCInfoCuts.h" 
53 #include "AliLog.h" 
54 //
55 #include "AliMathBase.h"
56 #include "AliTreeDraw.h"
57 //#include "TStatToolkit.h"
58
59 #include "AliMCInfo.h" 
60 #include "AliESDRecInfo.h" 
61 #include "AliComparisonDEdx.h" 
62
63 using namespace std;
64
65 ClassImp(AliComparisonDEdx)
66
67 //_____________________________________________________________________________
68 AliComparisonDEdx::AliComparisonDEdx():
69 //  TNamed("AliComparisonDEdx","AliComparisonDEdx"),
70   AliComparisonObject("AliComparisonDEdx"),
71
72   // dEdx 
73   fTPCSignalNormTan(0), 
74   fTPCSignalNormSPhi(0),
75   fTPCSignalNormTPhi(0), 
76   //
77   fTPCSignalNormTanSPhi(0),
78   fTPCSignalNormTanTPhi(0),
79   fTPCSignalNormTanSPt(0), 
80   
81   // Cuts 
82   fCutsRC(0), 
83   fCutsMC(0),
84   fMCPtMin(0),
85   fMCAbsTanThetaMax(0),
86   fMCPdgCode(0),
87
88   // histogram folder 
89   fAnalysisFolder(0)
90 {
91   Init();
92 }
93
94 //_____________________________________________________________________________
95 AliComparisonDEdx::~AliComparisonDEdx(){
96    
97   if(fTPCSignalNormTan)  delete fTPCSignalNormTan; fTPCSignalNormTan=0; 
98   if(fTPCSignalNormSPhi) delete fTPCSignalNormSPhi; fTPCSignalNormSPhi=0;
99   if(fTPCSignalNormTPhi) delete fTPCSignalNormTPhi; fTPCSignalNormTPhi=0;
100   //
101   if(fTPCSignalNormTanSPhi) delete fTPCSignalNormTanSPhi; fTPCSignalNormTanSPhi=0;
102   if(fTPCSignalNormTanTPhi) delete fTPCSignalNormTanTPhi; fTPCSignalNormTanTPhi=0;
103   if(fTPCSignalNormTanSPt)  delete fTPCSignalNormTanSPt; fTPCSignalNormTanSPt=0;
104
105   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
106 }
107
108 //_____________________________________________________________________________
109 void AliComparisonDEdx::Init()
110 {
111   // Init histograms
112   
113   // TPC dEdx
114   fTPCSignalNormTan = new TH2F("CdEdxTan","CdEdxTan",50, -2,2,  40,30,70); 
115   fTPCSignalNormTan->SetXTitle("tan(#theta)");
116   fTPCSignalNormTan->SetYTitle("rec. dE/dx / calc. dE/dx");
117
118   fTPCSignalNormSPhi   = new TH2F("CdEdxSPhi","CdEdxSPhi",10,0.0,1,40,30,70);
119   fTPCSignalNormSPhi->SetXTitle("sin(#phi)");
120   fTPCSignalNormSPhi->SetYTitle("rec. dE/dx / calc. dE/dx");
121
122   fTPCSignalNormTPhi   = new TH2F("CdEdxTPhi","CdEdxTPhi",10,0.0,2,40,30,70); 
123   fTPCSignalNormTPhi->SetXTitle("tan(#phi)");
124   fTPCSignalNormTPhi->SetYTitle("rec. dE/dx / calc. dE/dx");
125
126   fTPCSignalNormTanSPhi= new TH3F("CdEdxTanSPhi","CdEdxTanSPhi",20, -2,2, 10,0.0 ,1,  40,30,70);
127   fTPCSignalNormTanSPhi->SetXTitle("tan(#theta)");
128   fTPCSignalNormTanSPhi->SetYTitle("sin(#phi)");
129   fTPCSignalNormTanSPhi->SetZTitle("rec. dE/dx / calc. dE/dx");
130
131   fTPCSignalNormTanTPhi= new TH3F("CdEdxTanTPhi","CdEdxTanTPhi",20, -2,2, 10,0.0 ,1,  40,30,70);
132   fTPCSignalNormTanTPhi->SetXTitle("tan(#theta)");
133   fTPCSignalNormTanTPhi->SetYTitle("tan(#phi)");
134   fTPCSignalNormTanTPhi->SetZTitle("rec. dE/dx / calc. dE/dx");
135
136   fTPCSignalNormTanSPt= new TH3F("CdEdxTanSPt","CdEdxTanSPt",20, -2,2, 10,0.3 ,3, 40,30,70); 
137   fTPCSignalNormTanSPt->SetXTitle("tan(#theta)");
138   fTPCSignalNormTanSPt->SetYTitle("#sqrt{p_{t}}");
139   fTPCSignalNormTanSPt->SetZTitle("rec. dE/dx / calc. dE/dx");
140
141   // Init cuts
142   if(!fCutsMC) 
143     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
144   if(!fCutsRC) 
145     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
146
147     // init folder
148     fAnalysisFolder = CreateFolder("folderDEdx","Analysis de/dx Folder");
149 }
150
151 //_____________________________________________________________________________
152 void AliComparisonDEdx::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
153
154   // Fill dE/dx  comparison information
155   
156   Float_t mcpt = infoMC->GetParticle().Pt();
157   Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
158   Float_t mprim = infoMC->GetPrim();
159
160   // distance to Prim. vertex 
161   const Double_t* dv = infoMC->GetVDist(); 
162
163   Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
164   
165   // Check selection cuts 
166   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; 
167   if (!isPrim) return;
168   if (infoRC->GetStatus(1)!=3) return;
169   if (!infoRC->GetESDtrack()) return;  
170   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
171   if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
172   //if (mprim>1.4) return;
173   //if (mprim<0.5) return;
174   if (mprim > fCutsMC->GetMaxTPCSignal()) return;
175   if (mprim < fCutsMC->GetMinTPCSignal()) return;
176   if (infoRC->GetESDtrack()->GetTPCsignalN()<fCutsRC->GetMinTPCsignalN()) return;
177   //
178   Float_t ratio = infoRC->GetESDtrack()->GetTPCsignal()/infoMC->GetPrim();
179   Float_t sphi =  infoRC->GetESDtrack()->GetInnerParam()->GetSnp();
180   Float_t tphi =  sphi/TMath::Sqrt(1-sphi*sphi);
181
182   if (TMath::Abs(infoMC->GetParticle().GetPdgCode()) != GetMCPdgCode()) return;
183   //if (mcpt>0.5) {
184   if (mcpt > GetMCPtMin()) {
185     fTPCSignalNormTan->Fill(tantheta,ratio);    // only subset
186   }
187
188   //if (TMath::Abs(tantheta)<0.5){
189   if (TMath::Abs(tantheta) < GetMCAbsTanThetaMax()){
190     fTPCSignalNormSPhi->Fill(sphi,ratio);       // only subset
191     fTPCSignalNormTPhi->Fill(tphi,ratio);       // only subset
192   }
193   fTPCSignalNormTanSPhi->Fill(tantheta,sphi,ratio);    
194   fTPCSignalNormTanTPhi->Fill(tantheta,tphi,ratio);    
195   fTPCSignalNormTanSPt->Fill(tantheta,TMath::Sqrt(mcpt),ratio);
196 }
197
198 //_____________________________________________________________________________
199 Long64_t AliComparisonDEdx::Merge(TCollection* list) 
200 {
201   // Merge list of objects (needed by PROOF)
202
203   if (!list)
204   return 0;
205
206   if (list->IsEmpty())
207   return 1;
208
209   TIterator* iter = list->MakeIterator();
210   TObject* obj = 0;
211
212   // collection of generated histograms
213   Int_t count=0;
214   while((obj = iter->Next()) != 0) 
215   {
216     AliComparisonDEdx* entry = dynamic_cast<AliComparisonDEdx*>(obj);
217     if (entry == 0) continue;
218
219     fTPCSignalNormTan->Add(entry->fTPCSignalNormTan);
220     fTPCSignalNormSPhi->Add(entry->fTPCSignalNormSPhi);
221     fTPCSignalNormTPhi->Add(entry->fTPCSignalNormTPhi);
222     //
223     fTPCSignalNormTanSPhi->Add(entry->fTPCSignalNormTanSPhi);
224     fTPCSignalNormTanTPhi->Add(entry->fTPCSignalNormTanTPhi);
225     fTPCSignalNormTanSPt->Add(entry->fTPCSignalNormTanSPt);
226
227     count++;
228   }
229
230 return count;
231 }
232
233 //_____________________________________________________________________________
234 void AliComparisonDEdx::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
235 {
236   // Process comparison information
237   Process(infoMC,infoRC);
238 }
239
240 //_____________________________________________________________________________
241 TH1F* AliComparisonDEdx::MakeResol(TH2F * his, Int_t integ, Bool_t type)
242 {
243   // Make resolution histograms
244   TH1F *hisr, *hism;
245   if (!gPad) new TCanvas;
246   hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);
247   if (type) return hism;
248   else 
249     return hisr;
250 }
251
252 //_____________________________________________________________________________
253 void AliComparisonDEdx::Analyse()
254 {
255   // Analyse comparison information and store output histograms
256   // in the folder "folderDEdx"
257   //
258
259   TH1::AddDirectory(kFALSE);
260   
261   AliComparisonDEdx * comp=this;
262   TFolder *folder = comp->GetAnalysisFolder();
263
264   TH1F *hiss=0;
265   TGraph2D * gr=0;
266
267   // recreate folder every time
268   if(folder) delete folder;
269   folder = CreateFolder("folderDEdx","Analysis DEdx Folder");
270   folder->SetOwner();
271
272   // write results in the folder 
273   TCanvas * c = new TCanvas("can","TPC dedx");
274   c->cd();
275
276   hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,0);
277   hiss->SetXTitle("Tan(#theta)");
278   hiss->SetYTitle("#sigma_{dEdx}");
279   hiss->Draw();
280   hiss->SetName("TPCdEdxResolTan");
281
282   if(folder) folder->Add(hiss);
283   //
284   hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,1); 
285   hiss->SetXTitle("Tan(#theta)");
286   hiss->SetYTitle("<dEdx>");
287   hiss->Draw(); 
288   hiss->SetName("TPCdEdxMeanTan");
289
290   if(folder) folder->Add(hiss);
291   //
292   gr = AliMathBase::MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,4);
293   gr->GetXaxis()->SetTitle("Tan(#theta)");
294   gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
295   gr->GetZaxis()->SetTitle("<dEdx>");
296   gr->SetName("TPCdEdxMeanTanPt_1");
297   gr->GetHistogram()->Draw("colz"); 
298
299   if(folder) folder->Add(gr->GetHistogram());
300   //
301   gr = AliMathBase::MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,5);
302   gr->GetXaxis()->SetTitle("Tan(#theta)");
303   gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
304   gr->GetZaxis()->SetTitle("#sigma_{dEdx}");
305   gr->SetName("TPCdEdxMeanTanPt_2");
306   gr->GetHistogram()->Draw("colz"); 
307
308   if(folder) folder->Add(gr->GetHistogram());
309
310   // set pointer to fAnalysisFolder
311   fAnalysisFolder = folder;
312 }
313
314 //_____________________________________________________________________________
315 TFolder* AliComparisonDEdx::CreateFolder(TString name,TString title) { 
316 // create folder for analysed histograms
317 TFolder *folder = 0;
318   folder = new TFolder(name.Data(),title.Data());
319
320   return folder;
321 }
322