TPC/Upgrade/macros/AnaDelta.C

   1 #include <TStyle.h>
   2 #include <TROOT.h>
   3 #include <TAxis.h>
   4 #include <TF1.h>
   5 #include <TFile.h>
   6 #include <TH1.h>
   7 #include <THn.h>
   8 #include <TObjArray.h>
   9 #include <TObject.h>
  10 #include <TString.h>
  11 #include <TVectorT.h>
  12 #include <TCanvas.h>
  13 #include <TProfile2D.h>
  14 #include <TGraphErrors.h>
  15 #include <TTreeStream.h>
  16
  17 #include <AliExternalTrackParam.h>
  18 #include <AliTPCComposedCorrection.h>
  19 #include <AliTPCCorrectionLookupTable.h>
  20
  21 #include <AliToyMCEventGenerator.h>
  22
  23 /*
  24
  25 .L $ALICE_ROOT/TPC/Upgrade/macros/AnaDelta.C+g
  26
  27
  28 */
  29 TVectorD* MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax);
  30 TVectorD* MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax);
  31 TVectorD* MakeArbitraryBinning(const char* bins);
  32
  33 void DumpHn(THn *hn, TTreeSRedirector &stream);
  34 void AnaDeltaBase(TString file, TString outDir=".");
  35 void AnaDeltaTree(TString file, TString outFile="deltas_tree.root");
  36
  37 void AnaDelta(Int_t type, TString file, TString output="")
  38 {
  39   switch (type) {
  40     case 0:
  41       AnaDeltaBase(file,output);
  42       break;
  43     case 1:
  44       AnaDeltaTree(file,output);
  45       break;
  46   }
  47 }
  48
  49
  50 void AnaDeltaBase(TString file, TString outDir)
  51 {
  52   //
  53   //
  54   //
  55
  56   gStyle->SetOptFit();
  57
  58   TTreeSRedirector stream(Form("%s/deltas.root",outDir.Data()));
  59   gROOT->cd();
  60
  61   TFile f(file);
  62   THn *hn=(THn*)f.Get("hn");
  63
  64   DumpHn(hn, stream);
  65
  66   delete hn;
  67 }
  68
  69
  70 void AnaDeltaTree(TString file, TString outFile)
  71 {
  72   if (outFile.IsNull()) outFile="deltas_tree.root";
  73   TFile f(file);
  74   gROOT->cd();
  75   TTree *t = (TTree*)f.Get("delta");
  76   Float_t soneOverPt=0.;
  77   Float_t radius=0.;
  78   Float_t trackPhi=0.;
  79   Float_t trackY=0.;
  80   Float_t trackZ=0.;
  81   Float_t resRphi=0.;
  82   Double_t trackRes=0.;
  83   Float_t pointY=0.;
  84   Float_t pointZ=0.;
  85   Short_t npTRD=0.;
  86   Short_t event=0.;
  87
  88   t->SetBranchAddress("soneOverPt" , &soneOverPt);
  89   t->SetBranchAddress("r"          , &radius);
  90   t->SetBranchAddress("trackPhi"   , &trackPhi);
  91   t->SetBranchAddress("trackY"     , &trackY);
  92   t->SetBranchAddress("trackZ"     , &trackZ);
  93   t->SetBranchAddress("resRphi"    , &resRphi);
  94   t->SetBranchAddress("trackRes"   , &trackRes);
  95   t->SetBranchAddress("pointY"     , &pointY);
  96   t->SetBranchAddress("pointZ"     , &pointZ);
  97   t->SetBranchAddress("npTRD"      , &npTRD);
  98   t->SetBranchAddress("event"      , &event);
  99
 100   // make binning
 101   TVectorD *vR   = MakeLinBinning(10,86.,250.);
 102   TVectorD *vPhi = MakeLinBinning(18*8,0.,2*TMath::Pi());
 103   TVectorD *vZ   = MakeLinBinning(50,-250.,250.);
 104
 105   const Int_t nbins=4;
 106   Int_t bins[nbins]    = {vR->GetNrows()-1, vPhi->GetNrows()-1, vZ->GetNrows()-1, 80};
 107 //   Int_t bins[nbins]    = {16, 18*5, 50, 80};
 108   Double_t xmin[nbins] = {86. , 0.,           -250., -2.};
 109   Double_t xmax[nbins] = {250., 2*TMath::Pi(), 250.,  2.};
 110   THnF *hn = new THnF("hn", "hn", nbins, bins, xmin, xmax);
 111
 112   hn->GetAxis(0)->Set(vR  ->GetNrows()-1, vR  ->GetMatrixArray());
 113   hn->GetAxis(1)->Set(vPhi->GetNrows()-1, vPhi->GetMatrixArray());
 114   hn->GetAxis(2)->Set(vZ  ->GetNrows()-1, vZ  ->GetMatrixArray());
 115
 116   hn->GetAxis(0)->SetNameTitle("r","r (cm)");
 117   hn->GetAxis(1)->SetNameTitle("phi","#varphi");
 118   hn->GetAxis(2)->SetNameTitle("z","z (cm)");
 119   hn->GetAxis(3)->SetNameTitle("drphi","#Delta(r#varphi)");
 120
 121   for (Int_t iev=0; iev<t->GetEntries(); ++iev) {
 122     t->GetEntry(iev);
 123
 124     // cuts
 125     // -- on trd
 126     if (npTRD<2) continue;
 127     Double_t pt=1./TMath::Abs(soneOverPt);
 128     if (pt<0.8) continue;
 129
 130     Float_t resRphiRandom = resRphi*trackRes;
 131     Float_t deviation     = pointY-(trackY+resRphiRandom);
 132
 133     Double_t xx[4]={radius, trackPhi, trackZ ,deviation};
 134     hn->Fill(xx);
 135   }
 136
 137   // do fits and fill tree
 138   TTreeSRedirector stream(outFile.Data());
 139   gROOT->cd();
 140
 141   DumpHn(hn, stream);
 142
 143   stream.GetFile()->cd();
 144   hn->Write();
 145
 146   delete hn;
 147   delete vR;
 148   delete vPhi;
 149   delete vZ;
 150 }
 151
 152
 153 void AnaDeltaTree2(TString file/*, TString outDir="."*/)
 154 {
 155   //
 156   // NOTE: not finished
 157   //
 158   TFile f(file);
 159   gROOT->cd();
 160   TTree *t = (TTree*)f.Get("delta");
 161   Float_t soneOverPt=0.;
 162   Float_t radius=0.;
 163   Float_t trackPhi=0.;
 164   Float_t trackY=0.;
 165   Float_t trackZ=0.;
 166   Float_t resRphi=0.;
 167   Float_t trackRes=0.;
 168   Float_t pointY=0.;
 169   Float_t pointZ=0.;
 170   Float_t npTRD=0.;
 171   Float_t event=0.;
 172
 173   t->SetBranchAddress("soneOverPt" , &soneOverPt);
 174   t->SetBranchAddress("r"          , &radius);
 175   t->SetBranchAddress("trackPhi"   , &trackPhi);
 176   t->SetBranchAddress("trackY"     , &trackY);
 177   t->SetBranchAddress("trackZ"     , &trackZ);
 178   t->SetBranchAddress("resRphi"    , &resRphi);
 179   t->SetBranchAddress("trackRes"   , &trackRes);
 180   t->SetBranchAddress("pointY"     , &pointY);
 181   t->SetBranchAddress("pointZ"     , &pointZ);
 182   t->SetBranchAddress("npTRD"      , &npTRD);
 183   t->SetBranchAddress("event"      , &event);
 184
 185   // make binning
 186   TVectorD *vZ   = MakeLinBinning(50,-250.,250.);
 187   TVectorD *vPhi = MakeLinBinning(18*8,0.,TMath::Pi());
 188   TVectorD *vR   = MakeLinBinning(16,86.,250.);
 189
 190   TObjArray arrZ(vZ->GetNrows()-1);
 191   arrZ.SetOwner();
 192
 193   for (Int_t iev=0; iev<t->GetEntries(); ++iev) {
 194     t->GetEntry(iev);
 195
 196     // cuts
 197     // -- on trd
 198     if (npTRD<2) continue;
 199
 200     Float_t resRphiRandom=resRphi*trackRes;
 201
 202     Int_t binZ   = TMath::BinarySearch(vZ->GetNrows(),vZ->GetMatrixArray(),(Double_t)trackZ);
 203     Int_t binPhi = TMath::BinarySearch(vPhi->GetNrows(),vPhi->GetMatrixArray(),(Double_t)trackPhi);
 204     Int_t binR   = TMath::BinarySearch(vR->GetNrows(),vR->GetMatrixArray(),(Double_t)radius);
 205
 206     if (binZ<0)   binZ=0;
 207     if (binPhi<0) binPhi=0;
 208     if (binR<0)   binR=0;
 209
 210     TObjArray *arrPhi=(TObjArray*)arrZ.UncheckedAt(binZ);
 211     if (!arrPhi) {
 212       arrPhi=new TObjArray(vPhi->GetNrows()-1);
 213       arrZ.AddAt(arrPhi,binZ);
 214     }
 215
 216     TObjArray *arrR=(TObjArray*)arrPhi->UncheckedAt(binPhi);
 217     if (!arrR) {
 218       arrR=new TObjArray(vR->GetNrows()-1);
 219       arrPhi->AddAt(arrR,binPhi);
 220     }
 221
 222     TH1S *h = (TH1S*)arrR->UncheckedAt(binR);
 223     if (!h) {
 224       h = new TH1S(Form("h_%02d_%02d_%d02",binZ, binPhi, binR),
 225                    Form("z,phi,r: %02d,%02d,%d02; #Delta r#phi (cm)",binZ, binPhi, binR),
 226                    80, -2., 2.);
 227       arrR->AddAt(h, binR);
 228     }
 229
 230     h->Fill(trackY+resRphiRandom-pointY);
 231   }
 232
 233   // do fits and fill tree
 234 }
 235
 236 void AnaDeltaResiduals(TString fluctuationMap, TString averageMap, TString outFile="deltas_residuals.root")
 237 {
 238   //
 239   //
 240   //
 241
 242   TFile fFluct(fluctuationMap);
 243   AliTPCCorrectionLookupTable *corrFluct = (AliTPCCorrectionLookupTable*)fFluct.Get("map");
 244   fFluct.Close();
 245
 246   TFile fAverage(averageMap);
 247   AliTPCCorrectionLookupTable *corrAverage = (AliTPCCorrectionLookupTable*)fAverage.Get("map");
 248   fAverage.Close();
 249
 250 //   TObjArray *arrMaps = new TObjArray(2);
 251 //   arrMaps->Add(corrAverage); // correction with the average Map
 252 //   arrMaps->Add(corrFluct);   // distortion with the fluctuation Map
 253
 254   // create the composed correction
 255   // if the weight are set to +1 and -1, the first map will be responsible for the distortions
 256   // The second map for the corrections
 257   // !!!!! In AliTPCComposedCorrection::GetDistortion MakeInverseIterator is called !!!!
 258   // for this reason we have to add the maps in the wrong order
 259
 260 //   AliTPCComposedCorrection *residualDistortion = new AliTPCComposedCorrection(arrMaps, AliTPCComposedCorrection::kQueueResidual);
 261   Float_t dummy=0;
 262 //   TVectorD weights(2);
 263 //   weights(0)=+1.;
 264 //   weights(1)=-AliToyMCEventGenerator::GetSCScalingFactor(corrFluct, corrAverage,dummy);
 265 //   residualDistortion->SetWeights(&weights);
 266
 267   corrAverage->SetCorrScaleFactor(AliToyMCEventGenerator::GetSCScalingFactor(corrFluct, corrAverage,dummy));
 268
 269   TVectorD *vR   = MakeLinBinning(10,86.,250.);
 270   TVectorD *vPhi = MakeLinBinning(18*8,0.,2*TMath::Pi());
 271   TVectorD *vZ   = MakeLinBinning(50,-250.,250.);
 272
 273   const Int_t nbins=4;
 274   Int_t bins[nbins]    = {vR->GetNrows()-1, vPhi->GetNrows()-1, vZ->GetNrows()-1, 80};
 275   //   Int_t bins[nbins]    = {16, 18*5, 50, 80};
 276   Double_t xmin[nbins] = {86. , 0.,           -250., -2.};
 277   Double_t xmax[nbins] = {250., 2*TMath::Pi(), 250.,  2.};
 278   THnF *hn = new THnF("hn", "hn", nbins, bins, xmin, xmax);
 279
 280   hn->GetAxis(0)->Set(vR  ->GetNrows()-1, vR  ->GetMatrixArray());
 281   hn->GetAxis(1)->Set(vPhi->GetNrows()-1, vPhi->GetMatrixArray());
 282   hn->GetAxis(2)->Set(vZ  ->GetNrows()-1, vZ  ->GetMatrixArray());
 283
 284   AliExternalTrackParam vv;
 285
 286   for (Float_t iz=-245; iz<=245; iz+=2) {
 287     Short_t roc=(iz>=0)?0:18;
 288     for (Float_t ir=86; ir<250; ir+=1) {
 289       for (Float_t iphi=0; iphi<TMath::TwoPi(); iphi+=0.5*TMath::DegToRad()){
 290         Float_t x=ir*(Float_t)TMath::Cos(iphi);
 291         Float_t y=ir*(Float_t)TMath::Sin(iphi);
 292         Float_t x3[3]    = {x,y,iz};
 293         Float_t x3dc[3]    = {x,y,iz};
 294         Float_t dx3[3]   = {0.,0.,0.};
 295 //         residualDistortion->GetDistortion(x3,roc,dx3);
 296         corrFluct->DistortPoint(x3dc,roc);
 297         corrAverage->CorrectPoint(x3dc,roc);
 298         dx3[0]=x3dc[0]-x3[0];
 299         dx3[1]=x3dc[1]-x3[1];
 300         dx3[2]=x3dc[2]-x3[2];
 301
 302         Double_t ddx3[3]={dx3[0], dx3[1], dx3[2]};
 303         vv.Global2LocalPosition(ddx3,iphi);
 304
 305         Double_t xx[4]={ir, iphi, iz ,ddx3[1]};
 306         hn->Fill(xx);
 307
 308       }
 309     }
 310   }
 311
 312   TTreeSRedirector stream(outFile.Data());
 313   gROOT->cd();
 314
 315   DumpHn(hn, stream);
 316
 317   stream.GetFile()->cd();
 318   hn->Write();
 319
 320   delete hn;
 321   delete vR;
 322   delete vPhi;
 323   delete vZ;
 324
 325 //   delete residualDistortion;
 326 }
 327
 328 void DumpHn(THn *hn, TTreeSRedirector &stream)
 329 {
 330   TAxis *ar   = hn->GetAxis(0);
 331   TAxis *aphi = hn->GetAxis(1);
 332   TAxis *az   = hn->GetAxis(2);
 333
 334   // output Hn
 335   const Int_t nbins=3;
 336   Int_t bins[nbins]    = {1,1,1};
 337   Double_t xmin[nbins] = {0.,0.,0.};
 338   Double_t xmax[nbins] = {1.,1.,1.};
 339   THnF hnRes("hnRes", "hnRes", nbins, bins, xmin, xmax);
 340
 341   ar  ->Copy(*hnRes.GetAxis(0));
 342   aphi->Copy(*hnRes.GetAxis(1));
 343   az  ->Copy(*hnRes.GetAxis(2));
 344
 345
 346   for (Int_t iz=0; iz<az->GetNbins(); ++iz) {
 347     az->SetRange(iz+1,iz+1);
 348     TObjArray arrFits;
 349     arrFits.SetName(Form("z_%02d",iz));
 350     arrFits.SetOwner();
 351
 352     for (Int_t ir=0; ir<ar->GetNbins(); ++ir) {
 353       ar->SetRange(ir+1,ir+1);
 354       for (Int_t iphi=0; iphi<aphi->GetNbins(); ++iphi) {
 355         aphi->SetRange(iphi+1,iphi+1);
 356
 357         Float_t cr       = 0.;
 358         Float_t cphi     = 0.;
 359         Float_t cz       = 0.;
 360         Float_t mean     = 0.;
 361         Float_t meanErr  = 0.;
 362         Float_t sigma    = 0.;
 363         Float_t sigmaErr = 0.;
 364         Int_t   entries  = 0.;
 365         Float_t chi2ndf  = 0.;
 366         Float_t mean2    = 0.;
 367         Float_t meanErr2 = 0.;
 368         Float_t rms2     = 0.;
 369         Float_t rmsErr2  = 0.;
 370
 371         TH1 *hProj = hn->Projection(3);
 372         if (hProj->GetEntries()>1) {
 373           TF1 fg("fg","gaus",-2,2);
 374           cr   = ar->GetBinCenter(ir+1);
 375           cphi = aphi->GetBinCenter(iphi+1);
 376           cz   = az->GetBinCenter(iz+1);
 377           hProj->SetNameTitle(Form("h_%02d_%02d_%02d",iz, iphi, ir),
 378           Form("z,phi,r: %02d,%02d,%02d (%.2f, %.2f, %.2f)",iz,iphi,ir, cz, cphi, cr )
 379           );
 380           hProj->Fit(&fg,"LMQR");
 381           arrFits.Add(hProj);
 382
 383           mean     = fg.GetParameter(1);
 384           meanErr  = fg.GetParError(1);
 385           sigma    = fg.GetParameter(2);
 386           sigmaErr = fg.GetParError(2);
 387           entries  = hProj->GetEntries();
 388           chi2ndf  = fg.GetChisquare()/fg.GetNDF();
 389           mean2    = hProj->GetMean();
 390           meanErr2 = hProj->GetMeanError();
 391           rms2     = hProj->GetRMS();
 392           rmsErr2  = hProj->GetRMSError();
 393         } else {
 394           delete hProj;
 395         }
 396
 397         stream << "d" <<
 398         "ir="          << ir       <<
 399         "iphi="        << iphi     <<
 400         "iz="          << iz       <<
 401         "cr="          << cr       <<
 402         "cphi="        << cphi     <<
 403         "cz="          << cz       <<
 404         "mean="        << mean     <<
 405         "meanErr="     << meanErr  <<
 406         "sigma="       << sigma    <<
 407         "sigmaErr="    << sigmaErr <<
 408         "histMean="    << mean2    <<
 409         "histMeanErr=" << meanErr2 <<
 410         "histRMS="     << rms2    <<
 411         "histRMSErr="  << rmsErr2 <<
 412         "entries="     << entries  <<
 413         "chi2ndf="     << chi2ndf  <<
 414         "\n";
 415
 416 //        Double_t x[nbins]={cr, cphi, cz};
 417 //        if (meanErr<0.3) hnRes.Fill(x,mean);
 418       }
 419     }
 420     stream.GetFile()->cd();
 421     arrFits.Write(0x0,TObject::kSingleKey);
 422     gROOT->cd();
 423   }
 424
 425   stream.GetFile()->cd();
 426   hnRes.Write();
 427   gROOT->cd();
 428 }
 429
 430 //______________________________________________________________________________
 431 TVectorD* MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
 432 {
 433   //
 434   // Make logarithmic binning
 435   // the user has to delete the array afterwards!!!
 436   //
 437
 438   //check limits
 439   if (xmin<1e-20 || xmax<1e-20){
 440     printf("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
 441     return MakeLinBinning(nbinsX, xmin, xmax);
 442   }
 443   if (xmax<xmin){
 444     Double_t tmp=xmin;
 445     xmin=xmax;
 446     xmax=tmp;
 447   }
 448   TVectorD *binLim=new TVectorD(nbinsX+1);
 449   Double_t first=xmin;
 450   Double_t last=xmax;
 451   Double_t expMax=TMath::Log(last/first);
 452   for (Int_t i=0; i<nbinsX+1; ++i){
 453     (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
 454   }
 455   return binLim;
 456 }
 457
 458 //______________________________________________________________________________
 459 TVectorD* MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
 460 {
 461   //
 462   // Make linear binning
 463   // the user has to delete the array afterwards!!!
 464   //
 465   if (xmax<xmin){
 466     Double_t tmp=xmin;
 467     xmin=xmax;
 468     xmax=tmp;
 469   }
 470   TVectorD *binLim=new TVectorD(nbinsX+1);
 471   Double_t first=xmin;
 472   Double_t last=xmax;
 473   Double_t binWidth=(last-first)/nbinsX;
 474   for (Int_t i=0; i<nbinsX+1; ++i){
 475     (*binLim)[i]=first+binWidth*(Double_t)i;
 476   }
 477   return binLim;
 478 }
 479
 480 //_____________________________________________________________________________
 481 TVectorD* MakeArbitraryBinning(const char* bins)
 482 {
 483   //
 484   // Make arbitrary binning, bins separated by a ','
 485   //
 486   TString limits(bins);
 487   if (limits.IsNull()){
 488     printf("Bin Limit string is empty, cannot add the variable");
 489     return 0x0;
 490   }
 491
 492   TObjArray *arr=limits.Tokenize(",");
 493   Int_t nLimits=arr->GetEntries();
 494   if (nLimits<2){
 495     printf("Need at leas 2 bin limits, cannot add the variable");
 496     delete arr;
 497     return 0x0;
 498   }
 499
 500   TVectorD *binLimits=new TVectorD(nLimits);
 501   for (Int_t iLim=0; iLim<nLimits; ++iLim){
 502     (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
 503   }
 504
 505   delete arr;
 506   return binLimits;
 507 }
 508
 509 void PlotFromTree(TTree *d, TString outDir=".")
 510 {
 511   TCanvas *c=new TCanvas;
 512   gStyle->SetOptStat(0);
 513   d->SetMarkerStyle(20);
 514   d->SetMarkerSize(1);
 515
 516   TProfile2D pRZ("pRZ",";z (cm); r(cm)",50,-250,250,10,85,250);
 517   d->Draw("entries:cr:cz>>pRZ","","profcolz");
 518   pRZ.GetZaxis()->UnZoom();
 519   c->SaveAs(Form("%s/entries_average.png",outDir.Data()));
 520   d->Draw("entries:cr:cz>>pRZ","iphi==2","profcolz");
 521   c->SaveAs(Form("%s/entries_onePhi.png",outDir.Data()));
 522
 523   pRZ.SetMaximum(0.04);
 524   d->Draw("meanErr:cr:cz>>pRZ","","profcolz");
 525   c->SaveAs(Form("%s/meanErr_average.png",outDir.Data()));
 526   d->Draw("meanErr:cr:cz>>pRZ","iphi==2","profcolz");
 527   c->SaveAs(Form("%s/meanErr_onePhi.png",outDir.Data()));
 528
 529
 530   d->Draw("mean:cphi:cr","iz==25","colz");
 531   c->SaveAs(Form("%s/mean_oneZ_phi_allR.png",outDir.Data()));
 532   d->Draw("mean:meanErr:cphi","iz==25&&ir==2","goff");
 533   TGraphErrors *grmean_phi=new TGraphErrors(d->GetSelectedRows(),d->GetV3(),d->GetV1(),0,d->GetV2());
 534   grmean_phi->SetTitle(";#varphi;#LT#Delta r#varphi#GT");
 535   grmean_phi->SetMarkerStyle(20);
 536   grmean_phi->SetMarkerSize(1);
 537   grmean_phi->Draw("ap");
 538   c->SaveAs(Form("%s/mean_oneZ_phi_oneR.png",outDir.Data()));
 539
 540   d->Draw("mean:cr:cphi","iz==25","colz");
 541   c->SaveAs(Form("%s/mean_oneZ_r_allPhi.png",outDir.Data()));
 542
 543   d->Draw("mean:meanErr:cr","iz==25&&iphi==2","goff");
 544   TGraphErrors *grmean_r=new TGraphErrors(d->GetSelectedRows(),d->GetV3(),d->GetV1(),0,d->GetV2());
 545   grmean_r->SetTitle(";r (cm);#LT#Delta r#varphi#GT");
 546   grmean_r->SetMarkerStyle(20);
 547   grmean_r->SetMarkerSize(1);
 548   grmean_r->Draw("ap");
 549   c->SaveAs(Form("%s/mean_oneZ_r_onePhi.png",outDir.Data()));
 550
 551
 552   d->Draw("meanErr:cphi:cr","iz==25","colz");
 553   c->SaveAs(Form("%s/meanErr_oneZ_phi_allR.png",outDir.Data()));
 554   d->Draw("meanErr:cphi","iz==25&&ir==2");
 555   c->SaveAs(Form("%s/meanErr_oneZ_phi_oneR.png",outDir.Data()));
 556
 557   d->Draw("meanErr:cr:cphi","iz==25","colz");
 558   c->SaveAs(Form("%s/meanErr_oneZ_r_allPhi.png",outDir.Data()));
 559
 560   d->Draw("meanErr:cr","iz==25&&iphi==2");
 561   c->SaveAs(Form("%s/meanErr_oneZ_r_onePhi.png",outDir.Data()));
 562
 563 }
 564
 565