4 // 1. dump information to the tree
6 gSystem->Load("libANALYSIS");
7 gSystem->Load("libTPCcalib");
8 gSystem->Load("libSTAT");
10 .L $ALICE_ROOT/TPC/CalibMacros/CalibPID.C+
12 paramCl = AliTPCcalibDB::Instance()->GetClusterParam();
14 LookupHisto() // change SetRange in LookupHisto if needed !, check with pid->GetHistQtot()->Projection(0,1)->Draw("colz")
22 gSystem->Load("libANALYSIS");
23 gSystem->Load("libTPCcalib");
24 gSystem->Load("libSTAT");
26 .L $ALICE_ROOT/TPC/CalibMacros/CalibPID.C+
28 paramCl = AliTPCcalibDB::Instance()->GetClusterParam();
30 TFile fff("lookupdEdx.root")
33 treeDump = (TTree*)fff.Get("dumpdEdx");
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,"*")
39 StoreParam("local:///lustre/alice/akalweit/OCDBforMC") // specify corresponding location before !!
47 #include "TObjArray.h"
52 #include "THnSparse.h"
58 #include "AliCDBManager.h"
59 #include "AliCDBMetaData.h"
61 #include "AliCDBRunRange.h"
62 #include "AliCDBStorage.h"
63 #include "AliTPCClusterParam.h"
65 #include "AliTPCcalibPID.h"
66 #include "AliTPCcalibDB.h"
67 #include "TStatToolkit.h"
69 AliTPCClusterParam * paramCl=0;
70 AliTPCcalibPID * pid =0;
75 TTree * treeRatioQmax;
76 TTree * treeRatioQtot;
77 Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
78 Int_t kmimarkers[10]={21,22,23,24,25,26,27,28,29,30};
81 // dE/dx, z, phi, theta, p, bg, ncls type
82 //Int_t binsQA[7] = {150, 10, 10, 10, 50, 50, 8};
84 void Init(char* name="calibPID06"){
87 TFile fcalib("CalibObjectsTrain2.root");
88 //TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib"); // old interface
89 pid = ( AliTPCcalibPID *) fcalib.Get(name);
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}";
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.)");
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());
106 pid->GetHistQmax()->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
107 pid->GetHistQtot()->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
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());
116 void StoreParam(char* localStorage = "local:///lustre/alice/akalweit/OCDB"){
118 AliCDBMetaData *metaData= new AliCDBMetaData();
119 metaData->SetObjectClassName("AliTPCClusterParam");
120 metaData->SetResponsible("Alexander Kalweit");
121 metaData->SetBeamPeriod(1);
122 metaData->SetAliRootVersion("05-24-00");
123 metaData->SetComment("October runs calibration");
124 AliCDBId id1("TPC/Calib/ClusterParam", runNumber, AliCDBRunRange::Infinity());
125 AliCDBStorage *gStorage = AliCDBManager::Instance()->GetStorage(localStorage);
126 gStorage->Put(paramCl, id1, metaData);
133 void 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);
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);
145 void SetType(Int_t type){
146 pid->GetHistQmax()->GetAxis(7)->SetRange(type,type);
147 pid->GetHistQtot()->GetAxis(7)->SetRange(type,type);
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);
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");
174 void Fit(Bool_t updateParam=kFALSE){
178 TStatToolkit toolkit;
180 TVectorD paramTot[5], paramMax[5];
181 TString *strQT[4], *strQM[4];
182 TVectorD paramTotRatio[5], paramMaxRatio[5];
183 TString *strQTRatio[5], *strQMRatio[5];
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))");
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))");
207 strFit+="(1-(sy*sz)/0.5)++";
208 strFit+="(1-(sy/0.5))++";
212 TCut cutAll="(p>1&&dr>0.1&&val<3)";
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));
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);
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);
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());
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());
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];
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;
261 (*paramCl->fQNormCorr)(ipad+6,icorr) = paramTot[ipad][icorr+1];
262 (*paramCl->fQNormCorr)(ipad+9,icorr) = paramMax[ipad][icorr+1];
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];
285 TF1 f1("f1","([0]+[1]/x^2)",0.5,2);
287 (pid->GetHistRatioQtot()->Projection(0,4))->ProfileX()->Fit(&f1);
289 (pid->GetHistRatioQtot()->Projection(0,4))->ProfileX()->Fit(&f1);
291 (pid->GetHistRatioQtot()->Projection(0,4))->ProfileX()->Fit(&f1);
294 (pid->GetHistRatioQmax()->Projection(0,4))->ProfileX()->Fit(&f1);
296 (pid->GetHistRatioQmax()->Projection(0,4))->ProfileX()->Fit(&f1);
298 (pid->GetHistRatioQmax()->Projection(0,4))->ProfileX()->Fit(&f1);
303 void 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);
307 pid->GetHistQtot()->GetAxis(0)->SetRangeUser(0.1,6.); // important adaption to be done here ...
308 pid->GetHistQmax()->GetAxis(6)->SetRangeUser(100,160);
309 pid->GetHistQmax()->GetAxis(0)->SetRangeUser(0.1,6); // important adaption to be done here ...
311 pid->GetHistQmax()->GetAxis(4)->SetRangeUser(minp,maxp);
312 pid->GetHistQtot()->GetAxis(4)->SetRangeUser(minp,maxp);
315 THnSparse *hisTot[4];
317 Float_t meanTotPad[4];
318 Float_t rmsTotPad[4];
320 THnSparse *hisMax[4];
322 Float_t meanMaxPad[4];
323 Float_t rmsMaxPad[4];
325 TObjArray *array = new TObjArray(8);
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();
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);
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]);
364 TTreeSRedirector * cstream = new TTreeSRedirector("lookupdEdx.root");
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++){
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;
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);
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);
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;
400 Float_t meanTot = hisproTot->GetMean();
401 Float_t rmsTot = hisproTot->GetRMS();
402 Float_t meanMax = hisproMax->GetMean();
403 Float_t rmsMax = hisproMax->GetRMS();
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);
422 htemp = hisTot[ipad]->Projection(1);
423 Float_t drc = htemp->GetMean();
425 htemp = hisTot[ipad]->Projection(2);
426 Float_t p2c = htemp->GetMean();
428 htemp = hisTot[ipad]->Projection(3);
429 Float_t p3c = htemp->GetMean();
433 if (meanTot==0) {meanTot=meanTotPad[ipad]; isOK=kFALSE;}
434 if (meanMax==0) {meanMax=meanMaxPad[ipad]; isOK=kFALSE;}
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;}
439 hisTotMean[ipad]->SetBinContent(ibindr,ibinty,ibintz,meanTot);
440 hisMaxMean[ipad]->SetBinContent(ibindr,ibinty,ibintz,meanMax);
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);
446 (*cstream)<<"dumpdEdx"<<
451 "meanTotP="<<meanTotPad[ipad]<<
452 "meanMaxP="<<meanMaxPad[ipad]<<
453 "meanTot="<<meanTot<<
455 "meanMax="<<meanMax<<
458 "meangTot="<<meangTot<<
459 "rmsgTot="<<rmsgTot<<
460 "meangMax="<<meangMax<<
461 "rmsgMax="<<rmsgMax<<
471 TFile *fstream = cstream->GetFile();
473 array->Write("histos",TObject::kSingleKey);
477 void FitFit(Bool_t updateParam=kFALSE){
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];
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))");
498 //strFit+="sqrt(dri*ty)++";
499 //strFit+="sqrt(dri*tz)++";
500 //strFit+="sqrt(ty*tz)++";
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
504 cutPads[ipad]=Form("ipad==%d",ipad);
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());
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());
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());
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());
525 for (Int_t ipad=0; ipad<3; ipad++){
526 for (Int_t icorr=0; icorr<3; icorr++){
528 (*paramCl->fQNormCorr)(ipad+6,icorr) += paramTot[ipad][icorr+1];
529 (*paramCl->fQNormCorr)(ipad+9,icorr) += paramMax[ipad][icorr+1];
531 mat(ipad+0,icorr) = paramTot[ipad][icorr+1];
532 mat(ipad+3,icorr) = paramMax[ipad][icorr+1];
536 Float_t normShortTot;
541 Float_t normShortMax;
545 TH1D * dummyHist = new TH1D("dummyHist","absolute gain alignment of pads",100,0,10);
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.;
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.;
563 cout << "Tot: " << normShortTot << " " << normMedTot << " " << normLongTot << endl;
564 cout << "Max: " << normShortMax << " " << normMedMax << " " << normLongMax << endl;
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);
571 (*paramCl->fQNormCorr)(3,5) *= (normShortMax/normMaxMean);
572 (*paramCl->fQNormCorr)(4,5) *= (normMedMax/normMaxMean);
573 (*paramCl->fQNormCorr)(5,5) *= (normLongMax/normMaxMean);