// Comapre the MC information with the reconstructed
//
+/*
+ after running analysis, read the file, and get component
+ gSystem->Load("libPWG1.so");
+ TFile f("Output.root");
+ AliComparisonDraw * comp = (AliComparisonDraw*)f.Get("AliComparisonDraw");
+ TF1 fl("fl","((min(250./(abs(x+0.000001)),250)-90))",0,2); // length function
+ TF1 fl2("fl2","[0]/((min(250./(abs(x+0.000001)),250)-90))^[1]",0,2);
+ fl2.SetParameter(1,1);
+ fl2.SetParameter(0,1);
+
+
+
+
+
+
+*/
+
+
#include "TSystem.h"
#include "TTimeStamp.h"
#include "TProof.h"
#include "TTree.h"
+#include "TH3F.h"
#include "TH2F.h"
-#include "TH1F.h"
-
+#include "TF1.h"
+#include "TProfile.h"
+#include "TProfile2D.h"
+#include "TGraph2D.h"
+#include "TPad.h"
+#include "TCanvas.h"
+#include "TGraph.h"
//
#include "AliTracker.h"
#include "AliMagF.h"
#include "AliTPCseed.h"
#include "AliTPCclusterMI.h"
//
+#include "AliMathBase.h"
+#include "AliTreeDraw.h"
#include "AliComparisonDraw.h"
ClassImp(AliComparisonDraw)
-
+Bool_t AliComparisonDraw::fBDraw; //option draw temporary results
AliComparisonDraw::AliComparisonDraw():
TObject(),
void AliComparisonDraw::InitHisto(){
//
//
+ // EFFICIENCY
+ //
+ // Efficiency as function of pt
+ fEffTPCPt = new TProfile("Eff_pt","Eff_Pt",50,0.1,3); // physical
+ fEffTPCPtMC = new TProfile("MC_Eff_pt","MC_Eff_Pt",50,0.1,3); // MC - particles make more than 50 rowdigits
+ fEffTPCPtF = new TProfile("F_Eff_pt","F_Eff_Pt",50,0.1,3); // tracking - under condition more than 50 rdigits
+
+ // Efficiency as function of pt
+ fEffTPCTan = new TProfile("Eff_tan","Eff_tan",50,-2.5,2.5); // physical
+ fEffTPCTanMC = new TProfile("MC_Eff_tan","MC_Eff_tan",50,-2.5,2.5); // MC - particles make more than 50 rowdigits
+ fEffTPCTanF = new TProfile("F_Eff_tan","F_Eff_tan",50,-2.5,2.5); // tracking - under condition more than 50 rdigits
+
+ fEffTPCPtTan = new TProfile2D("Eff_pt","Eff_Pt",10,0.1,3,20,-2.,2.);
+ fEffTPCPtTanMC = new TProfile2D("MC_Eff_pt","MC Eff Pt",10,0.1,3,20, -2.,2.);
+ fEffTPCPtTanF = new TProfile2D("MC_Eff_pt","MC Eff Pt",10,0.1,3,20, -2.,2.);
+
+ //
+ // TPC dEdx
+ //
+ fTPCSignalNormTan = new TH2F("CdEdxTan","CdEdxTan",50, -2,2, 40,30,70); // tpc signal normalized to the MC
+ fTPCSignalNormSPhi = new TH2F("CdEdxSPhi","CdEdxSPhi",10,0.0,1,40,30,70); // tpc signal normalized to the MC
+ fTPCSignalNormTPhi = new TH2F("CdEdxTPhi","CdEdxTPhi",10,0.0,2,40,30,70); // tpc signal normalized to the MC
+
+ fTPCSignalNormTanSPhi= new TH3F("CdEdxTanSPhi","CdEdxTanSPhi",20, -2,2, 10,0.0 ,1, 40,30,70); // tpc signal normalized to the mean signal - MC
+ fTPCSignalNormTanTPhi= new TH3F("CdEdxTanTPhi","CdEdxTanTPhi",20, -2,2, 10,0.0 ,1, 40,30,70); // tpc signal normalized to the mean signal - MC
+ fTPCSignalNormTanSPt= new TH3F("CdEdxTanSPt","CdEdxTanSPt",20, -2,2, 10,0.3 ,3, 40,30,70); // tpc signal normalized to the mean signal - MC
+
+
+
+ //
+ // RESOLUTION
+ //
+ fCPtResolTan = new TH2F("C Pt resol","C pt resol",50, -2,2,200,-0.2,0.2);
+ fCPtPullTan = new TH2F("C Pt pull","C pt pull",50, -2,2,200,-5,5);
+ //
+ fCPhiResolTan = new TH2F("CPhiResolTan","CPhiResolTan",50, -2,2,200,-0.025,0.025);
+ // angular resolution - constrained
+ fCTanResolTan = new TH2F("CTanResolTan","CTanResolTan",50, -2,2,200,-0.025,0.025);
+ // angular resolution - constrained
+ fCPtResolTan=new TH2F("CPtResol","CPtResol",50, -2,2,200,-0.2,0.2);;
+ // pt resolution - constrained
+ fCPhiPullTan = new TH2F("CPhiPullTan","CPhiPullTan",50, -2,2,200,-5,5);
+ // angular resolution - constrained
+ fCTanPullTan = new TH2F("CTanPullTan","CTanPullTan",50, -2,2,200,-5,5);
+ // angular resolution - constrained
+ fCPtPullTan=new TH2F("CPtPull","CPtPull",50, -2,2,200,-5,5);
+ // pt resolution - constrained
//
fPtResolLPT = new TH2F("Pt resol","pt resol",10, 0.1,3,200,-0.2,0.2);
fPtResolHPT = new TH2F("Pt resol","pt resol",10, 2,100,200,-0.3,0.3);
+ fPtPullLPT = new TH2F("Pt pool","pt pool",10, 0.1,3,200,-6,6);
+ fPtPullHPT = new TH2F("Pt pool","pt pool",10, 2,100,200,-6,6);
//
- fPtPoolLPT = new TH2F("Pt pool","pt pool",10, 0.1,3,200,-6,6);
- fPtPoolHPT = new TH2F("Pt pool","pt pool",10, 2,100,200,-6,6);
+ fD0TanSPtB1 = new TH3F("DCAyTanSPt","DCAyTanSPt",20,1,2, 10,0.3,2, 100,-4,4);
+ fD1TanSPtB1 = new TH3F("DCAzTanSPt","DCAzTanSPt",20,1,2, 10,0.3,2, 100,-4,4);
+ fD0TanSPtL1 = new TH3F("DCAyTanSPt","DCAyTanSPt",20,0,1, 10,0.3,2, 100,-0.1,0.1);
+ fD1TanSPtL1 = new TH3F("DCAzTanSPt","DCAzTanSPt",20,0,1, 10,0.3,2, 100, -0.1,0.1);
+
+
+
}
-void AliComparisonDraw::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
+void AliComparisonDraw::ProcessEff(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
//
+ // make efficiencies histograms
+ //
+ Float_t kptcut = 0.15;
+ Float_t ktancut=1.;
+ Int_t kmincl =50;
+ Float_t mcpt = infoMC->GetParticle().Pt();
+ Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
+ Bool_t isPrim = infoMC->GetParticle().R()<0.1 && TMath::Abs(infoMC->GetParticle().Vz())<10;
+ //z diamond and
+
+ if (!isPrim) return;
+
+ //pt
+ if (TMath::Abs(tantheta)<ktancut){
+ fEffTPCPt->Fill(mcpt, infoRC->GetStatus(1)==3);
+ fEffTPCPtMC->Fill(mcpt, infoMC->GetRowsWithDigits()>kmincl);
+ if (infoMC->GetRowsWithDigits()>kmincl){
+ fEffTPCPtF->Fill(mcpt, infoRC->GetStatus(1)==3);
+ }
+ }
+
+ //theta
+ if (TMath::Abs(mcpt)>kptcut){
+ fEffTPCTan->Fill(tantheta, infoRC->GetStatus(1)==3);
+ fEffTPCTanMC->Fill(tantheta, infoMC->GetRowsWithDigits()>kmincl);
+ if (infoMC->GetRowsWithDigits()>kmincl){
+ fEffTPCTanF->Fill(tantheta, infoRC->GetStatus(1)==3);
+ }
+ }
//
+ // pt-theta
+ //
+ fEffTPCPtTan->Fill(mcpt,tantheta,infoRC->GetStatus(1)==3);
+ fEffTPCPtTanMC->Fill(mcpt,tantheta,infoMC->GetRowsWithDigits()>50);
+ if (infoMC->GetRowsWithDigits()>kmincl){
+ fEffTPCPtTanF->Fill(mcpt,tantheta,infoRC->GetStatus(1)==3);
+ }
+}
+
+
+void AliComparisonDraw::ProcessResolConstrained(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
+ //
+ //
//
Float_t mcpt = infoMC->GetParticle().Pt();
+ Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
+ Bool_t isPrim = infoMC->GetParticle().R()<0.1 && TMath::Abs(infoMC->GetParticle().Vz())<10;
+ //z diamond and
+
+ if (!isPrim) return;
+ if (infoRC->GetStatus(1)!=3) return;
+ if (!infoRC->GetESDtrack()) return;
+ if (infoRC->GetESDtrack()->GetTPCNcls()<10) return;
+ if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
+
+ //
+ // constrained parameters resolution
+ //
+ const AliExternalTrackParam * cparam = infoRC->GetESDtrack()->GetConstrainedParam();
+ Float_t deltaCPt= (mcpt-cparam->Pt())/mcpt;
+ Float_t pullCPt= (1/mcpt-cparam->OneOverPt())/
+ TMath::Sqrt(cparam->GetSigma1Pt2());
+ Float_t deltaPhi = TMath::ATan2(cparam->Py(),cparam->Px())-
+ TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
+ Float_t pullPhi = deltaPhi/TMath::Sqrt(cparam->GetSigmaSnp2());
+
+ Float_t deltaTan = TMath::ATan2(cparam->Pz(),cparam->Pt())-
+ TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
+ Float_t pullTan = deltaPhi/TMath::Sqrt(cparam->GetSigmaSnp2());
+ fCPtResolTan->Fill(tantheta,deltaCPt);
+ fCPtPullTan->Fill(tantheta,pullCPt);
+ fCPhiResolTan->Fill(tantheta,deltaPhi);
+ fCPhiPullTan->Fill(tantheta,pullPhi);
+ fCTanResolTan->Fill(tantheta,deltaTan);
+ fCTanPullTan->Fill(tantheta,pullTan);
+
+}
+
+
+
+void AliComparisonDraw::ProcessTPCdedx(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
//
//
- if (infoRC->GetStatus(1)==0) return;
- if (!infoRC->GetESDtrack()) return; //buggy line
+ //
+ Float_t mcpt = infoMC->GetParticle().Pt();
+ Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
+ Bool_t isPrim = infoMC->GetParticle().R()<0.1 && TMath::Abs(infoMC->GetParticle().Vz())<10;
+ //z diamond and
+
+ if (!isPrim) return;
+ if (infoRC->GetStatus(1)!=3) return;
+ if (!infoRC->GetESDtrack()) return;
+ if (infoRC->GetESDtrack()->GetTPCNcls()<10) return;
+ if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
+ Float_t mprim = infoMC->GetPrim();
+ if (mprim>1.4) return;
+ if (mprim<0.5) return;
+ if (infoRC->GetESDtrack()->GetTPCsignalN()<50) return;
+ //
+ Float_t ratio = infoRC->GetESDtrack()->GetTPCsignal()/infoMC->GetPrim();
+ Float_t sphi = infoRC->GetESDtrack()->GetInnerParam()->GetSnp();
+ Float_t tphi = sphi/TMath::Sqrt(1-sphi*sphi);
+ if (TMath::Abs(infoMC->GetParticle().GetPdgCode())!=211) return;
+ if (mcpt>0.5){
+ fTPCSignalNormTan->Fill(tantheta,ratio); //only subset
+ }
+ if (TMath::Abs(tantheta)<0.5){
+ fTPCSignalNormSPhi->Fill(sphi,ratio); // only subset
+ fTPCSignalNormTPhi->Fill(tphi,ratio); // only subset
+ }
+ fTPCSignalNormTanSPhi->Fill(tantheta,sphi,ratio);
+ fTPCSignalNormTanTPhi->Fill(tantheta,tphi,ratio);
+ fTPCSignalNormTanSPt->Fill(tantheta,TMath::Sqrt(mcpt),ratio);
+}
+
+void AliComparisonDraw::ProcessDCA(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
+ //
+ //
+ //
+ Float_t mcpt = infoMC->GetParticle().Pt();
+ Float_t tantheta = TMath::Tan(infoMC->GetParticle().Theta()-TMath::Pi()*0.5);
+ Bool_t isPrim = infoMC->GetParticle().R()<0.1 && TMath::Abs(infoMC->GetParticle().Vz())<10;
+ //z diamond and
+ if (!isPrim) return;
+ if (infoRC->GetStatus(1)!=3) return;
+ if (!infoRC->GetESDtrack()) return;
if (infoRC->GetESDtrack()->GetTPCNcls()<10) return;
+ if (!infoRC->GetESDtrack()->GetConstrainedParam()) return;
+ Float_t spt = TMath::Sqrt(mcpt);
+ Float_t dca[2],cov[3];
+ infoRC->GetESDtrack()->GetImpactParameters(dca,cov);
+ Int_t clusterITS[100];
+ if (infoRC->GetESDtrack()->GetITSclusters(clusterITS)==0){
+ fD0TanSPtB1->Fill(tantheta,spt,dca[0]);
+ fD1TanSPtB1->Fill(tantheta,spt,dca[1]);
+ }
+ fD0TanSPtL1->Fill(tantheta,spt,dca[0]);
+ fD1TanSPtL1->Fill(tantheta,spt,dca[1]);
+}
+
+
+
+
+void AliComparisonDraw::Process(AliMCInfo* infoMC, AliESDRecInfo *infoRC){
+ //
+ //
+ //
+ ProcessEff(infoMC,infoRC);
+ ProcessResolConstrained(infoMC,infoRC);
+ ProcessTPCdedx(infoMC, infoRC);
+ ProcessDCA(infoMC, infoRC);
+ Float_t mcpt = infoMC->GetParticle().Pt();
+ Bool_t isPrim = infoMC->GetParticle().R()<0.1 && TMath::Abs(infoMC->GetParticle().Vz())<10;
+ //z diamond and
+
+ if (!isPrim) return;
+ //
+ //
+ if (infoRC->GetStatus(1)==0) return;
+ if (!infoRC->GetESDtrack()) return;
+ if (infoRC->GetESDtrack()->GetTPCNcls()<10) return;
// printf("Pt\t%f\t%f\n",mcpt, infoRC->GetESDtrack()->Pt());
Float_t deltaPt= (mcpt-infoRC->GetESDtrack()->Pt())/mcpt;
Float_t poolPt= (1/mcpt-infoRC->GetESDtrack()->OneOverPt())/
TMath::Sqrt(infoRC->GetESDtrack()->GetSigma1Pt2());
+
fPtResolLPT->Fill(mcpt,deltaPt);
fPtResolHPT->Fill(mcpt,deltaPt);
- fPtPoolLPT->Fill(mcpt,poolPt);
- fPtPoolHPT->Fill(mcpt,poolPt);
-
+ fPtPullLPT->Fill(mcpt,poolPt);
+ fPtPullHPT->Fill(mcpt,poolPt);
+}
+
+
+
+TH1F* AliComparisonDraw::MakeResol(TH2F * his, Int_t integ, Bool_t type){
+ TH1F *hisr, *hism;
+ if (!gPad) new TCanvas;
+ hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);
+ if (type) return hism;
+ else
+ return hisr;
+}
+
+
+TGraph2D * AliComparisonDraw::MakeStat2D(TH3 * his, Int_t delta0, Int_t delta1, Int_t type){
+ //
+ //
+ //
+ // delta - number of bins to integrate
+ // type - 0 - mean value
+
+ TAxis * xaxis = his->GetXaxis();
+ TAxis * yaxis = his->GetYaxis();
+ // TAxis * zaxis = his->GetZaxis();
+ Int_t nbinx = xaxis->GetNbins();
+ Int_t nbiny = yaxis->GetNbins();
+ char name[1000];
+ Int_t icount=0;
+ TGraph2D *graph = new TGraph2D(nbinx*nbiny);
+ TF1 f1("f1","gaus");
+ for (Int_t ix=0; ix<nbinx;ix++)
+ for (Int_t iy=0; iy<nbiny;iy++){
+ Float_t xcenter = xaxis->GetBinCenter(ix);
+ Float_t ycenter = yaxis->GetBinCenter(iy);
+ sprintf(name,"%s_%d_%d",his->GetName(), ix,iy);
+ TH1 *projection = his->ProjectionZ(name,ix-delta0,ix+delta0,iy-delta1,iy+delta1);
+ Float_t stat= 0;
+ if (type==0) stat = projection->GetMean();
+ if (type==1) stat = projection->GetRMS();
+ if (type==2 || type==3){
+ TVectorD vec(3);
+ AliMathBase::LTM((TH1F*)projection,&vec,0.7);
+ if (type==2) stat= vec[1];
+ if (type==3) stat= vec[0];
+ }
+ if (type==4|| type==5){
+ projection->Fit(&f1);
+ if (type==4) stat= f1.GetParameter(1);
+ if (type==5) stat= f1.GetParameter(2);
+ }
+ //printf("%d\t%f\t%f\t%f\n", icount,xcenter, ycenter, stat);
+ graph->SetPoint(icount,xcenter, ycenter, stat);
+ icount++;
+ }
+ return graph;
+}
+
+TGraph * AliComparisonDraw::MakeStat1D(TH3 * his, Int_t delta1, Int_t type){
+ //
+ //
+ //
+ // delta - number of bins to integrate
+ // type - 0 - mean value
+
+ TAxis * xaxis = his->GetXaxis();
+ TAxis * yaxis = his->GetYaxis();
+ // TAxis * zaxis = his->GetZaxis();
+ Int_t nbinx = xaxis->GetNbins();
+ Int_t nbiny = yaxis->GetNbins();
+ char name[1000];
+ Int_t icount=0;
+ TGraph *graph = new TGraph(nbinx);
+ TF1 f1("f1","gaus");
+ for (Int_t ix=0; ix<nbinx;ix++){
+ Float_t xcenter = xaxis->GetBinCenter(ix);
+ // Float_t ycenter = yaxis->GetBinCenter(iy);
+ sprintf(name,"%s_%d",his->GetName(), ix);
+ TH1 *projection = his->ProjectionZ(name,ix-delta1,ix+delta1,0,nbiny);
+ Float_t stat= 0;
+ if (type==0) stat = projection->GetMean();
+ if (type==1) stat = projection->GetRMS();
+ if (type==2 || type==3){
+ TVectorD vec(3);
+ AliMathBase::LTM((TH1F*)projection,&vec,0.7);
+ if (type==2) stat= vec[1];
+ if (type==3) stat= vec[0];
+ }
+ if (type==4|| type==5){
+ projection->Fit(&f1);
+ if (type==4) stat= f1.GetParameter(1);
+ if (type==5) stat= f1.GetParameter(2);
+ }
+ //printf("%d\t%f\t%f\t%f\n", icount,xcenter, ycenter, stat);
+ graph->SetPoint(icount,xcenter, stat);
+ icount++;
+ }
+ return graph;
+}
+
+//
+// Make derived plots
+//
+
+void AliComparisonDraw::MakePlots(){
+ //
+ //
+ //
+ AliComparisonDraw * comp=this;
+
+ TFile *fp = new TFile("picutures.root","recreate");
+ TH1F *hiss=0, *hism=0;
+ TGraph2D * gr=0, gr2=0;
+ TGraph * gr0 = 0;
+ TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
+ //
+ //
+ //
+ hiss = comp->MakeResol(comp->fCPtResolTan,1,0);
+ hiss->SetXTitle("Tan(#theta)");
+ hiss->SetYTitle("#sigmap_{t}/p_{t}");
+ hiss->Draw();
+ hiss->Write("CptResolTan");
+ //
+ //
+ hiss = comp->MakeResol(comp->fCPhiResolTan,1,0);
+ c->cd();
+ hiss->SetXTitle("Tan(#theta)");
+ hiss->SetYTitle("#sigma#phi (rad)");
+ hiss->Draw();
+ fp->cd();
+ hiss->Write("PhiResolTan");
+ //
+ hiss = comp->MakeResol(comp->fCTanResolTan,1,0);
+ c->cd();
+ hiss->SetXTitle("Tan(#theta)");
+ hiss->SetYTitle("#sigma#theta (rad)");
+ hiss->Draw();
+ fp->cd();
+ hiss->Write("ThetaResolTan");
+ //
+ //
+ hiss = comp->MakeResol(comp->fCTanResolTan,1,0);
+ c->cd();
+ hiss->SetXTitle("Tan(#theta)");
+ hiss->SetYTitle("#sigmap_{t}/p_{t} ");
+ hiss->Draw();
+ fp->cd();
+ //
+ //
+ //
+ hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,0);
+ hiss->SetXTitle("Tan(#theta)");
+ hiss->SetYTitle("#sigma_{dEdx}");
+ hiss->Draw();
+ fp->cd();
+ hiss->Write("TPCdEdxResolTan");
+ //
+ //
+ //
+ hiss = comp->MakeResol(comp->fTPCSignalNormTan,4,1);
+ hiss->SetXTitle("Tan(#theta)");
+ hiss->SetYTitle("<dEdx>");
+ hiss->Draw();
+ hiss->Write("TPCdEdxMeanTan");
+ //
+ //
+ gr = comp->MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,4);
+ gr->GetXaxis()->SetTitle("Tan(#theta)");
+ gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
+ gr->GetZaxis()->SetTitle("<dEdx>");
+ gr->Draw("colz");
+ gr->GetHistogram()->Write("TPCdEdxMeanTanPt");
+ //
+ //
+ gr = comp->MakeStat2D(comp->fTPCSignalNormTanSPt,3,1,5);
+ gr->GetXaxis()->SetTitle("Tan(#theta)");
+ gr->GetYaxis()->SetTitle("#sqrt{p_{t}(GeV)}");
+ gr->GetZaxis()->SetTitle("#sigma_{dEdx}");
+ gr->Draw("colz");
+ gr->GetHistogram()->Write("TPCdEdxMeanTanPt");
+ //
+ //
+ //
+ comp->fEffTPCTanF->SetXTitle("Tan(#theta)");
+ comp->fEffTPCTanF->SetYTitle("eff_{findable}");
+ comp->fEffTPCTanF->Draw();
+ comp->fEffTPCTanF->Write("EffTanFindable");
+ //
+ //
+ comp->fEffTPCTan->SetXTitle("Tan(#theta)");
+ comp->fEffTPCTan->SetYTitle("eff_{all}");
+ comp->fEffTPCTan->Draw();
+ comp->fEffTPCTan->Write("EffTanAll");
+ //
+ //DCA resolution
+ //
+ gr0 = comp->MakeStat1D(comp->fD0TanSPtB1,2,5);
+ gr0->GetXaxis()->SetTitle("Tan(#theta)");
+ gr0->GetYaxis()->SetTitle("#sigmaDCA (cm)");
+ gPad->Clear();
+ gr0->Draw("al*");
+ gr->GetHistogram()->Write("DCAResolTan");
+ //
+ //
+ //
+ gr = comp->MakeStat2D(comp->fD0TanSPtB1,4,2,5);
+ gr0->GetXaxis()->SetTitle("Tan(#theta)");
+ gr0->GetYaxis()->SetTitle("#sigmaDCA (cm)");
+ gPad->Clear();
+ gr0->Draw("al*");
+ gr->GetHistogram()->Write("DCAResolSPTTan");
+
+ fp->Close();
+
+
}
+
+