]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/AliComparisonDCA.cxx
Histogram ranges changed to cut off saturation peak and noise
[u/mrichter/AliRoot.git] / PWG1 / AliComparisonDCA.cxx
CommitLineData
09b20ad1 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) are stored in the output picture_dca.root file.
8//
9// Author: J.Otwinowski 04/02/2008
10//------------------------------------------------------------------------------
11
12/*
13 // after running analysis, read the file, and get component
14 gSystem->Load("libPWG1.so");
15 TFile f("Output.root");
16 AliComparisonDCA * comp = (AliComparisonDCA*)f.Get("AliComparisonDCA");
17
18 // analyse comparison data (output stored in pictures_dca.root)
19 comp->Analyse();
20
21 // TPC track length parameterisation
22 TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2); // length function
23 TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);
24 fl2.SetParameter(1,1);
25 fl2.SetParameter(0,1);
26*/
27
28#include <iostream>
29
30#include "TFile.h"
31#include "TCint.h"
32#include "TH3F.h"
33#include "TH2F.h"
34#include "TF1.h"
35#include "TProfile.h"
36#include "TProfile2D.h"
37#include "TGraph2D.h"
38#include "TCanvas.h"
39#include "TGraph.h"
40//
41#include "AliTracker.h"
42#include "AliESDEvent.h"
43#include "AliESD.h"
44#include "AliESDfriend.h"
45#include "AliESDfriendTrack.h"
46#include "AliRecInfoCuts.h"
47#include "AliMCInfoCuts.h"
48#include "AliLog.h"
49#include "AliESDVertex.h"
50//
51#include "AliMathBase.h"
52#include "AliTreeDraw.h"
53
54#include "AliMCInfo.h"
55#include "AliESDRecInfo.h"
56#include "AliComparisonDCA.h"
57
58using namespace std;
59
60ClassImp(AliComparisonDCA)
61
62//_____________________________________________________________________________
63AliComparisonDCA::AliComparisonDCA():
64 TNamed("AliComparisonDCA","AliComparisonDCA"),
65
66 // DCA histograms
67 fD0TanSPtB1(0),
68 fD1TanSPtB1(0),
69 fD0TanSPtL1(0),
70 fD1TanSPtL1(0),
71 fD0TanSPtInTPC(0),
72 fD1TanSPtInTPC(0),
73 fVertex(0),
74
75 // Cuts
76 fCutsRC(0),
77 fCutsMC(0)
78{
79 InitHisto();
80 InitCuts();
81
82 // vertex (0,0,0)
83 fVertex = new AliESDVertex();
84 fVertex->SetXv(0.0);
85 fVertex->SetYv(0.0);
86 fVertex->SetZv(0.0);
87}
88
89//_____________________________________________________________________________
90AliComparisonDCA::~AliComparisonDCA()
91{
92 //
93 if(fD0TanSPtB1) delete fD0TanSPtB1; fD0TanSPtB1=0;
94 if(fD1TanSPtB1) delete fD1TanSPtB1; fD1TanSPtB1=0;
95 if(fD0TanSPtL1) delete fD0TanSPtL1; fD0TanSPtL1=0;
96 if(fD1TanSPtL1) delete fD1TanSPtL1; fD1TanSPtL1=0;
97 if(fD0TanSPtInTPC) delete fD0TanSPtInTPC; fD0TanSPtInTPC=0;
98 if(fD1TanSPtInTPC) delete fD1TanSPtInTPC; fD1TanSPtInTPC=0;
99 if(fVertex) delete fVertex; fVertex=0;
100}
101
102//_____________________________________________________________________________
103void AliComparisonDCA::InitHisto()
104{
105 // DCA histograms
106 fD0TanSPtB1 = new TH3F("DCAyTanSPtB1","DCAyTanSPt",20,1,2, 10,0.3,2, 100,-4,4);
107 fD0TanSPtB1->SetXTitle("tan(#theta)");
108 fD0TanSPtB1->SetYTitle("sqrt(p_{t})");
109 fD0TanSPtB1->SetZTitle("DCA_{xy}");
110
111 fD1TanSPtB1 = new TH3F("DCAzTanSPtB1","DCAzTanSPt",20,1,2, 10,0.3,2, 100,-4,4);
112 fD1TanSPtB1->SetXTitle("tan(#theta)");
113 fD1TanSPtB1->SetYTitle("sqrt(p_{t})");
114 fD1TanSPtB1->SetZTitle("DCA_{z}");
115
116 fD0TanSPtL1 = new TH3F("DCAyTanSPtL1","DCAyTanSPt",20,0,1, 10,0.3,2, 100,-0.1,0.1);
117 fD0TanSPtL1->SetXTitle("tan(#theta)");
118 fD0TanSPtL1->SetYTitle("sqrt(p_{t})");
119 fD0TanSPtL1->SetZTitle("DCA_{xy}");
120
121 fD1TanSPtL1 = new TH3F("DCAzTanSPtL1","DCAzTanSPt",20,0,1, 10,0.3,2, 100, -0.1,0.1);
122 fD1TanSPtL1->SetXTitle("tan(#theta)");
123 fD1TanSPtL1->SetYTitle("sqrt(p_{t})");
124 fD1TanSPtL1->SetZTitle("DCA_{z}");
125
126 fD0TanSPtInTPC = new TH3F("DCAyTanSPtInTPC","DCAyTanSPt",100,0,100, 10,0.3,2, 100,-0.1,0.1);
127 fD0TanSPtInTPC->SetXTitle("tan(#theta)");
128 fD0TanSPtInTPC->SetYTitle("sqrt(p_{t})");
129 fD0TanSPtInTPC->SetZTitle("DCA_{xy}");
130
131 fD1TanSPtInTPC = new TH3F("DCAzTanSPtInTPC","DCAzTanSPt",100,0,100, 10,0.3,2, 100, -0.1,0.1);
132 fD1TanSPtInTPC->SetXTitle("tan(#theta)");
133 fD1TanSPtInTPC->SetYTitle("sqrt(p_{t})");
134 fD1TanSPtInTPC->SetZTitle("DCA_{z}");
135}
136
137//_____________________________________________________________________________
138void AliComparisonDCA::InitCuts()
139{
140 if(!fCutsMC)
141 AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
142 if(!fCutsRC)
143 AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
144}
145
146//_____________________________________________________________________________
147void AliComparisonDCA::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC)
148{
149 // Fill DCA comparison information
150
151 Float_t mcpt = infoMC->GetParticle().Pt();
152 Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
153 Bool_t isPrim = infoMC->GetParticle().R()<fCutsMC->GetMaxR() && TMath::Abs(infoMC->GetParticle().Vz())<fCutsMC->GetMaxVz() ;
154 Float_t spt = TMath::Sqrt(mcpt);
155
156 AliExternalTrackParam *track = 0;
157 Double_t kRadius = 3.0; // beam pipe radius
158 Double_t kMaxStep = 5.0; // max step
159 Double_t field = 0.4; // mag. field
160 Double_t kMaxD = 123456.0; // max distance
161
162 Int_t clusterITS[200];
163 Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
164 const Double_t* dv;
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
174 // calculate and set prim. vertex
175 dv = infoMC->GetVDist(); // distance to Prim. vertex
176
177 fVertex->SetXv( infoMC->GetParticle().Vx() - dv[0] );
178 fVertex->SetYv( infoMC->GetParticle().Vy() - dv[1] );
179 fVertex->SetZv( infoMC->GetParticle().Vz() - dv[2] );
180
181 // calculate track parameters at vertex
182
183 if (infoRC->GetESDtrack()->GetTPCInnerParam())
184 {
185 if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )
186 {
187 AliTracker::PropagateTrackTo(track,kRadius,infoMC->GetMass(),kMaxStep,kTRUE);
188 track->PropagateToDCA(fVertex,field,kMaxD,dca,cov);
189
190 fD0TanSPtInTPC->Fill(tantheta,spt,dca[0]);
191 fD1TanSPtInTPC->Fill(tantheta,spt,dca[1]);
192
193 delete track;
194 }
195 }
196
197 if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)==0){
198 fD0TanSPtB1->Fill(tantheta,spt,dca[0]);
199 fD1TanSPtB1->Fill(tantheta,spt,dca[1]);
200 }
201 fD0TanSPtL1->Fill(tantheta,spt,dca[0]);
202 fD1TanSPtL1->Fill(tantheta,spt,dca[1]);
203}
204
205//_____________________________________________________________________________
206Long64_t AliComparisonDCA::Merge(TCollection* list)
207{
208 // Merge list of objects (needed by PROOF)
209
210 if (!list)
211 return 0;
212
213 if (list->IsEmpty())
214 return 1;
215
216 TIterator* iter = list->MakeIterator();
217 TObject* obj = 0;
218
219 // collection of generated histograms
220 Int_t count=0;
221 while((obj = iter->Next()) != 0)
222 {
223 AliComparisonDCA* entry = dynamic_cast<AliComparisonDCA*>(obj);
224 if (entry == 0) {
225 Error("Add","Attempt to add object of class: %s to a %s",
226 obj->ClassName(),this->ClassName());
227 return -1;
228 }
229
230 fD0TanSPtB1->Add(entry->fD0TanSPtB1);
231 fD1TanSPtB1->Add(entry->fD1TanSPtB1);
232 fD0TanSPtL1->Add(entry->fD0TanSPtL1);
233 fD1TanSPtL1->Add(entry->fD1TanSPtL1);
234 fD0TanSPtInTPC->Add(entry->fD0TanSPtInTPC);
235 fD1TanSPtInTPC->Add(entry->fD1TanSPtInTPC);
236
237 count++;
238 }
239
240return count;
241}
242
243//_____________________________________________________________________________
244void AliComparisonDCA::Exec(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
245 // Process comparison information
246 Process(infoMC,infoRC);
247}
248
249//_____________________________________________________________________________
250void AliComparisonDCA::Analyse()
251{
252 // Analyse output histograms
253
254 AliComparisonDCA * comp=this;
255 TGraph2D * gr=0;
256 TGraph * gr0=0;
257
258 TFile *fp = new TFile("pictures_dca.root","recreate");
259 fp->cd();
260
261 TCanvas * c = new TCanvas("DCA","DCA resloution");
262 c->cd();
263
264 // DCA resolution
265 gr0 = AliMathBase::MakeStat1D(comp->fD0TanSPtB1,2,5);
266 gr0->GetXaxis()->SetTitle("Tan(#theta)");
267 gr0->GetYaxis()->SetTitle("#sigmaDCA (cm)");
268 gr0->Write("DCAResolTan");
269 //
270 gr = AliMathBase::MakeStat2D(comp->fD0TanSPtB1,4,2,5);
271 gr->GetXaxis()->SetTitle("Tan(#theta)");
272 gr->GetYaxis()->SetTitle("#sigmaDCA (cm)");
273 gr->GetHistogram()->Write("DCAResolSPTTan");
274
275 gPad->Clear();
276 gr0->Draw("al*");
277
278 fp->Close();
279}