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