TPC/CalibMacros/CalibPID.C

   1 /*
   2
   3 //
   4 // 1. dump information to the tree
   5 //
   6 gSystem->Load("libANALYSIS");
   7 gSystem->Load("libTPCcalib");
   8 gSystem->Load("libSTAT.so");
   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 // exit aliroot
  17
  18 //
  19 // 2. update the OCDB
  20 //
  21
  22 gSystem->Load("libANALYSIS");
  23 gSystem->Load("libTPCcalib");
  24 gSystem->Load("libSTAT.so");
  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 */
  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
  69 AliTPCClusterParam * paramCl=0;
  70 AliTPCcalibPID * pid =0;
  71 TObjArray fitArr;
  72 TTree * treeDump =0;
  73 TTree * treeQtot;
  74 TTree * treeQmax;
  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};
  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
  84 void Init(char* name="calibPID06"){
  85   //
  86   //
  87   TFile fcalib("CalibObjectsTrain2.root");
  88   //TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib"); // old interface
  89   pid = ( AliTPCcalibPID *) fcalib.Get(name);
  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
 116 void StoreParam(char* localStorage = "local:///lustre/alice/akalweit/OCDB"){
 117   Int_t runNumber = 0;
 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);
 127
 128 }
 129
 130
 131
 132
 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);
 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
 145 void 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
 161 void 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
 174 void Fit(Bool_t updateParam=kFALSE){
 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;
 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
 265     }
 266   }
 267
 268
 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
 281 void 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
 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 ...
 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
 477 void 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
 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);
 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   }
 522
 523
 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
 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
 575
 576 }
 577
 578
 579