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