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