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 /*
  18
  19 .L $ALICE_ROOT/TPC/Upgrade/macros/AnaDelta.C+g
  20
  21
  22 */
  23 TVectorD* MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax);
  24 TVectorD* MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax);
  25 TVectorD* MakeArbitraryBinning(const char* bins);
  26
  27 void DumpHn(THn *hn, TTreeSRedirector &stream);
  28
  29 void AnaDelta(TString file, TString outDir=".")
  30 {
  31   //
  32   //
  33   //
  34
  35   gStyle->SetOptFit();
  36
  37   TTreeSRedirector stream(Form("%s/deltas.root",outDir.Data()));
  38   gROOT->cd();
  39
  40   TFile f(file);
  41   THn *hn=(THn*)f.Get("hn");
  42
  43   DumpHn(hn, stream);
  44
  45   delete hn;
  46 }
  47
  48
  49 void AnaDeltaTree(TString file, TString outFile="deltas_tree.root")
  50 {
  51   TFile f(file);
  52   gROOT->cd();
  53   TTree *t = (TTree*)f.Get("delta");
  54   Float_t soneOverPt=0.;
  55   Float_t radius=0.;
  56   Float_t trackPhi=0.;
  57   Float_t trackY=0.;
  58   Float_t trackZ=0.;
  59   Float_t resRphi=0.;
  60   Double_t trackRes=0.;
  61   Float_t pointY=0.;
  62   Float_t pointZ=0.;
  63   Short_t npTRD=0.;
  64   Short_t event=0.;
  65
  66   t->SetBranchAddress("soneOverPt" , &soneOverPt);
  67   t->SetBranchAddress("r"          , &radius);
  68   t->SetBranchAddress("trackPhi"   , &trackPhi);
  69   t->SetBranchAddress("trackY"     , &trackY);
  70   t->SetBranchAddress("trackZ"     , &trackZ);
  71   t->SetBranchAddress("resRphi"    , &resRphi);
  72   t->SetBranchAddress("trackRes"   , &trackRes);
  73   t->SetBranchAddress("pointY"     , &pointY);
  74   t->SetBranchAddress("pointZ"     , &pointZ);
  75   t->SetBranchAddress("npTRD"      , &npTRD);
  76   t->SetBranchAddress("event"      , &event);
  77
  78   // make binning
  79   TVectorD *vR   = MakeLinBinning(10,86.,250.);
  80   TVectorD *vPhi = MakeLinBinning(18*8,0.,2*TMath::Pi());
  81   TVectorD *vZ   = MakeLinBinning(50,-250.,250.);
  82
  83   const Int_t nbins=4;
  84   Int_t bins[nbins]    = {vR->GetNrows()-1, vPhi->GetNrows()-1, vZ->GetNrows()-1, 80};
  85 //   Int_t bins[nbins]    = {16, 18*5, 50, 80};
  86   Double_t xmin[nbins] = {86. , 0.,           -250., -2.};
  87   Double_t xmax[nbins] = {250., 2*TMath::Pi(), 250.,  2.};
  88   THnF *hn = new THnF("hn", "hn", nbins, bins, xmin, xmax);
  89
  90   hn->GetAxis(0)->Set(vR  ->GetNrows()-1, vR  ->GetMatrixArray());
  91   hn->GetAxis(1)->Set(vPhi->GetNrows()-1, vPhi->GetMatrixArray());
  92   hn->GetAxis(2)->Set(vZ  ->GetNrows()-1, vZ  ->GetMatrixArray());
  93
  94
  95   for (Int_t iev=0; iev<t->GetEntries(); ++iev) {
  96     t->GetEntry(iev);
  97
  98     // cuts
  99     // -- on trd
 100     if (npTRD<2) continue;
 101     Double_t pt=1./TMath::Abs(soneOverPt);
 102     if (pt<0.8) continue;
 103
 104     Float_t resRphiRandom = resRphi*trackRes;
 105     Float_t deviation     = trackY+resRphiRandom-pointY;
 106
 107     Double_t xx[4]={radius, trackPhi, trackZ ,deviation};
 108     hn->Fill(xx);
 109   }
 110
 111   // do fits and fill tree
 112   TTreeSRedirector stream(outFile.Data());
 113   gROOT->cd();
 114
 115   DumpHn(hn, stream);
 116
 117   stream.GetFile()->cd();
 118   hn->Write();
 119
 120   delete hn;
 121   delete vR;
 122   delete vPhi;
 123   delete vZ;
 124 }
 125
 126
 127 void AnaDeltaTree2(TString file/*, TString outDir="."*/)
 128 {
 129   //
 130   // NOTE: not finished
 131   //
 132   TFile f(file);
 133   gROOT->cd();
 134   TTree *t = (TTree*)f.Get("delta");
 135   Float_t soneOverPt=0.;
 136   Float_t radius=0.;
 137   Float_t trackPhi=0.;
 138   Float_t trackY=0.;
 139   Float_t trackZ=0.;
 140   Float_t resRphi=0.;
 141   Float_t trackRes=0.;
 142   Float_t pointY=0.;
 143   Float_t pointZ=0.;
 144   Float_t npTRD=0.;
 145   Float_t event=0.;
 146
 147   t->SetBranchAddress("soneOverPt" , &soneOverPt);
 148   t->SetBranchAddress("r"          , &radius);
 149   t->SetBranchAddress("trackPhi"   , &trackPhi);
 150   t->SetBranchAddress("trackY"     , &trackY);
 151   t->SetBranchAddress("trackZ"     , &trackZ);
 152   t->SetBranchAddress("resRphi"    , &resRphi);
 153   t->SetBranchAddress("trackRes"   , &trackRes);
 154   t->SetBranchAddress("pointY"     , &pointY);
 155   t->SetBranchAddress("pointZ"     , &pointZ);
 156   t->SetBranchAddress("npTRD"      , &npTRD);
 157   t->SetBranchAddress("event"      , &event);
 158
 159   // make binning
 160   TVectorD *vZ   = MakeLinBinning(50,-250.,250.);
 161   TVectorD *vPhi = MakeLinBinning(18*8,0.,TMath::Pi());
 162   TVectorD *vR   = MakeLinBinning(16,86.,250.);
 163
 164   TObjArray arrZ(vZ->GetNrows()-1);
 165   arrZ.SetOwner();
 166
 167   for (Int_t iev=0; iev<t->GetEntries(); ++iev) {
 168     t->GetEntry(iev);
 169
 170     // cuts
 171     // -- on trd
 172     if (npTRD<2) continue;
 173
 174     Float_t resRphiRandom=resRphi*trackRes;
 175
 176     Int_t binZ   = TMath::BinarySearch(vZ->GetNrows(),vZ->GetMatrixArray(),(Double_t)trackZ);
 177     Int_t binPhi = TMath::BinarySearch(vPhi->GetNrows(),vPhi->GetMatrixArray(),(Double_t)trackPhi);
 178     Int_t binR   = TMath::BinarySearch(vR->GetNrows(),vR->GetMatrixArray(),(Double_t)radius);
 179
 180     if (binZ<0)   binZ=0;
 181     if (binPhi<0) binPhi=0;
 182     if (binR<0)   binR=0;
 183
 184     TObjArray *arrPhi=(TObjArray*)arrZ.UncheckedAt(binZ);
 185     if (!arrPhi) {
 186       arrPhi=new TObjArray(vPhi->GetNrows()-1);
 187       arrZ.AddAt(arrPhi,binZ);
 188     }
 189
 190     TObjArray *arrR=(TObjArray*)arrPhi->UncheckedAt(binPhi);
 191     if (!arrR) {
 192       arrR=new TObjArray(vR->GetNrows()-1);
 193       arrPhi->AddAt(arrR,binPhi);
 194     }
 195
 196     TH1S *h = (TH1S*)arrR->UncheckedAt(binR);
 197     if (!h) {
 198       h = new TH1S(Form("h_%02d_%02d_%d02",binZ, binPhi, binR),
 199                    Form("z,phi,r: %02d,%02d,%d02; #Delta r#phi (cm)",binZ, binPhi, binR),
 200                    80, -2., 2.);
 201       arrR->AddAt(h, binR);
 202     }
 203
 204     h->Fill(trackY+resRphiRandom-pointY);
 205   }
 206
 207   // do fits and fill tree
 208 }
 209
 210
 211 void DumpHn(THn *hn, TTreeSRedirector &stream)
 212 {
 213   TAxis *ar   = hn->GetAxis(0);
 214   TAxis *aphi = hn->GetAxis(1);
 215   TAxis *az   = hn->GetAxis(2);
 216
 217
 218   for (Int_t iz=0; iz<az->GetNbins(); ++iz) {
 219     az->SetRange(iz+1,iz+1);
 220     TObjArray arrFits;
 221     arrFits.SetName(Form("z_%02d",iz));
 222     arrFits.SetOwner();
 223
 224     for (Int_t ir=0; ir<ar->GetNbins(); ++ir) {
 225       ar->SetRange(ir+1,ir+1);
 226       for (Int_t iphi=0; iphi<aphi->GetNbins(); ++iphi) {
 227         aphi->SetRange(iphi+1,iphi+1);
 228
 229         TH1 *hProj = hn->Projection(3);
 230         if (hProj->GetEntries()<1) {
 231           delete hProj;
 232           continue;
 233         }
 234
 235         TF1 fg("fg","gaus",-2,2);
 236         Double_t cr   = ar->GetBinCenter(ir+1);
 237         Double_t cphi = aphi->GetBinCenter(iphi+1);
 238         Double_t cz   = az->GetBinCenter(iz+1);
 239         hProj->SetNameTitle(Form("h_%02d_%02d_%d02",iz+1, iphi+1, ir+1),
 240                             Form("z,phi,r: %02d,%02d,%d02 (%.2f, %.2f, %.2f)",iz+1,iphi+1,ir+1, cr, cphi, cz )
 241         );
 242         hProj->Fit(&fg,"QR");
 243         arrFits.Add(hProj);
 244
 245         Double_t mean     = fg.GetParameter(1);
 246         Double_t meanErr  = fg.GetParError(1);
 247         Double_t sigma    = fg.GetParameter(2);
 248         Double_t sigmaErr = fg.GetParError(2);
 249         Int_t    entries  = hProj->GetEntries();
 250         Double_t chi2ndf  = fg.GetChisquare()/fg.GetNDF();
 251
 252         stream << "d" <<
 253         "ir="        << ir       <<
 254         "iphi="      << iphi     <<
 255         "iz="        << iz       <<
 256         "cr="        << cr       <<
 257         "cphi="      << cphi     <<
 258         "cz="        << cz       <<
 259         "mean="      << mean     <<
 260         "meanErr="   << meanErr  <<
 261         "sigma="     << sigma    <<
 262         "sigmaErr="  << sigmaErr <<
 263         "entries="   << entries  <<
 264         "chi2ndf="   << chi2ndf  <<
 265         "\n";
 266       }
 267     }
 268     stream.GetFile()->cd();
 269     arrFits.Write(0x0,TObject::kSingleKey);
 270     gROOT->cd();
 271   }
 272
 273 }
 274
 275 //______________________________________________________________________________
 276 TVectorD* MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
 277 {
 278   //
 279   // Make logarithmic binning
 280   // the user has to delete the array afterwards!!!
 281   //
 282
 283   //check limits
 284   if (xmin<1e-20 || xmax<1e-20){
 285     printf("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
 286     return MakeLinBinning(nbinsX, xmin, xmax);
 287   }
 288   if (xmax<xmin){
 289     Double_t tmp=xmin;
 290     xmin=xmax;
 291     xmax=tmp;
 292   }
 293   TVectorD *binLim=new TVectorD(nbinsX+1);
 294   Double_t first=xmin;
 295   Double_t last=xmax;
 296   Double_t expMax=TMath::Log(last/first);
 297   for (Int_t i=0; i<nbinsX+1; ++i){
 298     (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
 299   }
 300   return binLim;
 301 }
 302
 303 //______________________________________________________________________________
 304 TVectorD* MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
 305 {
 306   //
 307   // Make linear binning
 308   // the user has to delete the array afterwards!!!
 309   //
 310   if (xmax<xmin){
 311     Double_t tmp=xmin;
 312     xmin=xmax;
 313     xmax=tmp;
 314   }
 315   TVectorD *binLim=new TVectorD(nbinsX+1);
 316   Double_t first=xmin;
 317   Double_t last=xmax;
 318   Double_t binWidth=(last-first)/nbinsX;
 319   for (Int_t i=0; i<nbinsX+1; ++i){
 320     (*binLim)[i]=first+binWidth*(Double_t)i;
 321   }
 322   return binLim;
 323 }
 324
 325 //_____________________________________________________________________________
 326 TVectorD* MakeArbitraryBinning(const char* bins)
 327 {
 328   //
 329   // Make arbitrary binning, bins separated by a ','
 330   //
 331   TString limits(bins);
 332   if (limits.IsNull()){
 333     printf("Bin Limit string is empty, cannot add the variable");
 334     return 0x0;
 335   }
 336
 337   TObjArray *arr=limits.Tokenize(",");
 338   Int_t nLimits=arr->GetEntries();
 339   if (nLimits<2){
 340     printf("Need at leas 2 bin limits, cannot add the variable");
 341     delete arr;
 342     return 0x0;
 343   }
 344
 345   TVectorD *binLimits=new TVectorD(nLimits);
 346   for (Int_t iLim=0; iLim<nLimits; ++iLim){
 347     (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
 348   }
 349
 350   delete arr;
 351   return binLimits;
 352 }
 353
 354 void PlotFromTree(TTree *d, TString outDir=".")
 355 {
 356   TCanvas *c=new TCanvas;
 357   gStyle->SetOptStat(0);
 358   d->SetMarkerStyle(20);
 359   d->SetMarkerSize(1);
 360
 361   TProfile2D pRZ("pRZ",";z (cm); r(cm)",50,-250,250,10,85,250);
 362   d->Draw("entries:cr:cz>>pRZ","","profcolz");
 363   pRZ.GetZaxis()->UnZoom();
 364   c->SaveAs(Form("%s/entries_average.png",outDir.Data()));
 365   d->Draw("entries:cr:cz>>pRZ","iphi==2","profcolz");
 366   c->SaveAs(Form("%s/entries_onePhi.png",outDir.Data()));
 367
 368   pRZ.SetMaximum(0.04);
 369   d->Draw("meanErr:cr:cz>>pRZ","","profcolz");
 370   c->SaveAs(Form("%s/meanErr_average.png",outDir.Data()));
 371   d->Draw("meanErr:cr:cz>>pRZ","iphi==2","profcolz");
 372   c->SaveAs(Form("%s/meanErr_onePhi.png",outDir.Data()));
 373
 374
 375   d->Draw("mean:cphi:cr","iz==25","colz");
 376   c->SaveAs(Form("%s/mean_oneZ_phi_allR.png",outDir.Data()));
 377   d->Draw("mean:meanErr:cphi","iz==25&&ir==2","goff");
 378   TGraphErrors *grmean_phi=new TGraphErrors(d->GetSelectedRows(),d->GetV3(),d->GetV1(),0,d->GetV2());
 379   grmean_phi->SetTitle(";#varphi;#LT#Delta r#varphi#GT");
 380   grmean_phi->SetMarkerStyle(20);
 381   grmean_phi->SetMarkerSize(1);
 382   grmean_phi->Draw("ap");
 383   c->SaveAs(Form("%s/mean_oneZ_phi_oneR.png",outDir.Data()));
 384
 385   d->Draw("mean:cr:cphi","iz==25","colz");
 386   c->SaveAs(Form("%s/mean_oneZ_r_allPhi.png",outDir.Data()));
 387
 388   d->Draw("mean:meanErr:cr","iz==25&&iphi==2","goff");
 389   TGraphErrors *grmean_r=new TGraphErrors(d->GetSelectedRows(),d->GetV3(),d->GetV1(),0,d->GetV2());
 390   grmean_r->SetTitle(";r (cm);#LT#Delta r#varphi#GT");
 391   grmean_r->SetMarkerStyle(20);
 392   grmean_r->SetMarkerSize(1);
 393   grmean_r->Draw("ap");
 394   c->SaveAs(Form("%s/mean_oneZ_r_onePhi.png",outDir.Data()));
 395
 396
 397   d->Draw("meanErr:cphi:cr","iz==25","colz");
 398   c->SaveAs(Form("%s/meanErr_oneZ_phi_allR.png",outDir.Data()));
 399   d->Draw("meanErr:cphi","iz==25&&ir==2");
 400   c->SaveAs(Form("%s/meanErr_oneZ_phi_oneR.png",outDir.Data()));
 401
 402   d->Draw("meanErr:cr:cphi","iz==25","colz");
 403   c->SaveAs(Form("%s/meanErr_oneZ_r_allPhi.png",outDir.Data()));
 404
 405   d->Draw("meanErr:cr","iz==25&&iphi==2");
 406   c->SaveAs(Form("%s/meanErr_oneZ_r_onePhi.png",outDir.Data()));
 407
 408 }