]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/TPCupgrade/macros/fastToyMCMI.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
TPC module
[u/mrichter/AliRoot.git] / TPC / TPCupgrade / macros / fastToyMCMI.C
CommitLineData
d9502b8d 1/*
c6b97cea 2 Fast toy MC for different purposes. primary goal prepare the motivation figures for the TPC TDR and
3 expalantion internal note.
4
5 testExtrapolationError() - to show track extrapolation errror
6 testdEdxGEM(); -to show the dEdx respolution detoriation as function of the electron trasnparency
d9502b8d 7
c6b97cea 8
e1ed5ec5 9 .x $HOME/NimStyle.C
d9502b8d 10 .L $ALICE_ROOT/TPC/Upgrade/macros/fastToyMCMI.C+
e1ed5ec5 11
d9502b8d 12*/
2dd02504 13#include "TStyle.h"
14#include "TCut.h"
15#include "TCanvas.h"
d9502b8d 16#include "TMath.h"
17#include "TVectorD.h"
18#include "TLinearFitter.h"
19#include "TGeoGlobalMagField.h"
20#include "TRandom.h"
21#include "TGraphErrors.h"
22#include "TStatToolkit.h"
23
24#include "TTreeStream.h"
25#include "AliExternalTrackParam.h"
26#include "AliCDBManager.h"
27#include "AliMagF.h"
28#include "AliTrackerBase.h"
29#include "TAxis.h"
305347b9 30#include "TLegend.h"
31#include "AliGeomManager.h"
e1ed5ec5 32#include "TCut.h"
33#include "TCanvas.h"
34#include "TH2.h"
35#include "TF1.h"
36#include "TStyle.h"
37
c6b97cea 38
39
40void testdEdxGEM(const Int_t ntracks=10000, Double_t clNoise=2, Double_t corrNoiseAdd=0.15, Double_t sigmaLandau=0.26);
82194cd3 41void testExtrapolationError(Int_t ntracks, Int_t useGeo=kFALSE, Int_t seed=0);
c6b97cea 42
43
82194cd3 44void fastToyMCMI(Int_t action=1, Float_t arg0=200, Float_t arg1=0, Float_t arg2=0){
45 //
46 //
47 //
48 if (action==1) testExtrapolationError(arg0,arg1,arg2);
49}
c6b97cea 50
51
52void testdEdxGEM(const Int_t ntracks, Double_t clNoise, Double_t corrNoiseAdd, Double_t sigmaLandau){
e1ed5ec5 53 //
54 // test dEdx performance as function of the electron transparency.
55 //
56 // For simulation standard gRandom->Landau() generator with mean 15 is used
57 // Charge between pad-rows are corelated via diffusion - qm = (0.15*q-1+0.7*q0+0.15*q1)
58 // Electrons are randomly absorbed dependending on the transparency parameter
59 //
60 // Parameters:
61 // clNoise - noise integrated by cluster (snoise~1 x 4 bins ~ 2)
62 // corrNoiseAdd - additional correlated noise to be added
63 // sigmaLandau - relative sigma of Landau distirbution
64 // Setup emulation
65 // pp setup testdEdxGEM(20000,2,0.15,0.26)
66 // PbPb setup testdEdxGEM(20000,4,0.5,0.26)
67 // PbPb setup testdEdxGEM(20000,4,0.6,0.26)
68 // PbPb setup testdEdxGEM(20000,4,0.7,0.26)
69 //
70
71 TTreeSRedirector *pcstream = new TTreeSRedirector("testdEdxResolution.root","update");
72 TH2F * phisdEdx = (TH2F*)pcstream->GetFile()->Get("hisdEdx");
73 if (!phisdEdx){
d9502b8d 74
e1ed5ec5 75 TVectorD qVector(160);
76 TVectorD qVectorCorr(160);
77 TVectorD qVectorAcc(160);
78 TVectorD qVectorSorted(160);
79 Int_t indexes[160];
80 Int_t ntrans=20;
81 TVectorD vecdEdx(ntrans);
82 TVectorD vecT(ntrans);
83 //
84
85 for (Int_t itrack=0; itrack<ntracks; itrack++){
86 Double_t meanEl=15*(1.+1*gRandom->Rndm());
87 Double_t sigmaEl=sigmaLandau*meanEl*TMath::Power(15./meanEl,0.25);
88 for (Int_t irow=0; irow<160; irow++){
89 qVector[irow]=gRandom->Landau(meanEl, sigmaEl);
90 }
91 qVectorCorr[0]=qVector[0];
92 qVectorCorr[158]=qVector[158];
93 for (Int_t irow=1; irow<159; irow++){ //corralte measurement via diffusion
94 qVectorCorr[irow]= 0.15*(qVector[irow-1]+ qVector[irow+1])+0.7*qVector[irow];
95 }
96
97 for (Int_t itrans=0; itrans<ntrans; itrans++){
98 Double_t transparency=(itrans+1.)/ntrans;
99 vecT[itrans]=transparency;
100 for (Int_t irow=0; irow<160; irow++) {
101 qVectorAcc[irow]=0;
102 for (Int_t iel=0; iel<TMath::Nint(qVectorCorr[irow]); iel++) {
103 if (gRandom->Rndm()<transparency) qVectorAcc[irow]+=1./transparency;
104 }
105 }
106 for (Int_t irow=0; irow<160; irow++) {
107 qVectorAcc[irow]+=gRandom->Gaus(0, clNoise);
108 if (qVectorAcc[irow]<0) qVectorAcc[irow]=0;
109 }
110 TMath::Sort(160,qVectorAcc.GetMatrixArray(), indexes, kFALSE);
111 for (Int_t irow=0; irow<160; irow++) {
112 qVectorSorted[irow]=qVectorAcc[indexes[irow]];
113 }
114 vecdEdx[itrans]=TMath::Mean(0.6*160., qVectorSorted.GetMatrixArray());
115 vecdEdx[itrans]+=gRandom->Gaus(0, corrNoiseAdd);
116 }
117 (*pcstream)<<"dEdx"<<
118 "itrack="<<itrack<< // itrack
119 "meanEl="<<meanEl<< //
120 "sigmaEl="<<sigmaEl<< //
121 "vecT.="<<&vecT<< //
122 "vecdEdx.="<<&vecdEdx<<
123 "qVector.="<<&qVector<<
124 "\n";
125 }
126 TTree * tree = (TTree*)(pcstream->GetFile()->Get("dEdx"));
127 tree->Draw("vecdEdx.fElements/(meanEl):vecT.fElements>>hisdEdx(16,0.199,1,100,0.70,1.50)","meanEl<100","colz");
128 phisdEdx= (TH2F*)tree->GetHistogram()->Clone();
129 gStyle->SetOptFit(1);
130 gStyle->SetOptStat(0);
131 phisdEdx->Write("hisdEdx");
132 }
133 delete pcstream;
134 gStyle->SetOptStat(0);
135 gStyle->SetOptFit(1);
136 pcstream = new TTreeSRedirector("testdEdxResolution.root","update");
137 phisdEdx = (TH2F*)pcstream->GetFile()->Get("hisdEdx");
138 TObjArray *arrFit = new TObjArray(3);
139 phisdEdx->FitSlicesY(0,0,-1,0,"QNR",arrFit);
140 TH1D * hisRes = (TH1D*)arrFit->At(2);
141 hisRes->Divide((TH1D*)arrFit->At(1));
142 hisRes->Scale(100);
143 hisRes->SetTitle("dEdx resolution");
144 hisRes->SetMarkerStyle(21);
145 hisRes->GetXaxis()->SetTitle("Eff. electron transparency");
146 hisRes->GetYaxis()->SetTitle("#sigma_{dEdx}/dEdx (%)");
147 hisRes->SetMinimum(4);
148 hisRes->SetMaximum(8);
149 //
150 TF1 * ftrans = new TF1("ftrans","[0]*x**(-(abs([1])+0.000001))");
151 ftrans->SetParName(0,"#sigma_{T1}");
152 ftrans->SetParName(1,"#sigma slope");
153 ftrans->SetParameters(0.05,1,0.05);
154 hisRes->SetMarkerStyle(21);
155 hisRes->SetMarkerSize(0.9);
156 TCanvas * canvasTrans = new TCanvas("canvasTrans", "canvasTrans",600,500);
157 canvasTrans->SetTicks(1,1);
158 hisRes->Fit(ftrans,"","",0.2,1);
159 canvasTrans->SaveAs("canvasElTrans.pdf");
160 //TH
161
162
d9502b8d 163}
164
82194cd3 165void testExtrapolationError(Int_t ntracks, Int_t useGeo, Int_t seed){
d9502b8d 166 //
167 // check the extrapolation error
168 // 2 scenarios
169 // - ITS extrapolation
170 // - ITS + TRD interpolation
171 //
172 //
82194cd3 173 gRandom->SetSeed(seed);
d9502b8d 174 const char * ocdbpath = "local://$ALICE_ROOT/OCDB/";
175 AliCDBManager * man = AliCDBManager::Instance();
176 man->SetDefaultStorage(ocdbpath);
177 man->SetRun(0);
178 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG, AliMagF::kBeamTypepp, 2.76/2.));
305347b9 179 Double_t bz=AliTrackerBase::GetBz();
82194cd3 180 if ( useGeo) AliGeomManager::LoadGeometry("geometry.root");
181 if (!useGeo) AliGeomManager::LoadGeometry();
305347b9 182 //
82194cd3 183 // Double_t rpoints[13]={2.2, 2.8, 3.6, 20, 22, 41, 43, 300,315,330,345,360,375};
d9502b8d 184 Double_t spoints[13]={0.0004,0.0004,0.0004,0.0004,0.0004,0.0004, 0.0004, 0.02,0.02,0.02,0.02,0.02,0.02}; // ITS layers R poition (http://arxiv.org/pdf/1304.1306v3.pdf - pixel scenario)
82194cd3 185 //
186 Double_t rpoints[13]={2.32, 3.13, 3.91, 19.41, 24.71, 35.33, 40.53, 300,315,330,345,360,375};
187
d9502b8d 188 //
189 //
190 //
191 const Int_t nbins=40;
192 TFile * f = TFile::Open("testExtrapolationErr.root","update");
193 if (f->Get("extrapol")==0){
194 delete f;
195 f=0;
305347b9 196 Double_t cv[21]={0};
197 Double_t covarSeed[15]={0};
198 for (Int_t i=0; i<15; i++) covarSeed[i]=0;
199 covarSeed[0]=0.1;
200 covarSeed[2]=0.1;
201 covarSeed[5]=0.001;
202 covarSeed[9]=0.001;
c6b97cea 203 covarSeed[14]=0.03;
305347b9 204 //
205 Double_t vertex[3]={0,0,0};
d9502b8d 206 TTreeSRedirector *pcstream = new TTreeSRedirector("testExtrapolationErr.root","update");
207 //
208 TVectorD vecR(nbins);
209 TVectorD vecITSErr0(nbins);
210 TVectorD vecITSTRDErr0(nbins);
c6b97cea 211 //
d9502b8d 212 for (Int_t itrack=0; itrack<ntracks; itrack++){
213 //
214 //Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
82194cd3 215 if (itrack%20==0) printf("processing\t%d\n", itrack);
216 Double_t phi = (gRandom->Rndm()-1)*TMath::Pi()/18;
d9502b8d 217 Double_t eta = gRandom->Uniform(-1, 1);
82194cd3 218 Double_t pt = 1./(gRandom->Rndm()+0.00001);
d9502b8d 219 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
220 Double_t pxyz[3];
221 pxyz[0]=pt*TMath::Cos(phi);
222 pxyz[1]=pt*TMath::Sin(phi);
223 pxyz[2]=pt*TMath::Tan(theta);
224 Short_t psign=(gRandom->Rndm()>0.5)? -1.:1.;
d9502b8d 225 AliExternalTrackParam *track= new AliExternalTrackParam(vertex, pxyz, cv, psign);
82194cd3 226 Double_t alpha=TMath::ATan2(pxyz[1],pxyz[0]);
227 track->Rotate(alpha);
305347b9 228 //
c6b97cea 229 // 0.) Estimate the error using the ITS extrapolation and ITS+TRD interpolation - neglecting the MS -inf. momanta tracks
d9502b8d 230 //
231 for (Int_t irbin=0; irbin<nbins; irbin++){
232 Double_t x0 =85+Double_t(irbin)*(245.-85.)/Double_t(nbins);
233 Double_t x;
234 //
235 // parabolic fit
236 //
237 TLinearFitter fitterITS(3,"pol2");
238 TLinearFitter fitterITSTRD(3,"pol2");
239 vecR[irbin]=x0;
240 for (Int_t i=0; i<13; i++) {
241 Double_t y = 0;
242 if (track->GetYAt(rpoints[irbin],bz,y)){
243 x=rpoints[i]-x0;
244 if (i<7) fitterITS.AddPoint(&x,y,spoints[i]);
245 fitterITSTRD.AddPoint(&x,y,spoints[i]);
246 }
247 }
248 fitterITS.Eval();
249 fitterITSTRD.Eval();
250 vecITSErr0[irbin]=fitterITS.GetParError(0);
251 vecITSTRDErr0[irbin]=fitterITSTRD.GetParError(0);
d9502b8d 252 }
305347b9 253 //
c6b97cea 254 // 1.) estimate q/pt resolution for the ITS+TPC, ITS+TPC+TRD and ITS+TRD scenario
305347b9 255 //
256 AliExternalTrackParam * param = new AliExternalTrackParam(*track);
257 Double_t *covar = (Double_t*) param->GetCovariance();
258 for (Int_t i=0; i<15; i++) covar[i]=covarSeed[i];
259 AliTrackerBase::PropagateTrackToBxByBz(param, 370,0.13,1,kFALSE);
260
261 AliExternalTrackParam *trackITSTPC= new AliExternalTrackParam(*param);
262 AliExternalTrackParam *trackITSTPCTRD= new AliExternalTrackParam(*param);
263 AliExternalTrackParam *trackITSTRD= new AliExternalTrackParam(*param);
c6b97cea 264 AliExternalTrackParam *trackTPCTRD= new AliExternalTrackParam(*param);
305347b9 265 //
266 Bool_t tStatus=kTRUE;
267 for (Int_t idet=2;idet>=0; idet--){
268 Int_t nlayers=7;
269 if (idet==1) nlayers=159;
270 if (idet==2) nlayers=6;
271 for (Int_t ilayer=nlayers; ilayer>=0; ilayer--){
272 Double_t rlayer=245.-ilayer;
273 Double_t slayer=0.1;
274 if (idet==0) {rlayer=rpoints[ilayer]; slayer=spoints[ilayer];}
275 if (idet==2) {rlayer=rpoints[ilayer+7]; slayer=spoints[ilayer+7];}
276 param->PropagateTo(rlayer,bz);
277 tStatus&=!AliTrackerBase::PropagateTrackToBxByBz(param, rlayer,0.13,1,kFALSE);
278 tStatus&=!AliTrackerBase::PropagateTrackToBxByBz(trackITSTPC, rlayer,0.13,1,kFALSE);
279 tStatus&=!AliTrackerBase::PropagateTrackToBxByBz(trackITSTPCTRD, rlayer,0.13,1,kFALSE);
280 tStatus&=!AliTrackerBase::PropagateTrackToBxByBz(trackITSTRD, rlayer,0.13,1,kFALSE);
281 tStatus&=!AliTrackerBase::PropagateTrackToBxByBz(trackTPCTRD, rlayer,0.13,1,kFALSE);
282 //if (tStatus==kFALSE) break;
283 Double_t pointPos[2]={param->GetY(),param->GetZ()};
284 Double_t pointCov[3]= { slayer*slayer,0,slayer*slayer};
285 tStatus&=!trackITSTPCTRD->Update(pointPos,pointCov);
286 if (idet!=2) {
287 trackITSTPC->Update(pointPos,pointCov);
288 }
289 if (idet!=1) {
290 trackITSTRD->Update(pointPos,pointCov);
291 }
292 if (idet!=0) {
293 trackTPCTRD->Update(pointPos,pointCov);
294 }
295 }
296 }
c6b97cea 297 //
298 // 2.) Estimate propagation error at given radius
299 //
300 // 2.a) Fit ITS and TRD track (gauss smeared point error 0 not MS used)
301 AliExternalTrackParam *trackITS= new AliExternalTrackParam(*track);
302 AliExternalTrackParam *trackITS0= new AliExternalTrackParam(*track);
303 covar = (Double_t*) trackITS->GetCovariance();
304 for (Int_t i=0; i<15; i++) covar[i]=covarSeed[i];
305 AliExternalTrackParam *trackTRD= new AliExternalTrackParam(*param);
306 AliExternalTrackParam *trackTRD0= new AliExternalTrackParam(*param);
307 {for (Int_t ilayer=0; ilayer<7; ilayer++){
308 AliTrackerBase::PropagateTrackToBxByBz(trackITS, rpoints[ilayer],0.13,1,kFALSE);
309 AliTrackerBase::PropagateTrackToBxByBz(trackITS0, rpoints[ilayer],0.13,1,kFALSE);
310 Double_t pointPos[2]={trackITS0->GetY()+gRandom->Gaus(0,spoints[ilayer]),trackITS0->GetZ()+gRandom->Gaus(0,spoints[ilayer])};
311 Double_t pointCov[3]= {spoints[ilayer]*spoints[ilayer],0,spoints[ilayer]*spoints[ilayer]};
312 trackITS->Update(pointPos,pointCov);
313 }}
314 {for (Int_t ilayer=5; ilayer>=0; ilayer--){
315 AliTrackerBase::PropagateTrackToBxByBz(trackTRD, rpoints[ilayer+7],0.13,1,kFALSE);
316 AliTrackerBase::PropagateTrackToBxByBz(trackTRD0, rpoints[ilayer+7],0.13,1,kFALSE);
317 Double_t pointPos[2]={trackTRD0->GetY()+gRandom->Gaus(0,spoints[ilayer+7]),trackTRD0->GetZ()+gRandom->Gaus(0,spoints[ilayer+7])};
318 Double_t pointCov[3]= {spoints[ilayer+7]*spoints[ilayer+7],0,spoints[ilayer+7]*spoints[ilayer+7]};
319 trackTRD->Update(pointPos,pointCov);
320 }}
321 //
322 // 2.b) get ITS extrapolation and TRD interpolation errors at random radisues positions
323 //
324 const Int_t npoints=20;
325 TVectorD vecRKalman(npoints);
326 TVectorD vecITSDeltaY(npoints), vecITSTRDDeltaY(npoints);
327 TVectorD vecITSErrY(npoints), vecITSTRDErrY(npoints);
328 for (Int_t ipoint=0; ipoint<npoints; ipoint++){
329 Double_t rLayer= 85.+gRandom->Rndm()*(245.-85.);
330 AliExternalTrackParam *trackITSP= new AliExternalTrackParam(*trackITS);
331 AliExternalTrackParam *trackITSP0= new AliExternalTrackParam(*trackITS0);
332 AliExternalTrackParam *trackTRDP= new AliExternalTrackParam(*trackTRD);
333 AliTrackerBase::PropagateTrackToBxByBz(trackITSP, rLayer,0.13,1,kFALSE);
334 AliTrackerBase::PropagateTrackToBxByBz(trackITSP0, rLayer,0.13,1,kFALSE);
335 AliTrackerBase::PropagateTrackToBxByBz(trackTRDP, rLayer,0.13,1,kFALSE);
336 vecRKalman[ipoint]=rLayer;
337 trackTRDP->Rotate(trackITSP->GetAlpha());
338 vecITSDeltaY[ipoint]=trackITSP->GetY()-trackITSP0->GetY();
339 vecITSErrY[ipoint]=TMath::Sqrt(trackITSP->GetSigmaY2());
340 AliTrackerBase::UpdateTrack(*trackITSP, *trackTRDP);
341 vecITSTRDDeltaY[ipoint]=trackITSP->GetY()-trackITSP0->GetY();
342 vecITSTRDErrY[ipoint]=TMath::Sqrt(trackITSP->GetSigmaY2());
343 }
344 //
305347b9 345
d9502b8d 346 //vecITSErr0.Print();
347 (*pcstream)<<"extrapol"<<
348 "itrack="<<itrack<<
305347b9 349 "track.="<<track<<
d9502b8d 350 "vecR.="<<&vecR<<
c6b97cea 351 //
d9502b8d 352 "vecITSErr0.="<<&vecITSErr0<<
353 "vecITSTRDErr0.="<<&vecITSTRDErr0<<
305347b9 354 "tStatus="<<tStatus<<
355 "trackITSTPC.="<<trackITSTPC<<
356 "trackITSTPCTRD.="<<trackITSTPCTRD<<
357 "trackITSTRD.="<<trackITSTRD<<
358 "trackTPCTRD.="<<trackTPCTRD<<
c6b97cea 359 //
360 // kalman extapoltation einterplotaltion studies
361 // extrapolation and ITSTRDitnerpolation at different radiuses
362 "verRKalman.="<<&vecRKalman<<
363 "vecITSDeltaY.="<<&vecITSDeltaY<<
364 "vecITSErrY.="<<&vecITSErrY<<
365 "vecITSTRDDeltaY.="<<&vecITSTRDDeltaY<<
366 "vecITSTRDErrY.="<<&vecITSTRDErrY<<
d9502b8d 367 "\n";
368 }
369 delete pcstream;
370 }
371 delete f;
2dd02504 372
373 gStyle->SetOptTitle(0);
d9502b8d 374 f = TFile::Open("testExtrapolationErr.root","update");
375 TTree * tree = (TTree*)f->Get("extrapol");
82194cd3 376 tree->SetMarkerStyle(25);
377 tree->SetMarkerSize(0.3);
378
d9502b8d 379 //
380 TGraphErrors * grITS0 = TStatToolkit::MakeGraphErrors(tree,"10*vecITSErr0.fElements:vecR.fElements","",25,4,0);
381 TGraphErrors * grITSTRD0 = TStatToolkit::MakeGraphErrors(tree,"10*vecITSTRDErr0.fElements:vecR.fElements","",21,2,0);
382 //
383 grITS0->GetXaxis()->SetTitle("radius (cm)");
2dd02504 384 grITS0->GetYaxis()->SetTitle("#sigma_{r#varphi} (mm)");
d9502b8d 385 grITS0->Draw("ap");
386 grITSTRD0->Draw("p");
2dd02504 387 TLegend * legend = new TLegend(0.11,0.65,0.55,0.89,"Track residuals at TPC (q/p_{t}=0)");
d9502b8d 388 legend->AddEntry(grITS0,"ITS extrapolation","p");
305347b9 389 legend->AddEntry(grITSTRD0,"ITS-TRD interpolation","p");
2dd02504 390 legend->SetFillColor(10);
d9502b8d 391 legend->Draw();
305347b9 392 //
393 //
394 //
395 TCut cut="abs(track.fP[4])<0.25";
396 TGraphErrors * grITSTPCqPt = TStatToolkit::MakeGraphErrors(tree,"sqrt(trackITSTPC.fC[14]):abs(track.fP[4])",cut,20,1,0.8);
397 TGraphErrors * grITSTRDqPt = TStatToolkit::MakeGraphErrors(tree,"sqrt(trackITSTRD.fC[14]):abs(track.fP[4])",cut,21,4,0.8);
398 TGraphErrors * grITSTPCTRDqPt = TStatToolkit::MakeGraphErrors(tree,"sqrt(trackITSTPCTRD.fC[14]):abs(track.fP[4])",cut,24,2,0.8);
399 TGraphErrors * grTPCTRDqPt = TStatToolkit::MakeGraphErrors(tree,"sqrt(trackTPCTRD.fC[14]):abs(track.fP[4])",cut,25,3,0.8);
400 grITSTPCqPt->GetXaxis()->SetTitle("q/p_{t} (c/GeV)");
401 grITSTPCqPt->GetYaxis()->SetTitle("#sigma_{q/p_{t}} (c/GeV)");
402 grITSTPCqPt->SetMaximum(0.003);
403 TCanvas * canvasResol = new TCanvas("canvasResol","canvasResol",800,600);
404 canvasResol->cd();
405 canvasResol->SetLeftMargin(0.15);
406 grITSTPCqPt->GetYaxis()->SetTitleOffset(1.3);
d9502b8d 407
305347b9 408 {
409 grITSTPCqPt->Draw("ap");
410 grITSTRDqPt->Draw("p");
411 grITSTPCTRDqPt->Draw("p");
412 grTPCTRDqPt->Draw("p");
413 TLegend * legendQpt = new TLegend(0.41,0.11,0.89,0.4,"q/p_{t} resolution (from covariance matrix) ");
414 legendQpt->AddEntry(grITSTPCqPt,"ITS+TPC","p");
415 legendQpt->AddEntry(grITSTRDqPt,"ITS+TRD","p");
416 legendQpt->AddEntry(grITSTPCTRDqPt,"ITS+TPC+TRD","p");
417 legendQpt->AddEntry(grTPCTRDqPt,"TPC+TRD","p");
418 legendQpt->Draw();
419 }
c6b97cea 420 //
82194cd3 421 //
422 TCanvas *canvasResolution = new TCanvas("canvasResolution","canvasResolution",600,600);
c6b97cea 423 canvasResolution->Divide(1,3);
424 //
425 canvasResolution->cd(1);
426 tree->Draw("vecITSErrY.fElements:verRKalman.fElements:abs(track.fP[4])","","colz");
427 canvasResolution->cd(2);
428 tree->Draw("vecITSTRDErrY.fElements:verRKalman.fElements:abs(track.fP[4])","","colz");
429 canvasResolution->cd(3);
430 tree->Draw("vecITSTRDErrY.fElements/vecITSErrY.fElements:verRKalman.fElements:abs(track.fP[4])","","colz");
431 canvasResolution->SaveAs("canvasResolution.pdf");
432 canvasResolution->SaveAs("canvasResolution.png");
82194cd3 433 //
434 //
435 //
436 TStatToolkit toolkit;
437 Double_t chi2=0;
438 Int_t npoints=0;
439 TVectorD param;
440 TMatrixD covar;
441 //
442 TCut cut="";
443 TString fstringITS="";
444 fstringITS+="abs(track.fP[4])*(verRKalman.fElements-40)^2++";
445 fstringITS+="(verRKalman.fElements-40)^2++";
446 TString * fitErrITS = TStatToolkit::FitPlane(tree,"vecITSErrY.fElements:(0.1+track.fP[4]**2)", fstringITS.Data(),"verRKalman.fElements>0", chi2,npoints,param,covar,-1,0,180000 , kFALSE);
447 tree->SetAlias("fitErrITS",fitErrITS->Data());
448 fitErrITS->Tokenize("++")->Print();
449 tree->Draw("vecITSErrY.fElements/fitErrITS:verRKalman.fElements:abs(track.fP[4])","","colz");
450 //
451 //
452 TString fstringITSTRD="";
453 fstringITSTRD+="abs(track.fP[4])++";
454 fstringITSTRD+="abs(track.fP[4]*verRKalman.fElements)++";
455 fstringITSTRD+="abs(track.fP[4]*verRKalman.fElements**2)++";
456 for (Int_t iter=0; iter<3; iter++){
457 TCut cutOut="1";
458 if (iter>0) cutOut="abs(vecITSTRDErrY.fElements/fitErrITSTRD-1)<0.8";
459 TString * fitErrITSTRD = TStatToolkit::FitPlane(tree,"vecITSTRDErrY.fElements:(0.1+track.fP[4]**2)", fstringITSTRD.Data(),"verRKalman.fElements>0"+cutOut, chi2,npoints,param,covar,-1,0,180000 , kFALSE);
460 tree->SetAlias("fitErrITSTRD",fitErrITSTRD->Data());
461 fitErrITSTRD->Tokenize("++")->Print();
462 }
463 tree->Draw("vecITSTRDErrY.fElements/fitErrITSTRD:verRKalman.fElements:abs(track.fP[4])","","colz");
464
465
466 TCanvas *canvasResolutionFit = new TCanvas("canvasResolutionFit","canvasResolutionFit",700,700);
467 canvasResolutionFit->Divide(1,3);
468 canvasResolutionFit->cd(1);
469 tree->Draw("fitErrITS:verRKalman.fElements:abs(track.fP[4])>>hITS","","colz");
470 htemp = (TH2F*)gPad->GetPrimitive("hITS")->Clone();
471 htemp->GetXaxis()->SetTitle("#it{r} (cm)");
472 htemp->GetYaxis()->SetTitle("#sigma_{r#phi} (cm) ITS");
473 htemp->GetZaxis()->SetTitle("q/#it{p_{t}} (#it{c}/GeV)");
474 htemp->Draw("colz");
475 canvasResolutionFit->cd(2);
476 tree->Draw("fitErrITSTRD:verRKalman.fElements:abs(track.fP[4])>>hITSTRD","","colz");
477 htemp = (TH2F*)gPad->GetPrimitive("hITSTRD")->Clone();
478 htemp->GetXaxis()->SetTitle("#it{r} (cm)");
479 htemp->GetYaxis()->SetTitle("#sigma_{r#phi} (cm) ITS+TRD");
480 htemp->GetZaxis()->SetTitle("q/#it{p_{t}} (#it{c}/GeV)");
481 htemp->Draw("colz");
482 canvasResolutionFit->cd(3);
483 tree->Draw("fitErrITSTRD:verRKalman.fElements:abs(track.fP[4])>>hITSTRD","","colz");
484 htemp = (TH2F*)gPad->GetPrimitive("hITSTRD")->Clone();
485 htemp->GetXaxis()->SetTitle("#it{r} (cm)");
486 htemp->GetYaxis()->SetTitle("#sigma_{r#phi} (cm) ITS+TRD");
487 htemp->GetZaxis()->SetTitle("q/#it{p_{t}} (#it{c}/GeV)");
488 htemp->Draw("colz");
489 canvasResolutionFit->SaveAs("canvasResolutionFit.pdf");
490 canvasResolutionFit->SaveAs("canvasResolutionFit.png");
491
d9502b8d 492}
464fc663 493
494
495void DrawInerpolationResolution(){
496
497 // resRphi = 0.004390 + oneOverPt*(-0.136403) + oneOverPt*radius*(0.002266) + oneOverPt*radius*radius*(-0.000006);
498
499 //TF1 f1("f1","[0]+[1]*(
500}
501
502void MakeRobustFitTest(){
503 //
504 //
505 //
506 //
507 TH2F * his2D = new TH2F("his2D", "his2D",50,0,1., 100,-2.,2.);
508 TH2F * his2DW = new TH2F("his2DW", "his2DW",50,0,1., 100,-2.,2.);
509 Double_t probRND = 0.1;
510 Int_t ntracks = 20*50000;
511 Int_t nclusters = 16;
512 Double_t sigmaCluster=0.1;
513 //
514 for (Int_t itrack=0; itrack<ntracks; itrack++){
515 Double_t x = gRandom->Rndm();
516 Double_t widthTrack = 0.1/(0.5+gRandom->Exp(0.5));
517 Double_t y = 0;
518 y= gRandom->Gaus(0.0,widthTrack);
519 Bool_t isRandom=gRandom->Rndm()<probRND;
520 //if (gRandom->Rndm()<probRND) y= -2+4*gRandom->Rndm();
521 //
522 Double_t sigmaTrack= TMath::Sqrt(sigmaCluster*sigmaCluster/nclusters+widthTrack*widthTrack);
523 for (Int_t icl=0; icl<nclusters; icl++){
524 his2D->Fill(x,y+gRandom->Gaus(0,sigmaCluster));
525 his2DW->Fill(x,y+gRandom->Gaus(0,sigmaCluster),1/sigmaTrack);
526 }
527 }
528 {
529 his2D->Draw("colz");
530 his2DW->Draw("colz");
531 his2D->FitSlicesY();
532 his2DW->FitSlicesY();
533 his2D_1->Draw();
534 his2DW_1->Draw("same");
535 }
536 TMath::RMS(50, &(his2D_1->GetArray()[1]));
537
538 TH1D * phis1D = his2D->ProjectionY("aaa",5,5);
539 Int_t nbinsY= phis1D->GetXaxis()->GetNbins();
540 TVectorD vector(nbinsY, &(phis1D->GetArray()[1]));
541
542 TStatToolkit::MakeStat1D((TH2*)his2D, 1,0.8,0,25,1)->Draw("alp");
543 TStatToolkit::MakeStat1D((TH2*)his2D, 1,0.8,2,20,2)->Draw("lp");
544 TStatToolkit::MakeStat1D((TH2*)his2D, 1,0.8,4,21,4)->Draw("lp");
545
546}