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