175079fe6f4033b3438f26cf42d167547c041d73
[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   gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C");
17   LoadMyLibs();
18   TFile f("Output.root");
19   //AliComparisonDEdx * compObj = (AliComparisonDEdx*)f.Get("AliComparisonDEdx");
20   AliComparisonDEdx * compObj = (AliComparisonDEdx*)cOutput->FindObject("AliComparisonDEdx");
21
22   // Analyse comparison data
23   compObj->Analyse();
24
25   // the output histograms/graphs will be stored in the folder "folderDEdx" 
26   compObj->GetAnalysisFolder()->ls("*");
27
28   // user can save whole comparison object (or only folder with anlysed histograms) 
29   // in the seperate output file (e.g.)
30   TFile fout("Analysed_DEdx.root"."recreate");
31   compObj->Write(); // compObj->GetAnalysisFolder()->Write();
32   fout.Close();
33
34 */
35
36 #include <iostream>
37
38 #include "TFile.h"
39 #include "TCint.h"
40 #include "TH3F.h"
41 #include "TH2F.h"
42 #include "TF1.h"
43 #include "TProfile.h"
44 #include "TProfile2D.h"
45 #include "TGraph2D.h"
46 #include "TCanvas.h"
47 #include "TGraph.h"
48 //
49 #include "AliESDEvent.h"
50 #include "AliESD.h"
51 #include "AliESDfriend.h"
52 #include "AliESDfriendTrack.h"
53 #include "AliRecInfoCuts.h" 
54 #include "AliMCInfoCuts.h" 
55 #include "AliLog.h" 
56 //
57 #include "AliMathBase.h"
58 #include "AliTreeDraw.h"
59 //#include "TStatToolkit.h"
60
61 #include "AliMCInfo.h" 
62 #include "AliESDRecInfo.h" 
63 #include "AliComparisonDEdx.h" 
64
65 using namespace std;
66
67 ClassImp(AliComparisonDEdx)
68
69 //_____________________________________________________________________________
70 AliComparisonDEdx::AliComparisonDEdx():
71 //  TNamed("AliComparisonDEdx","AliComparisonDEdx"),
72   AliComparisonObject("AliComparisonDEdx"),
73
74   // dEdx 
75   fTPCSignalNormTan(0), 
76   fTPCSignalNormSPhi(0),
77   fTPCSignalNormTPhi(0), 
78   //
79   fTPCSignalNormTanSPhi(0),
80   fTPCSignalNormTanTPhi(0),
81   fTPCSignalNormTanSPt(0), 
82   
83   // Cuts 
84   fCutsRC(0), 
85   fCutsMC(0),
86   fMCPtMin(0),
87   fMCAbsTanThetaMax(0),
88   fMCPdgCode(0),
89
90   // histogram folder 
91   fAnalysisFolder(0)
92 {
93   Init();
94 }
95
96 //_____________________________________________________________________________
97 AliComparisonDEdx::~AliComparisonDEdx(){
98    
99   if(fTPCSignalNormTan)  delete fTPCSignalNormTan; fTPCSignalNormTan=0; 
100   if(fTPCSignalNormSPhi) delete fTPCSignalNormSPhi; fTPCSignalNormSPhi=0;
101   if(fTPCSignalNormTPhi) delete fTPCSignalNormTPhi; fTPCSignalNormTPhi=0;
102   //
103   if(fTPCSignalNormTanSPhi) delete fTPCSignalNormTanSPhi; fTPCSignalNormTanSPhi=0;
104   if(fTPCSignalNormTanTPhi) delete fTPCSignalNormTanTPhi; fTPCSignalNormTanTPhi=0;
105   if(fTPCSignalNormTanSPt)  delete fTPCSignalNormTanSPt; fTPCSignalNormTanSPt=0;
106
107   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
108 }
109
110 //_____________________________________________________________________________
111 void AliComparisonDEdx::Init()
112 {
113   // Init histograms
114   
115   // TPC dEdx
116   fTPCSignalNormTan = new TH2F("CdEdxTan","CdEdxTan",50, -2,2,  40,30,70); 
117   fTPCSignalNormTan->SetXTitle("tan(#theta)");
118   fTPCSignalNormTan->SetYTitle("rec. dE/dx / calc. dE/dx");
119
120   fTPCSignalNormSPhi   = new TH2F("CdEdxSPhi","CdEdxSPhi",10,0.0,1,40,30,70);
121   fTPCSignalNormSPhi->SetXTitle("sin(#phi)");
122   fTPCSignalNormSPhi->SetYTitle("rec. dE/dx / calc. dE/dx");
123
124   fTPCSignalNormTPhi   = new TH2F("CdEdxTPhi","CdEdxTPhi",10,0.0,2,40,30,70); 
125   fTPCSignalNormTPhi->SetXTitle("tan(#phi)");
126   fTPCSignalNormTPhi->SetYTitle("rec. dE/dx / calc. dE/dx");
127
128   fTPCSignalNormTanSPhi= new TH3F("CdEdxTanSPhi","CdEdxTanSPhi",20, -2,2, 10,0.0 ,1,  40,30,70);
129   fTPCSignalNormTanSPhi->SetXTitle("tan(#theta)");
130   fTPCSignalNormTanSPhi->SetYTitle("sin(#phi)");
131   fTPCSignalNormTanSPhi->SetZTitle("rec. dE/dx / calc. dE/dx");
132
133   fTPCSignalNormTanTPhi= new TH3F("CdEdxTanTPhi","CdEdxTanTPhi",20, -2,2, 10,0.0 ,1,  40,30,70);
134   fTPCSignalNormTanTPhi->SetXTitle("tan(#theta)");
135   fTPCSignalNormTanTPhi->SetYTitle("tan(#phi)");
136   fTPCSignalNormTanTPhi->SetZTitle("rec. dE/dx / calc. dE/dx");
137
138   fTPCSignalNormTanSPt= new TH3F("CdEdxTanSPt","CdEdxTanSPt",20, -2,2, 10,0.3 ,3, 40,30,70); 
139   fTPCSignalNormTanSPt->SetXTitle("tan(#theta)");
140   fTPCSignalNormTanSPt->SetYTitle("#sqrt{p_{t}}");
141   fTPCSignalNormTanSPt->SetZTitle("rec. dE/dx / calc. dE/dx");
142
143   // Init cuts
144   if(!fCutsMC) 
145     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
146   if(!fCutsRC) 
147     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
148
149     // init folder
150     fAnalysisFolder = CreateFolder("folderDEdx","Analysis de/dx Folder");
151 }
152
153 //_____________________________________________________________________________
154 void AliComparisonDEdx::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
155
156   // Fill dE/dx  comparison information
157   
158   Float_t mcpt = infoMC->GetParticle().Pt();
159   Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
160   Float_t mprim = infoMC->GetPrim();
161
162   // distance to Prim. vertex 
163   const Double_t* dv = infoMC->GetVDist(); 
164
165   Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
166   
167   // Check selection cuts 
168   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; 
169   if (!isPrim) return;
170   if (infoRC->GetStatus(1)!=3) return;
171   if (!infoRC->GetESDtrack()) return;  
172   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
173   if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
174   //if (mprim>1.4) return;
175   //if (mprim<0.5) return;
176   if (mprim > fCutsMC->GetMaxTPCSignal()) return;
177   if (mprim < fCutsMC->GetMinTPCSignal()) return;
178   if (infoRC->GetESDtrack()->GetTPCsignalN()<fCutsRC->GetMinTPCsignalN()) return;
179   //
180   Float_t ratio = infoRC->GetESDtrack()->GetTPCsignal()/infoMC->GetPrim();
181   Float_t sphi =  infoRC->GetESDtrack()->GetInnerParam()->GetSnp();
182   Float_t tphi =  sphi/TMath::Sqrt((1.-sphi)*(1.+sphi));
183
184   if (TMath::Abs(infoMC->GetParticle().GetPdgCode()) != GetMCPdgCode()) return;
185   //if (mcpt>0.5) {
186   if (mcpt > GetMCPtMin()) {
187     fTPCSignalNormTan->Fill(tantheta,ratio);    // only subset
188   }
189
190   //if (TMath::Abs(tantheta)<0.5){
191   if (TMath::Abs(tantheta) < GetMCAbsTanThetaMax()){
192     fTPCSignalNormSPhi->Fill(sphi,ratio);       // only subset
193     fTPCSignalNormTPhi->Fill(tphi,ratio);       // only subset
194   }
195   fTPCSignalNormTanSPhi->Fill(tantheta,sphi,ratio);    
196   fTPCSignalNormTanTPhi->Fill(tantheta,tphi,ratio);    
197   fTPCSignalNormTanSPt->Fill(tantheta,TMath::Sqrt(mcpt),ratio);
198 }
199
200 //_____________________________________________________________________________
201 Long64_t AliComparisonDEdx::Merge(TCollection* list) 
202 {
203   // Merge list of objects (needed by PROOF)
204
205   if (!list)
206   return 0;
207
208   if (list->IsEmpty())
209   return 1;
210
211   TIterator* iter = list->MakeIterator();
212   TObject* obj = 0;
213
214   // collection of generated histograms
215   Int_t count=0;
216   while((obj = iter->Next()) != 0) 
217   {
218     AliComparisonDEdx* entry = dynamic_cast<AliComparisonDEdx*>(obj);
219     if (entry == 0) continue;
220
221     fTPCSignalNormTan->Add(entry->fTPCSignalNormTan);
222     fTPCSignalNormSPhi->Add(entry->fTPCSignalNormSPhi);
223     fTPCSignalNormTPhi->Add(entry->fTPCSignalNormTPhi);
224     //
225     fTPCSignalNormTanSPhi->Add(entry->fTPCSignalNormTanSPhi);
226     fTPCSignalNormTanTPhi->Add(entry->fTPCSignalNormTanTPhi);
227     fTPCSignalNormTanSPt->Add(entry->fTPCSignalNormTanSPt);
228
229     count++;
230   }
231
232 return count;
233 }
234
235 //_____________________________________________________________________________
236 void AliComparisonDEdx::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
237 {
238   // Process comparison information
239   Process(infoMC,infoRC);
240 }
241
242 //_____________________________________________________________________________
243 TH1F* AliComparisonDEdx::MakeResol(TH2F * his, Int_t integ, Bool_t type)
244 {
245   // Make resolution histograms
246   TH1F *hisr, *hism;
247   if (!gPad) new TCanvas;
248   hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);
249   if (type) return hism;
250   else 
251     return hisr;
252 }
253
254 //_____________________________________________________________________________
255 void AliComparisonDEdx::Analyse()
256 {
257   // Analyze comparison information and store output histograms
258   // in the folder "folderDEdx"
259   //
260
261   TH1::AddDirectory(kFALSE);
262   
263   AliComparisonDEdx * comp=this;
264   TObjArray *aFolderObj = new TObjArray;
265
266   TH1F *hiss=0;
267   TGraph2D * gr=0;
268
269   // write results in the folder 
270   TCanvas * c = new TCanvas("can","TPC dedx");
271   c->cd();
272
273   hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,0);
274   hiss->SetXTitle("Tan(#theta)");
275   hiss->SetYTitle("#sigma_{dEdx}");
276   hiss->Draw();
277   hiss->SetName("TPCdEdxResolTan");
278
279   aFolderObj->Add(hiss);
280   //
281   hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,1); 
282   hiss->SetXTitle("Tan(#theta)");
283   hiss->SetYTitle("<dEdx>");
284   hiss->Draw(); 
285   hiss->SetName("TPCdEdxMeanTan");
286
287   aFolderObj->Add(hiss);
288   //
289   gr = AliMathBase::MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,4);
290   gr->GetXaxis()->SetTitle("Tan(#theta)");
291   gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
292   gr->GetZaxis()->SetTitle("<dEdx>");
293   gr->SetName("TPCdEdxMeanTanPt_1");
294   gr->GetHistogram()->Draw("colz"); 
295
296   gr->GetHistogram()->SetName("TPCdEdxMeanTanPt_1");
297   aFolderObj->Add(gr->GetHistogram());
298   //
299   gr = AliMathBase::MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,5);
300   gr->GetXaxis()->SetTitle("Tan(#theta)");
301   gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
302   gr->GetZaxis()->SetTitle("#sigma_{dEdx}");
303   gr->SetName("TPCdEdxMeanTanPt_2");
304   gr->GetHistogram()->Draw("colz"); 
305
306   gr->GetHistogram()->SetName("TPCdEdxMeanTanPt_2");
307   aFolderObj->Add(gr->GetHistogram());
308
309   // export objects to analysis folder
310   fAnalysisFolder = ExportToFolder(aFolderObj);
311
312   // delete only TObjrArray
313   if(aFolderObj) delete aFolderObj;
314 }
315
316 //_____________________________________________________________________________
317 TFolder* AliComparisonDEdx::ExportToFolder(TObjArray * array) 
318 {
319   // recreate folder avery time and export objects to new one
320   //
321   AliComparisonDEdx * comp=this;
322   TFolder *folder = comp->GetAnalysisFolder();
323
324   TString name, title;
325   TFolder *newFolder = 0;
326   Int_t i = 0;
327   Int_t size = array->GetSize();
328
329   if(folder) { 
330      // get name and title from old folder
331      name = folder->GetName();  
332      title = folder->GetTitle();  
333
334          // delete old one
335      delete folder;
336
337          // create new one
338      newFolder = CreateFolder(name.Data(),title.Data());
339      newFolder->SetOwner();
340
341          // add objects to folder
342      while(i < size) {
343            newFolder->Add(array->At(i));
344            i++;
345          }
346   }
347
348 return newFolder;
349 }
350
351
352 //_____________________________________________________________________________
353 TFolder* AliComparisonDEdx::CreateFolder(TString name,TString title) { 
354 // create folder for analysed histograms
355 TFolder *folder = 0;
356   folder = new TFolder(name.Data(),title.Data());
357
358   return folder;
359 }
360