AliRoot Core  a565103 (a565103)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CalibPID.C
Go to the documentation of this file.
1 
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 
68 AliTPCClusterParam * paramCl=0;
69 AliTPCcalibPID * pid =0;
70 TObjArray fitArr;
71 TTree * treeDump =0;
72 TTree * treeQtot;
73 TTree * treeQmax;
74 TTree * treeRatioQmax;
75 TTree * treeRatioQtot;
76 Int_t kmicolors[10]={1,2,3,4,6,7,8,9,10,11};
77 Int_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 
83 void Init(char* name="calibPID06"){
85 
86  TFile fcalib("CalibObjectsTrain2.root");
87  //TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib"); // old interface
88  pid = ( AliTPCcalibPID *) fcalib.Get(name);
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 
115 void StoreParam(char* localStorage = "local:///lustre/alice/akalweit/OCDB"){
116  Int_t runNumber = 0;
117  AliCDBMetaData *metaData= new AliCDBMetaData();
118  metaData->SetObjectClassName("AliTPCClusterParam");
119  metaData->SetResponsible("Alexander Kalweit");
120  metaData->SetBeamPeriod(1);
121  metaData->SetAliRootVersion("05-24-00");
122  metaData->SetComment("October runs calibration");
123  AliCDBId id1("TPC/Calib/ClusterParam", runNumber, AliCDBRunRange::Infinity());
124  AliCDBStorage *gStorage = AliCDBManager::Instance()->GetStorage(localStorage);
125  gStorage->Put(paramCl, id1, metaData);
126 
127 }
128 
129 
130 
131 
132 void 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 
144 void 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 
160 void 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 
173 void Fit(Bool_t updateParam=kFALSE){
175 
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;
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 
263  }
264  }
265 
266 
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 
279 void 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 
301 void 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);
305  pid->GetHistQtot()->GetAxis(0)->SetRangeUser(0.1,6.); // important adaption to be done here ...
306  pid->GetHistQmax()->GetAxis(6)->SetRangeUser(100,160);
307  pid->GetHistQmax()->GetAxis(0)->SetRangeUser(0.1,6); // important adaption to be done here ...
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 
475 void 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 
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
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  }
520 
521 
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 
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 
573 
574 }
575 
576 
577 
cutAll
TCut cutAll
Definition: CalibAlignKalman.C:104
treeQmax
TTree * treeQmax
Definition: CalibPID.C:73
printf
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
AliTPCClusterParam::fQNormCorr
TMatrixD * fQNormCorr
q norm correction for analytica correction
Definition: AliTPCClusterParam.h:178
AliTPCClusterParam.h
pid
AliTPCcalibPID * pid
Definition: CalibPID.C:69
fitArr
TObjArray fitArr
Definition: CalibPID.C:70
TObjArray
#define TObjArray
Definition: AliMUONReconstructor.h:19
Init
void Init(char *name="calibPID06")
Definition: CalibPID.C:83
ReadTrees
void ReadTrees()
Definition: CalibPID.C:160
mat
TMatrixD mat
Definition: AnalyzeLaser.C:9
Fit
void Fit(Bool_t updateParam=kFALSE)
Definition: CalibPID.C:173
AliTPCcalibDB.h
npoints
npoints
Definition: driftITSTPC.C:85
SetType
void SetType(Int_t type)
Definition: CalibPID.C:144
array
TObjArray * array
Definition: AnalyzeLaser.C:12
SetRange
void SetRange(Int_t index, Float_t min, Float_t max)
Definition: CalibPID.C:132
AliTPCClusterParam
TPC cluster error, shape and charge parameterization as function of drift length and inclination angl...
Definition: AliTPCClusterParam.h:21
chi2
Double_t chi2
Definition: AnalyzeLaser.C:7
kmicolors
Int_t kmicolors[10]
Definition: CalibPID.C:76
StoreParam
void StoreParam(char *localStorage="local:///lustre/alice/akalweit/OCDB")
Definition: CalibPID.C:115
FitFit
void FitFit(Bool_t updateParam=kFALSE)
Definition: CalibPID.C:475
treeQtot
TTree * treeQtot
Definition: CalibPID.C:72
treeRatioQtot
TTree * treeRatioQtot
Definition: CalibPID.C:75
toolkit
TStatToolkit toolkit
Definition: gainCalib.C:18
fitQdep
void fitQdep()
Definition: CalibPID.C:279
htemp
TObject * htemp
Definition: PlotSys.C:37
LookupHisto
void LookupHisto(Int_t minTracks=200, Float_t minp=20, Float_t maxp=10000)
Definition: CalibPID.C:301
kmimarkers
Int_t kmimarkers[10]
Definition: CalibPID.C:77
treeRatioQmax
TTree * treeRatioQmax
Definition: CalibPID.C:74
paramCl
AliTPCClusterParam * paramCl
Definition: CalibPID.C:68
treeDump
TTree * treeDump
Definition: CalibPID.C:71