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