TPC/CalibMacros/CalibPID.C

   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 /// ~~~
  41
  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"){
  84   ///
  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){
 174   /// align pads
 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