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(corrAverage); // correction with the average Map
 232   arrMaps->Add(corrFluct);   // distortion with the fluctuation 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   // !!!!! In AliTPCComposedCorrection::GetDistortion MakeInverseIterator is called !!!!
 238   // for this reason we have to add the maps in the wrong order
 239
 240   AliTPCComposedCorrection *residualDistortion = new AliTPCComposedCorrection(arrMaps, AliTPCComposedCorrection::kQueueResidual);
 241   Float_t dummy=0;
 242   TVectorD weights(2);
 243   weights(1)=+1.;
 244   weights(0)=-AliToyMCEventGenerator::GetSCScalingFactor(corrFluct, corrAverage,dummy);
 245   residualDistortion->SetWeights(&weights);
 246
 247   TVectorD *vR   = MakeLinBinning(10,86.,250.);
 248   TVectorD *vPhi = MakeLinBinning(18*8,0.,2*TMath::Pi());
 249   TVectorD *vZ   = MakeLinBinning(50,-250.,250.);
 250
 251   const Int_t nbins=4;
 252   Int_t bins[nbins]    = {vR->GetNrows()-1, vPhi->GetNrows()-1, vZ->GetNrows()-1, 80};
 253   //   Int_t bins[nbins]    = {16, 18*5, 50, 80};
 254   Double_t xmin[nbins] = {86. , 0.,           -250., -2.};
 255   Double_t xmax[nbins] = {250., 2*TMath::Pi(), 250.,  2.};
 256   THnF *hn = new THnF("hn", "hn", nbins, bins, xmin, xmax);
 257
 258   hn->GetAxis(0)->Set(vR  ->GetNrows()-1, vR  ->GetMatrixArray());
 259   hn->GetAxis(1)->Set(vPhi->GetNrows()-1, vPhi->GetMatrixArray());
 260   hn->GetAxis(2)->Set(vZ  ->GetNrows()-1, vZ  ->GetMatrixArray());
 261
 262   AliExternalTrackParam vv;
 263
 264   for (Float_t iz=-245; iz<=245; iz+=1) {
 265     Short_t roc=(iz>=0)?0:18;
 266     for (Float_t ir=86; ir<250; ir+=1) {
 267       for (Float_t iphi=0; iphi<TMath::TwoPi(); iphi+=1*TMath::DegToRad()){
 268         Float_t x=ir*(Float_t)TMath::Cos(iphi);
 269         Float_t y=ir*(Float_t)TMath::Sin(iphi);
 270         Float_t x3[3]    = {x,y,iz};
 271         Float_t dx3[3]   = {0.,0.,0.};
 272         residualDistortion->GetDistortion(x3,roc,dx3);
 273
 274         Double_t ddx3[3]={dx3[0], dx3[1], dx3[2]};
 275         vv.Global2LocalPosition(ddx3,iphi);
 276
 277         Double_t xx[4]={ir, iphi, iz ,ddx3[1]};
 278         hn->Fill(xx);
 279
 280       }
 281     }
 282   }
 283
 284   TTreeSRedirector stream(outFile.Data());
 285   gROOT->cd();
 286
 287   DumpHn(hn, stream);
 288
 289   stream.GetFile()->cd();
 290   hn->Write();
 291
 292   delete hn;
 293   delete vR;
 294   delete vPhi;
 295   delete vZ;
 296
 297   delete residualDistortion;
 298 }
 299
 300 void DumpHn(THn *hn, TTreeSRedirector &stream)
 301 {
 302   TAxis *ar   = hn->GetAxis(0);
 303   TAxis *aphi = hn->GetAxis(1);
 304   TAxis *az   = hn->GetAxis(2);
 305
 306
 307   for (Int_t iz=0; iz<az->GetNbins(); ++iz) {
 308     az->SetRange(iz+1,iz+1);
 309     TObjArray arrFits;
 310     arrFits.SetName(Form("z_%02d",iz));
 311     arrFits.SetOwner();
 312
 313     for (Int_t ir=0; ir<ar->GetNbins(); ++ir) {
 314       ar->SetRange(ir+1,ir+1);
 315       for (Int_t iphi=0; iphi<aphi->GetNbins(); ++iphi) {
 316         aphi->SetRange(iphi+1,iphi+1);
 317
 318         TH1 *hProj = hn->Projection(3);
 319         if (hProj->GetEntries()<1) {
 320           delete hProj;
 321           continue;
 322         }
 323
 324         TF1 fg("fg","gaus",-2,2);
 325         Double_t cr   = ar->GetBinCenter(ir+1);
 326         Double_t cphi = aphi->GetBinCenter(iphi+1);
 327         Double_t cz   = az->GetBinCenter(iz+1);
 328         hProj->SetNameTitle(Form("h_%02d_%02d_%02d",iz+1, iphi+1, ir+1),
 329                             Form("z,phi,r: %02d,%02d,%02d (%.2f, %.2f, %.2f)",iz+1,iphi+1,ir+1, cz, cphi, cr )
 330         );
 331         hProj->Fit(&fg,"QR");
 332         arrFits.Add(hProj);
 333
 334         Double_t mean     = fg.GetParameter(1);
 335         Double_t meanErr  = fg.GetParError(1);
 336         Double_t sigma    = fg.GetParameter(2);
 337         Double_t sigmaErr = fg.GetParError(2);
 338         Int_t    entries  = hProj->GetEntries();
 339         Double_t chi2ndf  = fg.GetChisquare()/fg.GetNDF();
 340
 341         stream << "d" <<
 342         "ir="        << ir       <<
 343         "iphi="      << iphi     <<
 344         "iz="        << iz       <<
 345         "cr="        << cr       <<
 346         "cphi="      << cphi     <<
 347         "cz="        << cz       <<
 348         "mean="      << mean     <<
 349         "meanErr="   << meanErr  <<
 350         "sigma="     << sigma    <<
 351         "sigmaErr="  << sigmaErr <<
 352         "entries="   << entries  <<
 353         "chi2ndf="   << chi2ndf  <<
 354         "\n";
 355       }
 356     }
 357     stream.GetFile()->cd();
 358     arrFits.Write(0x0,TObject::kSingleKey);
 359     gROOT->cd();
 360   }
 361
 362 }
 363
 364 //______________________________________________________________________________
 365 TVectorD* MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
 366 {
 367   //
 368   // Make logarithmic binning
 369   // the user has to delete the array afterwards!!!
 370   //
 371
 372   //check limits
 373   if (xmin<1e-20 || xmax<1e-20){
 374     printf("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
 375     return MakeLinBinning(nbinsX, xmin, xmax);
 376   }
 377   if (xmax<xmin){
 378     Double_t tmp=xmin;
 379     xmin=xmax;
 380     xmax=tmp;
 381   }
 382   TVectorD *binLim=new TVectorD(nbinsX+1);
 383   Double_t first=xmin;
 384   Double_t last=xmax;
 385   Double_t expMax=TMath::Log(last/first);
 386   for (Int_t i=0; i<nbinsX+1; ++i){
 387     (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
 388   }
 389   return binLim;
 390 }
 391
 392 //______________________________________________________________________________
 393 TVectorD* MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
 394 {
 395   //
 396   // Make linear binning
 397   // the user has to delete the array afterwards!!!
 398   //
 399   if (xmax<xmin){
 400     Double_t tmp=xmin;
 401     xmin=xmax;
 402     xmax=tmp;
 403   }
 404   TVectorD *binLim=new TVectorD(nbinsX+1);
 405   Double_t first=xmin;
 406   Double_t last=xmax;
 407   Double_t binWidth=(last-first)/nbinsX;
 408   for (Int_t i=0; i<nbinsX+1; ++i){
 409     (*binLim)[i]=first+binWidth*(Double_t)i;
 410   }
 411   return binLim;
 412 }
 413
 414 //_____________________________________________________________________________
 415 TVectorD* MakeArbitraryBinning(const char* bins)
 416 {
 417   //
 418   // Make arbitrary binning, bins separated by a ','
 419   //
 420   TString limits(bins);
 421   if (limits.IsNull()){
 422     printf("Bin Limit string is empty, cannot add the variable");
 423     return 0x0;
 424   }
 425
 426   TObjArray *arr=limits.Tokenize(",");
 427   Int_t nLimits=arr->GetEntries();
 428   if (nLimits<2){
 429     printf("Need at leas 2 bin limits, cannot add the variable");
 430     delete arr;
 431     return 0x0;
 432   }
 433
 434   TVectorD *binLimits=new TVectorD(nLimits);
 435   for (Int_t iLim=0; iLim<nLimits; ++iLim){
 436     (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
 437   }
 438
 439   delete arr;
 440   return binLimits;
 441 }
 442
 443 void PlotFromTree(TTree *d, TString outDir=".")
 444 {
 445   TCanvas *c=new TCanvas;
 446   gStyle->SetOptStat(0);
 447   d->SetMarkerStyle(20);
 448   d->SetMarkerSize(1);
 449
 450   TProfile2D pRZ("pRZ",";z (cm); r(cm)",50,-250,250,10,85,250);
 451   d->Draw("entries:cr:cz>>pRZ","","profcolz");
 452   pRZ.GetZaxis()->UnZoom();
 453   c->SaveAs(Form("%s/entries_average.png",outDir.Data()));
 454   d->Draw("entries:cr:cz>>pRZ","iphi==2","profcolz");
 455   c->SaveAs(Form("%s/entries_onePhi.png",outDir.Data()));
 456
 457   pRZ.SetMaximum(0.04);
 458   d->Draw("meanErr:cr:cz>>pRZ","","profcolz");
 459   c->SaveAs(Form("%s/meanErr_average.png",outDir.Data()));
 460   d->Draw("meanErr:cr:cz>>pRZ","iphi==2","profcolz");
 461   c->SaveAs(Form("%s/meanErr_onePhi.png",outDir.Data()));
 462
 463
 464   d->Draw("mean:cphi:cr","iz==25","colz");
 465   c->SaveAs(Form("%s/mean_oneZ_phi_allR.png",outDir.Data()));
 466   d->Draw("mean:meanErr:cphi","iz==25&&ir==2","goff");
 467   TGraphErrors *grmean_phi=new TGraphErrors(d->GetSelectedRows(),d->GetV3(),d->GetV1(),0,d->GetV2());
 468   grmean_phi->SetTitle(";#varphi;#LT#Delta r#varphi#GT");
 469   grmean_phi->SetMarkerStyle(20);
 470   grmean_phi->SetMarkerSize(1);
 471   grmean_phi->Draw("ap");
 472   c->SaveAs(Form("%s/mean_oneZ_phi_oneR.png",outDir.Data()));
 473
 474   d->Draw("mean:cr:cphi","iz==25","colz");
 475   c->SaveAs(Form("%s/mean_oneZ_r_allPhi.png",outDir.Data()));
 476
 477   d->Draw("mean:meanErr:cr","iz==25&&iphi==2","goff");
 478   TGraphErrors *grmean_r=new TGraphErrors(d->GetSelectedRows(),d->GetV3(),d->GetV1(),0,d->GetV2());
 479   grmean_r->SetTitle(";r (cm);#LT#Delta r#varphi#GT");
 480   grmean_r->SetMarkerStyle(20);
 481   grmean_r->SetMarkerSize(1);
 482   grmean_r->Draw("ap");
 483   c->SaveAs(Form("%s/mean_oneZ_r_onePhi.png",outDir.Data()));
 484
 485
 486   d->Draw("meanErr:cphi:cr","iz==25","colz");
 487   c->SaveAs(Form("%s/meanErr_oneZ_phi_allR.png",outDir.Data()));
 488   d->Draw("meanErr:cphi","iz==25&&ir==2");
 489   c->SaveAs(Form("%s/meanErr_oneZ_phi_oneR.png",outDir.Data()));
 490
 491   d->Draw("meanErr:cr:cphi","iz==25","colz");
 492   c->SaveAs(Form("%s/meanErr_oneZ_r_allPhi.png",outDir.Data()));
 493
 494   d->Draw("meanErr:cr","iz==25&&iphi==2");
 495   c->SaveAs(Form("%s/meanErr_oneZ_r_onePhi.png",outDir.Data()));
 496
 497 }
 498
 499