When calculating a*a-b*b the form (a-b)*(a+b) is usually more numerically stable.
[u/mrichter/AliRoot.git] / PWG1 / AliComparisonDEdx.cxx
CommitLineData
3baa4bfd 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
b4126c69 16 gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C");
17 LoadMyLibs();
3baa4bfd 18 TFile f("Output.root");
19 AliComparisonDEdx * compObj = (AliComparisonDEdx*)f.Get("AliComparisonDEdx");
20
b4126c69 21 // Analyse comparison data
3baa4bfd 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
64using namespace std;
65
66ClassImp(AliComparisonDEdx)
67
68//_____________________________________________________________________________
69AliComparisonDEdx::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//_____________________________________________________________________________
96AliComparisonDEdx::~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//_____________________________________________________________________________
110void 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//_____________________________________________________________________________
153void 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();
60e55aee 181 Float_t tphi = sphi/TMath::Sqrt((1.-sphi)*(1.+sphi));
3baa4bfd 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//_____________________________________________________________________________
200Long64_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
231return count;
232}
233
234//_____________________________________________________________________________
235void AliComparisonDEdx::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
236{
237 // Process comparison information
238 Process(infoMC,infoRC);
239}
240
241//_____________________________________________________________________________
242TH1F* 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//_____________________________________________________________________________
254void AliComparisonDEdx::Analyse()
255{
b4126c69 256 // Analyze comparison information and store output histograms
3baa4bfd 257 // in the folder "folderDEdx"
258 //
259
260 TH1::AddDirectory(kFALSE);
261
262 AliComparisonDEdx * comp=this;
b4126c69 263 TObjArray *aFolderObj = new TObjArray;
3baa4bfd 264
265 TH1F *hiss=0;
266 TGraph2D * gr=0;
267
3baa4bfd 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
b4126c69 278 aFolderObj->Add(hiss);
3baa4bfd 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
b4126c69 286 aFolderObj->Add(hiss);
3baa4bfd 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
decf0997 295 gr->GetHistogram()->SetName("TPCdEdxMeanTanPt_1");
b4126c69 296 aFolderObj->Add(gr->GetHistogram());
3baa4bfd 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
decf0997 305 gr->GetHistogram()->SetName("TPCdEdxMeanTanPt_2");
b4126c69 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//_____________________________________________________________________________
316TFolder* 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;
3baa4bfd 335
b4126c69 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
347return newFolder;
3baa4bfd 348}
349
b4126c69 350
3baa4bfd 351//_____________________________________________________________________________
352TFolder* AliComparisonDEdx::CreateFolder(TString name,TString title) {
353// create folder for analysed histograms
354TFolder *folder = 0;
355 folder = new TFolder(name.Data(),title.Data());
356
357 return folder;
358}
359