]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/CalibMacros/CalibPID.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
fixing warnings (Federico C.)
[u/mrichter/AliRoot.git] / TPC / CalibMacros / CalibPID.C
CommitLineData
8a92e133 1/*
2 .x ConfigOCDB.C
3
4 .x ~/rootlogon.C
5 .x ~/UliStyle.C
6 gSystem->Load("libANALYSIS");
7 gSystem->Load("libTPCcalib");
8 gSystem->Load("libSTAT.so");
9
10 .L $ALICE_ROOT/TPC/CalibMacros/CalibPID.C+
11 paramCl = AliTPCcalibDB::Instance()->GetClusterParam();
12 // init
13 Init("calibPID06");
14 SetRange(6,100,160);
15 SetRange(0,0.5,10);
16 SetRange(4,0,200);
17 SetType(1);
18 pid->MakeReport();
19
20*/
21#include "TMath.h"
22#include "TString.h"
23#include "TFile.h"
24#include "TTree.h"
25#include "TObjArray.h"
26#include "TF1.h"
27#include "TH1.h"
28#include "TH2.h"
29#include "TH3F.h"
30#include "THnSparse.h"
31#include "TProfile.h"
32#include "TCut.h"
33#include "TMatrixD.h"
34#include "TVectorD.h"
35//
36#include "AliCDBManager.h"
37#include "AliCDBMetaData.h"
38#include "AliCDBId.h"
39#include "AliCDBRunRange.h"
40#include "AliCDBStorage.h"
41#include "AliTPCClusterParam.h"
42//
43#include "AliTPCcalibPID.h"
44#include "AliTPCcalibDB.h"
45#include "TStatToolkit.h"
46
47AliTPCClusterParam * paramCl=0;
48AliTPCcalibPID * pid =0;
49TObjArray fitArr;
50TTree * treeDump =0;
51TTree * treeQtot;
52TTree * treeQmax;
53TTree * treeRatioQmax;
54TTree * treeRatioQtot;
55Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
56Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
57
58// 0 1 2 3 4 5 6 7
59// dE/dx, z, phi, theta, p, bg, ncls type
60//Int_t binsQA[7] = {150, 10, 10, 10, 50, 50, 8};
61
62void Init(char* name="calibPID06"){
63 //
64 //
65 TFile fcalib("CalibObjects.root");
66 TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
67 pid = ( AliTPCcalibPID *)array->FindObject(name);
68 TString axisName[9];
69 axisName[0] ="dE/dx"; axisName[1] ="z (cm)";
70 axisName[2] ="sin(#phi)"; axisName[3] ="tan(#theta)";
71 axisName[4] ="p (GeV)"; axisName[5] ="#beta#gamma";
72 axisName[6] ="N_{cl}";
73
74 pid->GetHistQtot()->SetTitle("Q_{tot};(z,sin(#phi),tan(#theta),p,betaGamma,ncls); TPC signal Q_{tot} (a.u.)");
75 pid->GetHistQmax()->SetTitle("Q_{max};(z,sin(#phi),tan(#theta),p,betaGamma,ncls); TPC signal Q_{max} (a.u.)");
76
77 for (Int_t ivar2=0;ivar2<7;ivar2++){
78 pid->GetHistQmax()->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
79 pid->GetHistQtot()->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
80 pid->GetHistRatioMaxTot()->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
81 pid->GetHistRatioQmax()->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
82 pid->GetHistRatioQtot()->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
83 //
84 pid->GetHistQmax()->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
85 pid->GetHistQtot()->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
86
87 pid->GetHistRatioMaxTot()->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
88 pid->GetHistRatioQmax()->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
89 pid->GetHistRatioQtot()->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
90 }
91}
92
93
94void StoreParam(){
95 Int_t runNumber = 0;
96 AliCDBMetaData *metaData= new AliCDBMetaData();
97 metaData->SetObjectClassName("AliTPCClusterParam");
98 metaData->SetResponsible("Marian Ivanov");
99 metaData->SetBeamPeriod(1);
100 metaData->SetAliRootVersion("05-23-02");
101 metaData->SetComment("October runs calibration");
102 AliCDBId id1("TPC/Calib/ClusterParam", runNumber, AliCDBRunRange::Infinity());
103 AliCDBStorage *gStorage = AliCDBManager::Instance()->GetStorage("local:///u/miranov/OCDB/TPCmc");
104 gStorage->Put(paramCl, id1, metaData);
105
106}
107
108
109
110
111void SetRange(Int_t index, Float_t min, Float_t max){
112 pid->GetHistQmax()->GetAxis(index)->SetRangeUser(min,max);
113 pid->GetHistQtot()->GetAxis(index)->SetRangeUser(min,max);
114
115 pid->GetHistRatioMaxTot()->GetAxis(index)->SetRangeUser(min,max);
116 pid->GetHistRatioQtot()->GetAxis(index)->SetRangeUser(min,max);
117 pid->GetHistRatioQmax()->GetAxis(index)->SetRangeUser(min,max);
118 pid->GetHistRatioTruncQtot()->GetAxis(index)->SetRangeUser(min,max);
119 pid->GetHistRatioTruncQmax()->GetAxis(index)->SetRangeUser(min,max);
120}
121
122
123void SetType(Int_t type){
124 pid->GetHistQmax()->GetAxis(7)->SetRange(type,type);
125 pid->GetHistQtot()->GetAxis(7)->SetRange(type,type);
126 //
127 //
128 pid->GetHistRatioMaxTot()->GetAxis(7)->SetRange(type,type);
129 pid->GetHistRatioQtot()->GetAxis(7)->SetRange(type,type);
130 pid->GetHistRatioQmax()->GetAxis(7)->SetRange(type,type);
131 pid->GetHistRatioTruncQtot()->GetAxis(7)->SetRange(type,type);
132 pid->GetHistRatioTruncQmax()->GetAxis(7)->SetRange(type,type);
133
134}
135
136
137
138
139void ReadTrees(){
140 TFile f0("dumpQtot.root");
141 treeQtot = (TTree*)f0.Get("Dump");
142 TFile f1("dumpQmax.root");
143 treeQmax = (TTree*)f1.Get("Dump");
144 TFile f2("dumpRatioQtot.root");
145 treeRatioQtot = (TTree*)f2.Get("Dump");
146 TFile f3("dumpRatioQmax.root");
147 treeRatioQmax = (TTree*)f3.Get("Dump");
148}
149
150
151
152void Fit(){
153 //
154 // align pads
155 //
156 TStatToolkit toolkit;
157 Double_t chi2;
158 TVectorD paramTot[5], paramMax[5];
159 TString *strQT[4], *strQM[4];
160 TVectorD paramTotRatio[5], paramMaxRatio[5];
161 TString *strQTRatio[5], *strQMRatio[5];
162
163 TMatrixD covar;
164 Int_t npoints;
165 treeQmax->SetAlias("cdr","(1-1/(1+dr))");
166 treeQmax->SetAlias("cty","(1-1/sqrt(1+ty^2))");
167 treeQmax->SetAlias("ctz","(1-1/sqrt(1+p3^2))");
168 treeQtot->SetAlias("cdr","(1-1/(1+dr))");
169 treeQtot->SetAlias("cty","(1-1/sqrt(1+ty^2))");
170 treeQtot->SetAlias("ctz","(1-1/sqrt(1+p3^2))");
171 //
172 treeRatioQmax->SetAlias("cdr","(1-1/(1+dr))");
173 treeRatioQmax->SetAlias("cty","(1-1/sqrt(1+ty^2))");
174 treeRatioQmax->SetAlias("ctz","(1-1/sqrt(1+p3^2))");
175 treeRatioQtot->SetAlias("cdr","(1-1/(1+dr))");
176 treeRatioQtot->SetAlias("cty","(1-1/sqrt(1+ty^2))");
177 treeRatioQtot->SetAlias("ctz","(1-1/sqrt(1+p3^2))");
178
179
180 TString strFit="";
181 //
182 strFit+="dr++";
183 strFit+="abs(p3)++";
184 //
185 strFit+="(1-(sy*sz)/0.5)++";
186 strFit+="(1-(sy/0.5))++";
187 strFit+="(1-sz)++";
188
189 //
190 TCut cutAll="(p>1&&dr>0.1&&val<3)";
191 TCut cutPads[4];
192 TCut cutPadsW[4];
193 for (Int_t ipad=0;ipad<4;ipad++){
194 printf("fitting pad\t%d\n",ipad);
195 cutPads[ipad]=Form("abs(type-%f)<0.2&&(p>1&&dr>0.1&&val<3)",Double_t(ipad+1.5));
196 cutPadsW[ipad]=Form("(abs(type-%f)<0.2&&(p>1&&dr>0.1&&val<3))*bincont",Double_t(ipad+1.5));
197 //
198 strQT[ipad] = toolkit.FitPlane(treeQtot,"val:1/sqrt(bincont)",strFit.Data(),cutAll+cutPads[ipad],chi2,npoints,paramTot[ipad],covar);
199 strQM[ipad] = toolkit.FitPlane(treeQmax,"val:1/sqrt(bincont)",strFit.Data(),cutAll+cutPads[ipad],chi2,npoints,paramMax[ipad],covar);
200
201 strQTRatio[ipad] = toolkit.FitPlane(treeRatioQtot,"val:1/sqrt(bincont)",strFit.Data(),cutAll+cutPads[ipad],chi2,npoints,paramTotRatio[ipad],covar);
202 strQMRatio[ipad] = toolkit.FitPlane(treeRatioQmax,"val:1/sqrt(bincont)",strFit.Data(),cutAll+cutPads[ipad],chi2,npoints,paramMaxRatio[ipad],covar);
203
204 treeRatioQmax->SetAlias(Form("fitQM%d",ipad),strQM[ipad]->Data());
205 treeRatioQmax->SetAlias(Form("fitQT%d",ipad),strQT[ipad]->Data());
206 treeRatioQtot->SetAlias(Form("fitQM%d",ipad),strQM[ipad]->Data());
207 treeRatioQtot->SetAlias(Form("fitQT%d",ipad),strQT[ipad]->Data());
208 treeQmax->SetAlias(Form("fitQM%d",ipad),strQM[ipad]->Data());
209 treeQmax->SetAlias(Form("fitQT%d",ipad),strQT[ipad]->Data());
210 treeQtot->SetAlias(Form("fitQM%d",ipad),strQM[ipad]->Data());
211 treeQtot->SetAlias(Form("fitQT%d",ipad),strQT[ipad]->Data());
212
213 treeQmax->SetAlias(Form("fitQMR%d",ipad),strQMRatio[ipad]->Data());
214 treeQmax->SetAlias(Form("fitQTR%d",ipad),strQTRatio[ipad]->Data());
215 treeQtot->SetAlias(Form("fitQMR%d",ipad),strQMRatio[ipad]->Data());
216 treeQtot->SetAlias(Form("fitQTR%d",ipad),strQTRatio[ipad]->Data());
217 treeRatioQmax->SetAlias(Form("fitQMR%d",ipad),strQMRatio[ipad]->Data());
218 treeRatioQmax->SetAlias(Form("fitQTR%d",ipad),strQTRatio[ipad]->Data());
219 treeRatioQtot->SetAlias(Form("fitQMR%d",ipad),strQMRatio[ipad]->Data());
220 treeRatioQtot->SetAlias(Form("fitQTR%d",ipad),strQTRatio[ipad]->Data());
221 }
222 //
223 //
224 //
225 for (Int_t ipad=0;ipad<4;ipad++){
226 for (Int_t icorr=1;icorr<4;icorr++) {
227 paramMax[ipad][icorr]/=paramMax[ipad][0];
228 paramTot[ipad][icorr]/=paramTot[ipad][0];
229 paramMaxRatio[ipad][icorr]/=paramMaxRatio[ipad][0];
230 paramTotRatio[ipad][icorr]/=paramTotRatio[ipad][0];
231 }
232 }
233
234 for (Int_t ipad=0;ipad<3;ipad++){
235 for (Int_t icorr=1;icorr<4;icorr++) {
236 paramMax[ipad][icorr]+=(paramMax[3][icorr]+paramMax[ipad][icorr])*0.5;
237 paramTot[ipad][icorr]+=(paramTot[3][icorr]+paramTot[ipad][icorr])*0.5;
238 }
239 }
240
241// for (Int_t ipad=0;ipad<3;ipad++){
242// TVectorD *vecMax = (TVectorD*)(paramCl->fQNormGauss->At(3*1+ipad));
243// TVectorD *vecTot = (TVectorD*)(paramCl->fQNormGauss->At(3*0+ipad));
244// for (Int_t icorr=1;icorr<4;icorr++){
245// (*vecMax)[icorr]+=paramMax[ipad][icorr]/(*vecMax)[0];
246// (*vecTot)[icorr]+=paramTot[ipad][icorr]/(*vecTot)[0];
247// }
248// }
249
250
251}
252
253void fitQdep(){
254 SetRange(6,100,160);
255 SetRange(4,0.5,2);
256 SetRange(0,0.5,1.5);
257 TF1 f1("f1","([0]+[1]/x^2)",0.5,2);
258 SetType(2);
259 (pid->GetHistRatioQtot()->Projection(0,4))->ProfileX()->Fit(&f1);
260 SetType(3);
261 (pid->GetHistRatioQtot()->Projection(0,4))->ProfileX()->Fit(&f1);
262 SetType(4);
263 (pid->GetHistRatioQtot()->Projection(0,4))->ProfileX()->Fit(&f1);
264 //
265 SetType(2);
266 (pid->GetHistRatioQmax()->Projection(0,4))->ProfileX()->Fit(&f1);
267 SetType(3);
268 (pid->GetHistRatioQmax()->Projection(0,4))->ProfileX()->Fit(&f1);
269 SetType(4);
270 (pid->GetHistRatioQmax()->Projection(0,4))->ProfileX()->Fit(&f1);
271
272
273}
274
275void LookupHisto(Int_t minTracks=200, Float_t minp=20, Float_t maxp=10000){
276 TF1 f1("myg","gaus",0,10);
277 Int_t dim[4]={0,1,2,3};
278 pid->GetHistQtot()->GetAxis(6)->SetRangeUser(100,160);
279 pid->GetHistQtot()->GetAxis(0)->SetRangeUser(0.1,2.5);
280 pid->GetHistQmax()->GetAxis(6)->SetRangeUser(100,160);
281 pid->GetHistQmax()->GetAxis(0)->SetRangeUser(0.1,2.5);
282
283 pid->GetHistQmax()->GetAxis(4)->SetRangeUser(minp,maxp);
284 pid->GetHistQtot()->GetAxis(4)->SetRangeUser(minp,maxp);
285
286
287 THnSparse *hisTot[4];
288 TH3F *hisTotMean[4];
289 Float_t meanTotPad[4];
290 Float_t rmsTotPad[4];
291 //
292 THnSparse *hisMax[4];
293 TH3F *hisMaxMean[4];
294 Float_t meanMaxPad[4];
295 Float_t rmsMaxPad[4];
296
297 TObjArray *array = new TObjArray(8);
298
299 for (Int_t ipad=0;ipad<4;ipad++){
300 pid->GetHistQtot()->GetAxis(7)->SetRange(2+ipad,2+ipad);
301 pid->GetHistQmax()->GetAxis(7)->SetRange(2+ipad,2+ipad);
302 hisTot[ipad] = pid->GetHistQtot()->Projection(4,dim);
303 hisMax[ipad] = pid->GetHistQmax()->Projection(4,dim);
304 Float_t drmin = hisTot[ipad]->GetAxis(1)->GetXmin();
305 Float_t drmax = hisTot[ipad]->GetAxis(1)->GetXmax();
306 Float_t p2min = hisTot[ipad]->GetAxis(2)->GetXmin();
307 Float_t p2max = hisTot[ipad]->GetAxis(2)->GetXmax();
308 Float_t p3min = hisTot[ipad]->GetAxis(3)->GetXmin();
309 Float_t p3max = hisTot[ipad]->GetAxis(3)->GetXmax();
310 meanTotPad[ipad] =hisTot[ipad]->Projection(0)->GetMean();
311 meanMaxPad[ipad] =hisMax[ipad]->Projection(0)->GetMean();
312 rmsTotPad[ipad] =hisTot[ipad]->Projection(0)->GetRMS();
313 rmsMaxPad[ipad] =hisMax[ipad]->Projection(0)->GetRMS();
314
315 hisTotMean[ipad]=new TH3F(Form("htot%d",ipad),Form("htot%d",ipad),
316 hisTot[ipad]->GetAxis(1)->GetNbins(), drmin,drmax,
317 hisTot[ipad]->GetAxis(2)->GetNbins(), p2min,p2max,
318 hisTot[ipad]->GetAxis(3)->GetNbins(), p3min,p3max);
319 hisMaxMean[ipad]=new TH3F(Form("hmax%d",ipad),Form("hmax%d",ipad),
320 hisMax[ipad]->GetAxis(1)->GetNbins(), drmin,drmax,
321 hisMax[ipad]->GetAxis(2)->GetNbins(), p2min,p2max,
322 hisMax[ipad]->GetAxis(3)->GetNbins(), p3min,p3max);
323 array->AddAt(hisTotMean[ipad],ipad);
324 array->AddAt(hisMaxMean[ipad],ipad+4);
325 }
326
327 for (Int_t ibindr=1;ibindr<hisTot[0]->GetAxis(1)->GetNbins()+1; ibindr++)
328 for (Int_t ibinty=1;ibinty<hisTot[0]->GetAxis(2)->GetNbins()+1; ibinty++)
329 for (Int_t ibintz=1;ibintz<hisTot[0]->GetAxis(3)->GetNbins()+1; ibintz++)
330 for (Int_t ipad=0;ipad<4; ipad++){
331 hisTotMean[ipad]->SetBinContent(ibindr,ibinty,ibintz,meanTotPad[ipad]);
332 hisMaxMean[ipad]->SetBinContent(ibindr,ibinty,ibintz,meanMaxPad[ipad]);
333 }
334
335
336 TTreeSRedirector * cstream = new TTreeSRedirector("lookupdEdx.root");
337
338 for (Int_t ibindr=1;ibindr<hisTot[0]->GetAxis(1)->GetNbins()+1; ibindr+=1)
339 for (Int_t ibinty=1;ibinty<hisTot[0]->GetAxis(2)->GetNbins()+1; ibinty+=1)
340 for (Int_t ibintz=1;ibintz<hisTot[0]->GetAxis(3)->GetNbins()+1; ibintz+=1)
341 for (Int_t ipad=0;ipad<4; ipad++){
342 //
343 Float_t dr = hisTot[ipad]->GetAxis(1)->GetBinCenter(ibindr);
344 Float_t p2 = hisTot[ipad]->GetAxis(2)->GetBinCenter(ibinty);
345 Float_t p3 = hisTot[ipad]->GetAxis(3)->GetBinCenter(ibintz);
346 Int_t sumTot=0, sumMax=0;
347 TH1D * hisproTot =0;
348 TH1D * hisproMax =0;
349 Float_t delta=0.1;
350 while(1){
351 hisTot[ipad]->GetAxis(1)->SetRangeUser(dr-delta,dr+delta);
352 hisTot[ipad]->GetAxis(2)->SetRangeUser(p2-delta,p2+delta);
353 hisTot[ipad]->GetAxis(3)->SetRangeUser(p3-delta,p3+delta);
354 //
355 hisMax[ipad]->GetAxis(1)->SetRangeUser(dr-delta,dr+delta);
356 hisMax[ipad]->GetAxis(2)->SetRangeUser(p2-delta,p2+delta);
357 hisMax[ipad]->GetAxis(3)->SetRangeUser(p3-delta,p3+delta);
358 // 4 sigma cut
359 hisTot[ipad]->GetAxis(0)->SetRangeUser(meanTotPad[ipad]-4.*rmsTotPad[ipad],meanTotPad[ipad]+4.*rmsTotPad[ipad]);
360 hisMax[ipad]->GetAxis(0)->SetRangeUser(meanMaxPad[ipad]-4.*rmsMaxPad[ipad],meanMaxPad[ipad]+4.*rmsMaxPad[ipad]);
361 hisproTot = hisTot[ipad]->Projection(0);
362 hisproMax = hisMax[ipad]->Projection(0);
363 sumTot = hisproTot->GetSum();
364 sumMax = hisproMax->GetSum();
365 if (sumTot>minTracks) break;
366 if (delta>0.35) break;
367 delete hisproTot;
368 delete hisproMax;
369 delta+=0.1;
370 }
371 //
372 Float_t meanTot = hisproTot->GetMean();
373 Float_t rmsTot = hisproTot->GetRMS();
374 Float_t meanMax = hisproMax->GetMean();
375 Float_t rmsMax = hisproMax->GetRMS();
376 //
377
378 Float_t meangTot = -1;
379 Float_t rmsgTot = -1;
380 Float_t meangMax = -1;
381 Float_t rmsgMax = -1;
382 if(sumTot>minTracks/2 && rmsTot>0.02&&rmsMax>0.02){
383 f1.SetParameters(hisproTot->GetMaximum(),meanTot,rmsTot);
384 hisproTot->Fit(&f1,"QNR","",0.1,10);
385 meangTot = f1.GetParameter(1);
386 rmsgTot = f1.GetParameter(2);
387 f1.SetParameters(hisproMax->GetMaximum(),meanMax,rmsMax);
388 hisproMax->Fit(&f1,"QNR","",0.1,10);
389 meangMax = f1.GetParameter(1);
390 rmsgMax = f1.GetParameter(2);
391 }
392
393 TH1 *htemp = 0;
394 htemp = hisTot[ipad]->Projection(1);
395 Float_t drc = htemp->GetMean();
396 delete htemp;
397 htemp = hisTot[ipad]->Projection(2);
398 Float_t p2c = htemp->GetMean();
399 delete htemp;
400 htemp = hisTot[ipad]->Projection(3);
401 Float_t p3c = htemp->GetMean();
402 delete htemp;
403 //
404 Bool_t isOK=kTRUE;
405 if (meanTot==0) {meanTot=meanTotPad[ipad]; isOK=kFALSE;}
406 if (meanMax==0) {meanMax=meanMaxPad[ipad]; isOK=kFALSE;}
407 //
408 if (TMath::Abs(rmsTot/meanTot-0.12)>0.04) {meanTot=meanTotPad[ipad]; isOK=kFALSE;}
409 if (TMath::Abs(rmsMax/meanMax-0.12)>0.04) {meanMax=meanMaxPad[ipad]; isOK=kFALSE;}
410 //
411 hisTotMean[ipad]->SetBinContent(ibindr,ibinty,ibintz,meanTot);
412 hisMaxMean[ipad]->SetBinContent(ibindr,ibinty,ibintz,meanMax);
413 //
414 printf("%d\t%f\t%f\t%f\t%f\t%f\t%f\n", ipad, dr,p2,p3, meanTot, rmsTot, meangTot);
415 printf("%d\t%f\t%f\t%f\t%f\t%f\t%f\n", ipad, dr,p2,p3, meanMax, rmsMax, meangMax);
416 delete hisproTot;
417 delete hisproMax;
418 (*cstream)<<"dumpdEdx"<<
419 "ipad="<<ipad<<
420 "isOK="<<isOK<<
421 "sumTot="<<sumTot<<
422 "sumMax="<<sumMax<<
423 "meanTotP="<<meanTotPad[ipad]<<
424 "meanMaxP="<<meanMaxPad[ipad]<<
425 "meanTot="<<meanTot<<
426 "rmsTot="<<rmsTot<<
427 "meanMax="<<meanMax<<
428 "rmsMax="<<rmsMax<<
429 //
430 "meangTot="<<meangTot<<
431 "rmsgTot="<<rmsgTot<<
432 "meangMax="<<meangMax<<
433 "rmsgMax="<<rmsgMax<<
434 //
435 "dr="<<dr<<
436 "p2="<<p2<<
437 "p3="<<p3<<
438 "drc="<<drc<<
439 "p2c="<<p2c<<
440 "p3c="<<p3c<<
441 "\n";
442 }
443 TFile *fstream = cstream->GetFile();
444 fstream->cd();
445 array->Write("histos",TObject::kSingleKey);
446 delete cstream;
447}
448
449void FitFit(Bool_t updateParam=kFALSE){
450
451 TStatToolkit toolkit;
452 Double_t chi2Tot[5],chi2Max[5];
453 TVectorD paramTot[5], paramMax[5];
454 TString *strQT[4], *strQM[4];
455 TVectorD paramTotRatio[5], paramMaxRatio[5];
456 TString *strQTRatio[5], *strQMRatio[5];
457 TMatrixD covar;
458 Int_t npoints;
459 TCut cutPads[5];
460 treeDump->SetAlias("drm","(0.5-abs(drc))");
461 treeDump->SetAlias("dri","(1-abs(drc))");
462 treeDump->SetAlias("ty","tan(asin(p2c))**2");
463 treeDump->SetAlias("tz","(abs(p3c)**2*sqrt(1+ty))");
464 //
465 TString strFit="";
466 strFit+="drm++";
467 strFit+="ty++";
468 strFit+="tz++";
469 //
470 //strFit+="sqrt(dri*ty)++";
471 //strFit+="sqrt(dri*tz)++";
472 //strFit+="sqrt(ty*tz)++";
473
474 TCut cutAll = "meangTot>0.0&&sumMax>150&&sumTot>150&&rmsgMax/rmsMax<1.5&&abs(p3)<1";
475 for (Int_t ipad=0;ipad<4;ipad++){
476 cutPads[ipad]=Form("ipad==%d",ipad);
477 //
478 strQT[ipad] = toolkit.FitPlane(treeDump,"meangTot/meanTotP",strFit.Data(),cutAll+cutPads[ipad],chi2Tot[ipad],npoints,paramTot[ipad],covar);
479 chi2Tot[ipad]=TMath::Sqrt(chi2Tot[ipad]/npoints);
480 printf("Tot%d\t%f\t%s\t\n",ipad,chi2Tot[ipad],strQT[ipad]->Data());
481 //
482 strQM[ipad] = toolkit.FitPlane(treeDump,"meangMax/meanMaxP",strFit.Data(),cutAll+cutPads[ipad],chi2Max[ipad],npoints,paramMax[ipad],covar);
483 chi2Max[ipad]=TMath::Sqrt(chi2Max[ipad]/npoints);
484 printf("Max%d\t%f\t%s\t\n",ipad,chi2Max[ipad],strQM[ipad]->Data());
485 //
486 strQTRatio[ipad] = toolkit.FitPlane(treeDump,"meangTot/meangMax",strFit.Data(),cutAll+cutPads[ipad],chi2Max[ipad],npoints,paramTotRatio[ipad],covar);
487 chi2Max[ipad]=TMath::Sqrt(chi2Max[ipad]/npoints);
488 printf("Ratio%d\t%f\t%s\t\n\n",ipad,chi2Max[ipad],strQTRatio[ipad]->Data());
489 //
490 treeDump->SetAlias(Form("fitQT%d",ipad),strQT[ipad]->Data());
491 treeDump->SetAlias(Form("fitQM%d",ipad),strQM[ipad]->Data());
492 treeDump->SetAlias(Form("fitQTM%d",ipad),strQM[ipad]->Data());
493 }
494 TMatrixD mat(6,6);
495 for (Int_t ipad=0; ipad<3; ipad++){
496 for (Int_t icorr=0; icorr<3; icorr++){
497 if (updateParam){
498 (*paramCl->fQNormCorr)(ipad+6,icorr) += paramTot[ipad][icorr+1];
499 (*paramCl->fQNormCorr)(ipad+9,icorr) += paramMax[ipad][icorr+1];
500 }
501 mat(ipad+0,icorr) = paramTot[ipad][icorr+1];
502 mat(ipad+3,icorr) = paramMax[ipad][icorr+1];
503 }
504 }
505
506
507}
508
509
510