+void BinScan(Int_t npoints){
+
+
+ TTreeSRedirector *pcstream = new TTreeSRedirector("localBins.root","update");
+ TTree * resolScan = (TTree*)pcstream->GetFile()->Get("resolScan");
+ if (!resolScan){
+ TTree * meanNormZ = (TTree*) pcstream->GetFile()->Get("meanNormZ");
+ TTree * meanDistortion = (TTree*) pcstream->GetFile()->Get("meanDistortion");
+ // meanNormZ->SetMarkerStyle(25); meanNormZ->SetMarkerSize(0.5);
+ // meanDistortion->SetMarkerStyle(25); meanDistortion->SetMarkerSize(0.5);
+ // Int_t colors[5]={1,2,3,4,6};
+ // Int_t markers[5]={21,20,23,24,25};
+ // TH2 * his2D;
+ // TObjArray arrFit(3);
+
+ {
+ //Int_t npoints=50000;
+ TCut cutDist="abs(ddeltaR0-vecDMeanRRND.fElements)<1&&abs(ddeltaRPhi0-vecDMeanRPhiRND.fElements)<1";
+ TCut cutGeom="abs(z0)<85&&r0<120.";
+ //
+ Int_t entries1 = meanDistortion->Draw("vecZRND.fElements:vecRRND.fElements:vecPhiRND.fElements","izNorm==1"+cutGeom+cutDist,"goff",npoints);
+ TVectorD vecBR1(entries1,meanDistortion->GetV1());
+ TVectorD vecBZ1(entries1,meanDistortion->GetV2());
+ TVectorD vecBPhi1(entries1,meanDistortion->GetV3());
+ meanDistortion->Draw("ddeltaR0-vecDMeanRRND.fElements:ddeltaRPhi0-vecDMeanRPhiRND.fElements","izNorm==1"+cutGeom+cutDist,"goff",npoints);
+ TVectorD vecDR1(entries1,meanDistortion->GetV1());
+ TVectorD vecDRPhi1(entries1,meanDistortion->GetV2());
+ //
+ Int_t entries0 = meanDistortion->Draw("vecZRND.fElements:vecRRND.fElements:vecPhiRND.fElements","izNorm==0"+cutGeom+cutDist,"goff",npoints);
+ TVectorD vecBR0(entries0,meanDistortion->GetV1());
+ TVectorD vecBZ0(entries0,meanDistortion->GetV2());
+ TVectorD vecBPhi0(entries0,meanDistortion->GetV3());
+ meanDistortion->Draw("ddeltaR0-vecDMeanRRND.fElements:ddeltaRPhi0-vecDMeanRPhiRND.fElements","izNorm==0"+cutGeom+cutDist,"goff",npoints);
+ TVectorD vecDR0(entries0,meanDistortion->GetV1());
+ TVectorD vecDRPhi0(entries0,meanDistortion->GetV2());
+ //
+ TVectorD vecSelR(TMath::Max(entries0,entries1));
+ TVectorD vecSelRPhi(TMath::Max(entries0,entries1));
+ //
+ for (Int_t iz=1; iz<10; iz+=1){
+ for (Int_t ir=1; ir<10; ir+=1){
+ for (Int_t iphi=1; iphi<10; iphi++){
+ Double_t zbin=3*iz;
+ Double_t rbin=3*ir;
+ Double_t phibin=0.025*iphi;
+ Int_t counter=0;
+ //
+ counter=0;
+ for (Int_t ipoint=0; ipoint<entries1; ipoint++){
+ Bool_t isOK=TMath::Abs(vecBZ1[ipoint]/zbin-1)<0.2;
+ isOK&=TMath::Abs(vecBR1[ipoint]/rbin-1.)<0.2;
+ isOK&=TMath::Abs(vecBPhi1[ipoint]/phibin-1.)<0.2;
+ if (isOK) {
+ vecSelRPhi[counter]=vecDRPhi1[ipoint];
+ vecSelR[counter]=vecDR1[ipoint];
+ counter++;
+ }
+ }
+ Double_t meanR1=0,rmsR1=0;
+ Double_t meanRPhi1=0,rmsRPhi1=0;
+ if (counter>3) AliMathBase::EvaluateUni(counter,vecSelR.GetMatrixArray(),meanR1,rmsR1,0.9*counter);
+ if (counter>3) AliMathBase::EvaluateUni(counter,vecSelRPhi.GetMatrixArray(),meanRPhi1,rmsRPhi1,0.9*counter);
+ //
+ counter=0;
+ for (Int_t ipoint=0; ipoint<entries0; ipoint++){
+ Bool_t isOK=TMath::Abs(vecBZ0[ipoint]/zbin-1)<0.2;
+ isOK&=TMath::Abs(vecBR0[ipoint]/rbin-1.)<0.2;
+ isOK&=TMath::Abs(vecBPhi0[ipoint]/phibin-1.)<0.2;
+ if (isOK) {
+ vecSelRPhi[counter]=vecDRPhi0[ipoint];
+ vecSelR[counter]=vecDR0[ipoint];
+ counter++;
+ }
+ }
+ Double_t meanR0=0, rmsR0=0;
+ Double_t meanRPhi0=0, rmsRPhi0=0;
+ if (counter>3) AliMathBase::EvaluateUni(counter,vecSelR.GetMatrixArray(),meanR0,rmsR0,0.9*counter);
+ if (counter>3) AliMathBase::EvaluateUni(counter,vecSelRPhi.GetMatrixArray(),meanRPhi0,rmsRPhi0,0.9*counter);
+ //
+ printf("%f\t%f\t%f\t%f\t%f\n",zbin,rbin,phibin,rmsR0/(rmsR1+0.0001), rmsRPhi0/(rmsRPhi1+0.00001));
+ (*pcstream)<<"resolScan"<<
+ "counter="<<counter<<
+ //
+ "iz="<<iz<<
+ "ir="<<ir<<
+ "iphi="<<iphi<<
+ //
+ "zbin="<<zbin<<
+ "rbin="<<rbin<<
+ "phibin="<<phibin<<
+ //
+ "meanR0="<<meanR0<<
+ "rmsR0="<<rmsR0<<
+ "meanRPhi0="<<meanRPhi0<<
+ "rmsRPhi0="<<rmsRPhi0<<
+ //
+ "meanR1="<<meanR1<<
+ "rmsR1="<<rmsR1<<
+ "meanRPhi1="<<meanRPhi1<<
+ "rmsRPhi1="<<rmsRPhi1<<
+ "\n";
+ }
+ }
+ }
+ }
+ delete pcstream;
+ }
+ //
+ pcstream = new TTreeSRedirector("localBins.root","update");
+ resolScan = (TTree*)pcstream->GetFile()->Get("resolScan");
+ resolScan->SetMarkerStyle(25);
+ //
+ Int_t colors[5]={1,2,3,4,6};
+ Int_t markers[5]={21,20,23,24,25};
+ gStyle->SetTitleFontSize(32);
+ gStyle->SetTitleFontSize(35);
+ //
+ //
+ //
+ for (Int_t itype=0; itype<2; itype++){
+ TCanvas * canvasRes = new TCanvas(TString::Format("canvasRes%d",itype),"canvasRes",800,800);
+ canvasRes->SetRightMargin(0.05);
+ canvasRes->SetLeftMargin(0.2);
+ canvasRes->SetBottomMargin(0.18);
+ canvasRes->Divide(2,3,0,0);
+ TLatex latexDraw;
+ latexDraw.SetTextSize(0.08);
+ //
+ for (Int_t iz=1; iz<6; iz+=2){
+ TLegend * legend0 = new TLegend(0.17,0.3,0.80,0.6,TString::Format("Residuals after mean correction"));
+ TLegend * legend1 = new TLegend(0.07,0.3,0.90,0.6,TString::Format("Residual after applying #it{q(z)} correction"));
+ legend0->SetBorderSize(0);
+ legend1->SetBorderSize(0);
+ for (Int_t ir=1; ir<8; ir+=2){
+ TCut cutR(TString::Format("ir==%d",ir));
+ TCut cutZ(TString::Format("iz==%d",iz));
+ TGraphErrors * gr0=0, *gr1=0;
+ if (itype==0){
+ gr0=TStatToolkit::MakeGraphErrors(resolScan,"10*rmsR0:phibin:10*rmsR0/sqrt(counter)",cutR+cutZ,markers[ir/2],colors[ir/2],0.75);
+ gr1=TStatToolkit::MakeGraphErrors(resolScan,"10*rmsR1:phibin:10*rmsR1/sqrt(counter)",cutR+cutZ,markers[ir/2],colors[ir/2],0.75);
+ gr0->GetYaxis()->SetTitle("#it{#sigma(#Delta_{R}-#bar{#Delta_{R}})} (mm)");
+ }
+ if (itype==1){
+ gr0=TStatToolkit::MakeGraphErrors(resolScan,"10*rmsRPhi0:phibin:10*rmsRPhi0/sqrt(counter)",cutR+cutZ,markers[ir/2],colors[ir/2],0.75);
+ gr1=TStatToolkit::MakeGraphErrors(resolScan,"10*rmsPhi1:phibin:10*rmsPhi1/sqrt(counter)",cutR+cutZ,markers[ir/2],colors[ir/2],0.75);
+ gr0->GetYaxis()->SetTitle("#it{#sigma(#Delta_{R#phi}-#bar{#Delta_{R#phi}})} (mm)");
+ }
+ gr0->GetXaxis()->SetTitle("#Delta#phi bin width");
+ gr1->GetXaxis()->SetTitle("#Delta#phi bin width");
+ SetGraphTDRStyle(gr0);
+ SetGraphTDRStyle(gr1);
+ gr0->GetXaxis()->SetLimits(0,0.25);
+ gr1->GetXaxis()->SetLimits(0,0.25);
+ gr0->SetMinimum(-0.5);
+ gr0->SetMaximum(0.7);
+ gr1->SetMinimum(-0.5);
+ gr1->SetMaximum(0.7);
+ canvasRes->cd(iz)->SetTicks(3,3);
+ canvasRes->cd(iz)->SetGrid(0,3);
+ canvasRes->cd(iz)->SetLeftMargin(0.15);
+ if (ir==1) gr0->Draw("alp");
+ gr0->Draw("lp");
+ canvasRes->cd(iz+1)->SetTicks(3,3);
+ canvasRes->cd(iz+1)->SetGrid(0,3);
+ if (ir==1) gr1->Draw("alp");
+ gr1->Draw("lp");
+ legend0->AddEntry(gr0,TString::Format("#it{#Delta_{R}}=%1.1f (cm)",ir*3.),"p");
+ legend1->AddEntry(gr1,TString::Format("#it{#Delta_{R}}=%1.1f (cm)",ir*3.),"p");
+ //
+ canvasRes->cd(iz);
+ latexDraw.DrawLatex(0.01,-0.45,TString::Format("#Delta_{Z}=%1.0f cm, R<120 cm, |Z|<85 cm ",iz*3.));
+ canvasRes->cd(iz+1);
+ latexDraw.DrawLatex(0.01,-0.45,TString::Format("#Delta_{Z}=%1.0f cm, R<120 cm, |Z|<85 cm ",iz*3.));
+ }
+ if (iz==5){
+ legend0->SetTextSize(0.06);
+ legend1->SetTextSize(0.06);
+ canvasRes->cd(iz);
+ legend0->Draw();
+ canvasRes->cd(iz+1);
+ legend1->Draw();
+ }
+ }
+ if (itype==0){
+ canvasRes->SaveAs("canvasSpaceChargeBinFlucR.pdf");
+ canvasRes->SaveAs("canvasSpaceChargeBinFlucR.png");
+ }
+ if (itype==1){
+ canvasRes->SaveAs("canvasSpaceChargeBinFlucRPhi.pdf");
+ canvasRes->SaveAs("canvasSpaceChargeBinFlucRPhi.png"); //~/hera/alice/miranov/SpaceCharge/Fluctuations/PbPbWithGain/dirmergeAll/dEpsilon20/canvasSpaceChargeBinFlucRPhi.pdf
+ }
+ }