[u/mrichter/AliRoot.git] / TPC / Upgrade / macros / AnaDelta.C

#include <TStyle.h>
#include <TROOT.h>
#include <TAxis.h>
#include <TF1.h>
#include <TFile.h>
#include <TH1.h>
#include <THn.h>
#include <TObjArray.h>
#include <TObject.h>
#include <TString.h>
#include <TVectorT.h>
#include <TCanvas.h>
#include <TProfile2D.h>
#include <TGraphErrors.h>
#include <TTreeStream.h>

#include <AliExternalTrackParam.h>
#include <AliTPCComposedCorrection.h>
#include <AliTPCCorrectionLookupTable.h>

#include <AliToyMCEventGenerator.h>

/*

.L $ALICE_ROOT/TPC/Upgrade/macros/AnaDelta.C+g


*/
TVectorD* MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax);
TVectorD* MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax);
TVectorD* MakeArbitraryBinning(const char* bins);

void DumpHn(THn *hn, TTreeSRedirector &stream);

void AnaDelta(TString file, TString outDir=".")
{
  //
  //
  //

  gStyle->SetOptFit();
  
  TTreeSRedirector stream(Form("%s/deltas.root",outDir.Data()));
  gROOT->cd();
  
  TFile f(file);
  THn *hn=(THn*)f.Get("hn");

  DumpHn(hn, stream);

  delete hn;
}


void AnaDeltaTree(TString file, TString outFile="deltas_tree.root")
{
  TFile f(file);
  gROOT->cd();
  TTree *t = (TTree*)f.Get("delta");
  Float_t soneOverPt=0.;
  Float_t radius=0.;
  Float_t trackPhi=0.;
  Float_t trackY=0.;
  Float_t trackZ=0.;
  Float_t resRphi=0.;
  Double_t trackRes=0.;
  Float_t pointY=0.;
  Float_t pointZ=0.;
  Short_t npTRD=0.;
  Short_t event=0.;
  
  t->SetBranchAddress("soneOverPt" , &soneOverPt);
  t->SetBranchAddress("r"          , &radius);
  t->SetBranchAddress("trackPhi"   , &trackPhi);
  t->SetBranchAddress("trackY"     , &trackY);
  t->SetBranchAddress("trackZ"     , &trackZ);
  t->SetBranchAddress("resRphi"    , &resRphi);
  t->SetBranchAddress("trackRes"   , &trackRes);
  t->SetBranchAddress("pointY"     , &pointY);
  t->SetBranchAddress("pointZ"     , &pointZ);
  t->SetBranchAddress("npTRD"      , &npTRD);
  t->SetBranchAddress("event"      , &event);
  
  // make binning
  TVectorD *vR   = MakeLinBinning(10,86.,250.);
  TVectorD *vPhi = MakeLinBinning(18*8,0.,2*TMath::Pi());
  TVectorD *vZ   = MakeLinBinning(50,-250.,250.);

  const Int_t nbins=4;
  Int_t bins[nbins]    = {vR->GetNrows()-1, vPhi->GetNrows()-1, vZ->GetNrows()-1, 80};
//   Int_t bins[nbins]    = {16, 18*5, 50, 80};
  Double_t xmin[nbins] = {86. , 0.,           -250., -2.};
  Double_t xmax[nbins] = {250., 2*TMath::Pi(), 250.,  2.};
  THnF *hn = new THnF("hn", "hn", nbins, bins, xmin, xmax);

  hn->GetAxis(0)->Set(vR  ->GetNrows()-1, vR  ->GetMatrixArray());
  hn->GetAxis(1)->Set(vPhi->GetNrows()-1, vPhi->GetMatrixArray());
  hn->GetAxis(2)->Set(vZ  ->GetNrows()-1, vZ  ->GetMatrixArray());
  
  
  for (Int_t iev=0; iev<t->GetEntries(); ++iev) {
    t->GetEntry(iev);
    
    // cuts
    // -- on trd
    if (npTRD<2) continue;
    Double_t pt=1./TMath::Abs(soneOverPt);
    if (pt<0.8) continue;
    
    Float_t resRphiRandom = resRphi*trackRes;
    Float_t deviation     = trackY+resRphiRandom-pointY;
    
    Double_t xx[4]={radius, trackPhi, trackZ ,deviation};
    hn->Fill(xx);
  }
  
  // do fits and fill tree
  TTreeSRedirector stream(outFile.Data());
  gROOT->cd();

  DumpHn(hn, stream);

  stream.GetFile()->cd();
  hn->Write();
  
  delete hn;
  delete vR;
  delete vPhi;
  delete vZ;
}


void AnaDeltaTree2(TString file/*, TString outDir="."*/)
{
  //
  // NOTE: not finished
  //
  TFile f(file);
  gROOT->cd();
  TTree *t = (TTree*)f.Get("delta");
  Float_t soneOverPt=0.;
  Float_t radius=0.;
  Float_t trackPhi=0.;
  Float_t trackY=0.;
  Float_t trackZ=0.;
  Float_t resRphi=0.;
  Float_t trackRes=0.;
  Float_t pointY=0.;
  Float_t pointZ=0.;
  Float_t npTRD=0.;
  Float_t event=0.;
  
  t->SetBranchAddress("soneOverPt" , &soneOverPt);
  t->SetBranchAddress("r"          , &radius);
  t->SetBranchAddress("trackPhi"   , &trackPhi);
  t->SetBranchAddress("trackY"     , &trackY);
  t->SetBranchAddress("trackZ"     , &trackZ);
  t->SetBranchAddress("resRphi"    , &resRphi);
  t->SetBranchAddress("trackRes"   , &trackRes);
  t->SetBranchAddress("pointY"     , &pointY);
  t->SetBranchAddress("pointZ"     , &pointZ);
  t->SetBranchAddress("npTRD"      , &npTRD);
  t->SetBranchAddress("event"      , &event);

  // make binning
  TVectorD *vZ   = MakeLinBinning(50,-250.,250.);
  TVectorD *vPhi = MakeLinBinning(18*8,0.,TMath::Pi());
  TVectorD *vR   = MakeLinBinning(16,86.,250.);

  TObjArray arrZ(vZ->GetNrows()-1);
  arrZ.SetOwner();
  
  for (Int_t iev=0; iev<t->GetEntries(); ++iev) {
    t->GetEntry(iev);
    
    // cuts
    // -- on trd
    if (npTRD<2) continue;

    Float_t resRphiRandom=resRphi*trackRes;

    Int_t binZ   = TMath::BinarySearch(vZ->GetNrows(),vZ->GetMatrixArray(),(Double_t)trackZ);
    Int_t binPhi = TMath::BinarySearch(vPhi->GetNrows(),vPhi->GetMatrixArray(),(Double_t)trackPhi);
    Int_t binR   = TMath::BinarySearch(vR->GetNrows(),vR->GetMatrixArray(),(Double_t)radius);

    if (binZ<0)   binZ=0;
    if (binPhi<0) binPhi=0;
    if (binR<0)   binR=0;

    TObjArray *arrPhi=(TObjArray*)arrZ.UncheckedAt(binZ);
    if (!arrPhi) {
      arrPhi=new TObjArray(vPhi->GetNrows()-1);
      arrZ.AddAt(arrPhi,binZ);
    }

    TObjArray *arrR=(TObjArray*)arrPhi->UncheckedAt(binPhi);
    if (!arrR) {
      arrR=new TObjArray(vR->GetNrows()-1);
      arrPhi->AddAt(arrR,binPhi);
    }

    TH1S *h = (TH1S*)arrR->UncheckedAt(binR);
    if (!h) {
      h = new TH1S(Form("h_%02d_%02d_%d02",binZ, binPhi, binR),
                   Form("z,phi,r: %02d,%02d,%d02; #Delta r#phi (cm)",binZ, binPhi, binR),
                   80, -2., 2.);
      arrR->AddAt(h, binR);
    }

    h->Fill(trackY+resRphiRandom-pointY);
  }

  // do fits and fill tree
}

void AnaDeltaResiduals(TString fluctuationMap, TString averageMap, TString outFile="deltas_residuals.root")
{
  //
  //
  //

  TFile fFluct(fluctuationMap);
  AliTPCCorrectionLookupTable *corrFluct = (AliTPCCorrectionLookupTable*)fFluct.Get("map");
  fFluct.Close();
  
  TFile fAverage(averageMap);
  AliTPCCorrectionLookupTable *corrAverage = (AliTPCCorrectionLookupTable*)fAverage.Get("map");
  fAverage.Close();
  
  TObjArray *arrMaps = new TObjArray(2);
  arrMaps->Add(corrAverage); // correction with the average Map
  arrMaps->Add(corrFluct);   // distortion with the fluctuation Map
  
  // create the composed correction
  // if the weight are set to +1 and -1, the first map will be responsible for the distortions
  // The second map for the corrections
  // !!!!! In AliTPCComposedCorrection::GetDistortion MakeInverseIterator is called !!!!
  // for this reason we have to add the maps in the wrong order
  
  AliTPCComposedCorrection *residualDistortion = new AliTPCComposedCorrection(arrMaps, AliTPCComposedCorrection::kQueueResidual);
  Float_t dummy=0;
  TVectorD weights(2);
  weights(1)=+1.;
  weights(0)=-AliToyMCEventGenerator::GetSCScalingFactor(corrFluct, corrAverage,dummy);
  residualDistortion->SetWeights(&weights);

  TVectorD *vR   = MakeLinBinning(10,86.,250.);
  TVectorD *vPhi = MakeLinBinning(18*8,0.,2*TMath::Pi());
  TVectorD *vZ   = MakeLinBinning(50,-250.,250.);
  
  const Int_t nbins=4;
  Int_t bins[nbins]    = {vR->GetNrows()-1, vPhi->GetNrows()-1, vZ->GetNrows()-1, 80};
  //   Int_t bins[nbins]    = {16, 18*5, 50, 80};
  Double_t xmin[nbins] = {86. , 0.,           -250., -2.};
  Double_t xmax[nbins] = {250., 2*TMath::Pi(), 250.,  2.};
  THnF *hn = new THnF("hn", "hn", nbins, bins, xmin, xmax);
  
  hn->GetAxis(0)->Set(vR  ->GetNrows()-1, vR  ->GetMatrixArray());
  hn->GetAxis(1)->Set(vPhi->GetNrows()-1, vPhi->GetMatrixArray());
  hn->GetAxis(2)->Set(vZ  ->GetNrows()-1, vZ  ->GetMatrixArray());
  
  AliExternalTrackParam vv;
  
  for (Float_t iz=-245; iz<=245; iz+=1) {
    Short_t roc=(iz>=0)?0:18;
    for (Float_t ir=86; ir<250; ir+=1) {
      for (Float_t iphi=0; iphi<TMath::TwoPi(); iphi+=1*TMath::DegToRad()){
        Float_t x=ir*(Float_t)TMath::Cos(iphi);
        Float_t y=ir*(Float_t)TMath::Sin(iphi);
        Float_t x3[3]    = {x,y,iz};
        Float_t dx3[3]   = {0.,0.,0.};
        residualDistortion->GetDistortion(x3,roc,dx3);

        Double_t ddx3[3]={dx3[0], dx3[1], dx3[2]};
        vv.Global2LocalPosition(ddx3,iphi);

        Double_t xx[4]={ir, iphi, iz ,ddx3[1]};
        hn->Fill(xx);
        
      }
    }
  }

  TTreeSRedirector stream(outFile.Data());
  gROOT->cd();
  
  DumpHn(hn, stream);
  
  stream.GetFile()->cd();
  hn->Write();
  
  delete hn;
  delete vR;
  delete vPhi;
  delete vZ;
  
  delete residualDistortion;
}

void DumpHn(THn *hn, TTreeSRedirector &stream)
{
  TAxis *ar   = hn->GetAxis(0);
  TAxis *aphi = hn->GetAxis(1);
  TAxis *az   = hn->GetAxis(2);
  
  
  for (Int_t iz=0; iz<az->GetNbins(); ++iz) {
    az->SetRange(iz+1,iz+1);
    TObjArray arrFits;
    arrFits.SetName(Form("z_%02d",iz));
    arrFits.SetOwner();
    
    for (Int_t ir=0; ir<ar->GetNbins(); ++ir) {
      ar->SetRange(ir+1,ir+1);
      for (Int_t iphi=0; iphi<aphi->GetNbins(); ++iphi) {
        aphi->SetRange(iphi+1,iphi+1);
        
        TH1 *hProj = hn->Projection(3);
        if (hProj->GetEntries()<1) {
          delete hProj;
          continue;
        }
        
        TF1 fg("fg","gaus",-2,2);
        Double_t cr   = ar->GetBinCenter(ir+1);
        Double_t cphi = aphi->GetBinCenter(iphi+1);
        Double_t cz   = az->GetBinCenter(iz+1);
        hProj->SetNameTitle(Form("h_%02d_%02d_%02d",iz+1, iphi+1, ir+1),
                            Form("z,phi,r: %02d,%02d,%02d (%.2f, %.2f, %.2f)",iz+1,iphi+1,ir+1, cz, cphi, cr )
        );
        hProj->Fit(&fg,"QR");
        arrFits.Add(hProj);

        Double_t mean     = fg.GetParameter(1);
        Double_t meanErr  = fg.GetParError(1);
        Double_t sigma    = fg.GetParameter(2);
        Double_t sigmaErr = fg.GetParError(2);
        Int_t    entries  = hProj->GetEntries();
        Double_t chi2ndf  = fg.GetChisquare()/fg.GetNDF();

        stream << "d" <<
        "ir="        << ir       <<
        "iphi="      << iphi     <<
        "iz="        << iz       <<
        "cr="        << cr       <<
        "cphi="      << cphi     <<
        "cz="        << cz       <<
        "mean="      << mean     <<
        "meanErr="   << meanErr  <<
        "sigma="     << sigma    <<
        "sigmaErr="  << sigmaErr <<
        "entries="   << entries  <<
        "chi2ndf="   << chi2ndf  <<
        "\n";
      }
    }
    stream.GetFile()->cd();
    arrFits.Write(0x0,TObject::kSingleKey);
    gROOT->cd();
  }
  
}

//______________________________________________________________________________
TVectorD* MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
{
  //
  // Make logarithmic binning
  // the user has to delete the array afterwards!!!
  //
  
  //check limits
  if (xmin<1e-20 || xmax<1e-20){
    printf("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
    return MakeLinBinning(nbinsX, xmin, xmax);
  }
  if (xmax<xmin){
    Double_t tmp=xmin;
    xmin=xmax;
    xmax=tmp;
  }
  TVectorD *binLim=new TVectorD(nbinsX+1);
  Double_t first=xmin;
  Double_t last=xmax;
  Double_t expMax=TMath::Log(last/first);
  for (Int_t i=0; i<nbinsX+1; ++i){
    (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
  }
  return binLim;
}

//______________________________________________________________________________
TVectorD* MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
{
  //
  // Make linear binning
  // the user has to delete the array afterwards!!!
  //
  if (xmax<xmin){
    Double_t tmp=xmin;
    xmin=xmax;
    xmax=tmp;
  }
  TVectorD *binLim=new TVectorD(nbinsX+1);
  Double_t first=xmin;
  Double_t last=xmax;
  Double_t binWidth=(last-first)/nbinsX;
  for (Int_t i=0; i<nbinsX+1; ++i){
    (*binLim)[i]=first+binWidth*(Double_t)i;
  }
  return binLim;
}

//_____________________________________________________________________________
TVectorD* MakeArbitraryBinning(const char* bins)
{
  //
  // Make arbitrary binning, bins separated by a ','
  //
  TString limits(bins);
  if (limits.IsNull()){
    printf("Bin Limit string is empty, cannot add the variable");
    return 0x0;
  }
  
  TObjArray *arr=limits.Tokenize(",");
  Int_t nLimits=arr->GetEntries();
  if (nLimits<2){
    printf("Need at leas 2 bin limits, cannot add the variable");
    delete arr;
    return 0x0;
  }
  
  TVectorD *binLimits=new TVectorD(nLimits);
  for (Int_t iLim=0; iLim<nLimits; ++iLim){
    (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
  }
  
  delete arr;
  return binLimits;
}

void PlotFromTree(TTree *d, TString outDir=".")
{
  TCanvas *c=new TCanvas;
  gStyle->SetOptStat(0);
  d->SetMarkerStyle(20);
  d->SetMarkerSize(1);
  
  TProfile2D pRZ("pRZ",";z (cm); r(cm)",50,-250,250,10,85,250);
  d->Draw("entries:cr:cz>>pRZ","","profcolz");
  pRZ.GetZaxis()->UnZoom();
  c->SaveAs(Form("%s/entries_average.png",outDir.Data()));
  d->Draw("entries:cr:cz>>pRZ","iphi==2","profcolz");
  c->SaveAs(Form("%s/entries_onePhi.png",outDir.Data()));
  
  pRZ.SetMaximum(0.04);
  d->Draw("meanErr:cr:cz>>pRZ","","profcolz");
  c->SaveAs(Form("%s/meanErr_average.png",outDir.Data()));
  d->Draw("meanErr:cr:cz>>pRZ","iphi==2","profcolz");
  c->SaveAs(Form("%s/meanErr_onePhi.png",outDir.Data()));
  
  
  d->Draw("mean:cphi:cr","iz==25","colz");
  c->SaveAs(Form("%s/mean_oneZ_phi_allR.png",outDir.Data()));
  d->Draw("mean:meanErr:cphi","iz==25&&ir==2","goff");
  TGraphErrors *grmean_phi=new TGraphErrors(d->GetSelectedRows(),d->GetV3(),d->GetV1(),0,d->GetV2());
  grmean_phi->SetTitle(";#varphi;#LT#Delta r#varphi#GT");
  grmean_phi->SetMarkerStyle(20);
  grmean_phi->SetMarkerSize(1);
  grmean_phi->Draw("ap");
  c->SaveAs(Form("%s/mean_oneZ_phi_oneR.png",outDir.Data()));
  
  d->Draw("mean:cr:cphi","iz==25","colz");
  c->SaveAs(Form("%s/mean_oneZ_r_allPhi.png",outDir.Data()));
  
  d->Draw("mean:meanErr:cr","iz==25&&iphi==2","goff");
  TGraphErrors *grmean_r=new TGraphErrors(d->GetSelectedRows(),d->GetV3(),d->GetV1(),0,d->GetV2());
  grmean_r->SetTitle(";r (cm);#LT#Delta r#varphi#GT");
  grmean_r->SetMarkerStyle(20);
  grmean_r->SetMarkerSize(1);
  grmean_r->Draw("ap");
  c->SaveAs(Form("%s/mean_oneZ_r_onePhi.png",outDir.Data()));
  
  
  d->Draw("meanErr:cphi:cr","iz==25","colz");
  c->SaveAs(Form("%s/meanErr_oneZ_phi_allR.png",outDir.Data()));
  d->Draw("meanErr:cphi","iz==25&&ir==2");
  c->SaveAs(Form("%s/meanErr_oneZ_phi_oneR.png",outDir.Data()));
  
  d->Draw("meanErr:cr:cphi","iz==25","colz");
  c->SaveAs(Form("%s/meanErr_oneZ_r_allPhi.png",outDir.Data()));
  
  d->Draw("meanErr:cr","iz==25&&iphi==2");
  c->SaveAs(Form("%s/meanErr_oneZ_r_onePhi.png",outDir.Data()));
  
}
Commit	Line	Data
1a18cb2b	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>
03de6136	11	#include <TVectorT.h>
	12	#include <TCanvas.h>
	13	#include <TProfile2D.h>
9bdb0975	14	#include <TGraphErrors.h>
1a18cb2b	15	#include <TTreeStream.h>
1a18cb2b	16
6a37efb5	17	#include <AliExternalTrackParam.h>
	18	#include <AliTPCComposedCorrection.h>
	19	#include <AliTPCCorrectionLookupTable.h>
	20
	21	#include <AliToyMCEventGenerator.h>
	22
1a18cb2b	23	/*
	24
	25	.L $ALICE_ROOT/TPC/Upgrade/macros/AnaDelta.C+g
	26
	27
	28	*/
03de6136	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);
1a18cb2b	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
03de6136	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(188,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
9bdb0975	133	void AnaDeltaTree2(TString file/, TString outDir="."/)
03de6136	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(188,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
6a37efb5	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);
6a37efb5	231	arrMaps->Add(corrAverage); // correction with the average Map
22814dc7	232	arrMaps->Add(corrFluct); // distortion with the fluctuation Map
6a37efb5	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
22814dc7	237	// !!!!! In AliTPCComposedCorrection::GetDistortion MakeInverseIterator is called !!!!
	238	// for this reason we have to add the maps in the wrong order
	239
6a37efb5	240	AliTPCComposedCorrection *residualDistortion = new AliTPCComposedCorrection(arrMaps, AliTPCComposedCorrection::kQueueResidual);
	241	Float_t dummy=0;
	242	TVectorD weights(2);
22814dc7	243	weights(1)=+1.;
22814dc7	244	weights(0)=-AliToyMCEventGenerator::GetSCScalingFactor(corrFluct, corrAverage,dummy);
6a37efb5	245	residualDistortion->SetWeights(&weights);
	246
	247	TVectorD *vR = MakeLinBinning(10,86.,250.);
	248	TVectorD vPhi = MakeLinBinning(188,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	}
03de6136	299
	300	void DumpHn(THn *hn, TTreeSRedirector &stream)
	301	{
1a18cb2b	302	TAxis *ar = hn->GetAxis(0);
	303	TAxis *aphi = hn->GetAxis(1);
	304	TAxis *az = hn->GetAxis(2);
03de6136	305
1a18cb2b	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);
03de6136	317
1a18cb2b	318	TH1 *hProj = hn->Projection(3);
	319	if (hProj->GetEntries()<1) {
	320	delete hProj;
	321	continue;
	322	}
03de6136	323
1a18cb2b	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);
ee3eac14	328	hProj->SetNameTitle(Form("h_%02d_%02d_%02d",iz+1, iphi+1, ir+1),
ee3eac14	329	Form("z,phi,r: %02d,%02d,%02d (%.2f, %.2f, %.2f)",iz+1,iphi+1,ir+1, cz, cphi, cr )
03de6136	330	);
1a18cb2b	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();
03de6136	340
1a18cb2b	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	}
03de6136	361
03de6136	362	}
1a18cb2b	363
03de6136	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]=firstTMath::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;
1a18cb2b	441	}
1a18cb2b	442
03de6136	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");
9bdb0975	467	TGraphErrors *grmean_phi=new TGraphErrors(d->GetSelectedRows(),d->GetV3(),d->GetV1(),0,d->GetV2());
03de6136	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");
9bdb0975	478	TGraphErrors *grmean_r=new TGraphErrors(d->GetSelectedRows(),d->GetV3(),d->GetV1(),0,d->GetV2());
03de6136	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
6a37efb5	497	}
	498
	499