Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-04-04 08:12:08 

0001 #include <cmath>
0002 
0003 #include <uspin/SpinDBOutput.h>
0004 #include <uspin/SpinDBInput.h>
0005 #include <uspin/SpinDBContentv1.h>
0006 #include <TFile.h>
0007 
0008 float sqrtasym(float N1, float N2, float N3, float N4);
0009 float sqrtasymerr(float N1, float N2, float N3, float N4, float E1, float E2, float E3, float E4);
0010 void FitAsym(TGraphErrors *gAsym, double &eps, double &epserr, double &chi2);
0011 void DrawAsym(int runnumber, TGraphErrors *gSqrtAsymBlue, TGraphErrors *gSqrtAsymYellow);
0012 void DrawCanvas(int runnumber, TH2 *hits_up, TH2 *hits_down, TGraphErrors *gRawSimpleAsym, TGraphErrors *gRawSqrtAsym);
0013 
0014 
0015 R__LOAD_LIBRARY(libuspin.so)
0016 
0017 float ZDC1CUT = 100; // keep above (nominal 65)
0018 float ZDC2CUT = 15; // keep above (nominal 25)
0019 float VETOCUT = 150; // keep below (nominal 150)
0020 float RMINCUT = 2; // keep above (nominal 2)
0021 float RMAXCUT = 4; // keep below (nominal 4)
0022 
0023 
0024 //void drawAsym(const std::string infile = "fakeAsymmetry/FakeAsymmetry.root", int storenumber = 34485, int runnumber = 42797)
0025 //void drawAsym(const std::string infile = "store34485/42797/smdmerge.root", int storenumber = 34485, int runnumber = 42797)
0026 void drawAsym(const std::string infile = "47972/0/zdcneutronlocpol_47972_all.root", int storenumber = 34782, int runnumber = 47972)
0027 //void drawAsym(const std::string infile = "store34492/42836/smdmerge_old.root", int storenumber = 34492, int runnumber = 42836)
0028 {
0029 
0030   TFile *f = new TFile(infile.c_str());
0031 
0032   //========================== Hists/Graphs ==============================//
0033   TString beam[2] = {"Blue", "Yellow"};
0034 
0035   //== 1: Blue North, 2: Yellow North, 3: Blue South, 4: Yellow South
0036   TGraphErrors *gSqrtAsym[4];
0037   
0038   TH2D *nxy_hits_up[2];
0039   TH2D *nxy_hits_down[2];
0040   TH2D *sxy_hits_up[2];
0041   TH2D *sxy_hits_down[2];
0042 
0043   TH1D *Nup[2];
0044   TH1D *Ndown[2];
0045   TH1D *Sup[2];
0046   TH1D *Sdown[2];
0047   
0048   int nasymbins = 16; // defined over entire [-pi,pi]: For sqrt asym there are nasymbins/2 total bins so this must be divisible by 2
0049 
0050   for (int i = 0; i < 2; i++)
0051   {
0052     gSqrtAsym[i%2] = new TGraphErrors();
0053     gSqrtAsym[i%2]->SetTitle(Form("Square root asymmetry %s North",beam[i].Data()));
0054     gSqrtAsym[i%2+2] = new TGraphErrors();
0055     gSqrtAsym[i%2+2]->SetTitle(Form("Square root asymmetry %s South",beam[i].Data()));
0056 
0057     nxy_hits_up[i] = new TH2D(Form("nxy_hits_up_%s",beam[i].Data()),Form("nxy_hits_up_%s",beam[i].Data()),48,-5.5,5.5,48, -5.92, 5.92);
0058     nxy_hits_down[i] = new TH2D(Form("nxy_hits_down_%s",beam[i].Data()),Form("nxy_hits_down_%s",beam[i].Data()),48,-5.5,5.5,48, -5.92, 5.92);
0059     sxy_hits_up[i] = new TH2D(Form("sxy_hits_up_%s",beam[i].Data()),Form("sxy_hits_up_%s",beam[i].Data()),48,-5.5,5.5,48, -5.92, 5.92);
0060     sxy_hits_down[i] = new TH2D(Form("sxy_hits_down_%s",beam[i].Data()),Form("sxy_hits_down_%s",beam[i].Data()),48,-5.5,5.5,48, -5.92, 5.92);
0061 
0062     Nup[i] = new TH1D(Form("Nup_%s",beam[i].Data()),Form("Nup_%s",beam[i].Data()),nasymbins,-TMath::Pi(),TMath::Pi()); Nup[i]->Sumw2();
0063     Ndown[i] = new TH1D(Form("Ndown_%s",beam[i].Data()),Form("Ndown_%s",beam[i].Data()),nasymbins,-TMath::Pi(),TMath::Pi()); Ndown[i]->Sumw2();
0064     Sup[i] = new TH1D(Form("Sup_%s",beam[i].Data()),Form("Sup_%s",beam[i].Data()),nasymbins,-TMath::Pi(),TMath::Pi()); Sup[i]->Sumw2();
0065     Sdown[i] = new TH1D(Form("Sdown_%s",beam[i].Data()),Form("Sdown_%s",beam[i].Data()),nasymbins,-TMath::Pi(),TMath::Pi()); Sdown[i]->Sumw2();
0066 
0067   }
0068   //======================================================================//
0069 
0070   //========================== SMD Hit Tree ==============================//
0071   TTree *smdHits = (TTree*)f->Get("smdHits");
0072   int bunchnumber, showerCutN, showerCutS;
0073   float n_x, n_y, s_x, s_y;
0074   float zdcN1_adc, zdcN2_adc, veto_NF, veto_NB;
0075   float zdcS1_adc, zdcS2_adc, veto_SF, veto_SB;
0076   smdHits -> SetBranchAddress ("bunchnumber",       &bunchnumber);
0077   smdHits -> SetBranchAddress ("showerCutN",       &showerCutN);
0078   smdHits -> SetBranchAddress ("showerCutS",       &showerCutS);
0079   smdHits -> SetBranchAddress ("n_x",       &n_x);
0080   smdHits -> SetBranchAddress ("n_y",       &n_y);
0081   smdHits -> SetBranchAddress ("s_x",       &s_x);
0082   smdHits -> SetBranchAddress ("s_y",       &s_y);
0083   smdHits -> SetBranchAddress ("zdcN1_adc",       &zdcN1_adc);
0084   smdHits -> SetBranchAddress ("zdcN2_adc",       &zdcN2_adc);
0085   smdHits -> SetBranchAddress ("zdcS1_adc",       &zdcS1_adc);
0086   smdHits -> SetBranchAddress ("zdcS2_adc",       &zdcS2_adc);
0087   smdHits -> SetBranchAddress ("veto_NF",       &veto_NF);
0088   smdHits -> SetBranchAddress ("veto_NB",       &veto_NB);
0089   smdHits -> SetBranchAddress ("veto_SF",       &veto_SF);
0090   smdHits -> SetBranchAddress ("veto_SB",       &veto_SB);
0091   int nentries = smdHits->GetEntries();
0092   //======================================================================//
0093   
0094   //============================= Option A: Spin DB ================================//
0095   unsigned int qa_level = 0xffff;
0096   SpinDBOutput spin_out("phnxrc");
0097   SpinDBContentv1 spin_cont;
0098   spin_out.StoreDBContent(runnumber,runnumber,qa_level);
0099   SpinDBContent* spin_ptr = &spin_cont;
0100   spin_out.GetDBContentStore(spin_ptr,runnumber);
0101   
0102   int crossingshift = spin_cont.GetCrossingShift();
0103   //crossingshift = 0;
0104   //crossingshift = ixs;
0105   //if (ixs < 0){crossingshift = 120+ixs;}
0106   std::cout << crossingshift << std::endl;
0107   
0108   int bspinpat[120] = {0};
0109   int yspinpat[120] = {0};
0110   for (int i = 0; i < 120; i++)
0111   {
0112     bspinpat[i] = spin_cont.GetSpinPatternBlue(i);
0113     yspinpat[i] = spin_cont.GetSpinPatternYellow(i);
0114 
0115   }
0116   //======================================================================//
0117 
0118   
0119   
0120   //============================= Option B: Get spin pattern manually ================================//
0121   int bpat[120] = {0};
0122   int ypat[120] = {0};
0123   std::string preset_pattern_blue[8];
0124   std::string preset_pattern_yellow[8];
0125   preset_pattern_blue[0] = "+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+*********";
0126   preset_pattern_blue[1] = "-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-*********";
0127   preset_pattern_blue[2] = "+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+*********";
0128   preset_pattern_blue[3] = "-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-*********";
0129   preset_pattern_blue[4] = "++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++-*********";
0130   preset_pattern_blue[5] = "--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--+*********";
0131   preset_pattern_blue[6] = "++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++-*********";
0132   preset_pattern_blue[7] = "--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--+*********";
0133 
0134   
0135   preset_pattern_yellow[0] = "++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++-*********";
0136   preset_pattern_yellow[1] = "++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++-*********";
0137   preset_pattern_yellow[2] = "--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--+*********";
0138   preset_pattern_yellow[3] = "--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--+*********";
0139   preset_pattern_yellow[4] = "+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+*********";
0140   preset_pattern_yellow[5] = "+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+*********";
0141   preset_pattern_yellow[6] = "-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-*********";
0142   preset_pattern_yellow[7] = "-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-*********";
0143 
0144   int spinpatternNo = 3;
0145   if (storenumber == 34485){spinpatternNo = 3;} //pattern names start at '111x111_P1' so index 3  means '111x111_P4' and so on
0146   
0147   for (int i = 0; i < 120; i++)
0148   {
0149     if (preset_pattern_blue[spinpatternNo].at(i) == '+') {bpat[i] = 1;}
0150     else if (preset_pattern_blue[spinpatternNo].at(i) == '-'){bpat[i] = -1;}
0151     else if (preset_pattern_blue[spinpatternNo].at(i) == '*'){bpat[i] = 10;}
0152 
0153     if (preset_pattern_yellow[spinpatternNo].at(i) == '+'){ypat[i] = 1;}
0154     else if (preset_pattern_yellow[spinpatternNo].at(i) == '-'){ypat[i] = -1;}
0155     else if (preset_pattern_yellow[spinpatternNo].at(i) == '*'){ypat[i] = 10;}
0156 
0157   }
0158   //=====================================================================================//
0159 
0160   //======================= process event ==================================//
0161   std::cout << nentries << std::endl;
0162   
0163   int clockOffset = 0; //this is the offset between GL1 and ZDC events
0164 
0165   for (int i = 0; i < nentries - clockOffset; i++)
0166   {
0167     if (i % 1000000 == 0){std::cout << "Entry: " << i << std::endl;}
0168     smdHits->GetEntry(i+clockOffset);
0169    
0170     int sphenix_cross = (bunchnumber + crossingshift) % 120;
0171 
0172     
0173     smdHits->GetEntry(i);
0174 
0175     int bspin = bspinpat[sphenix_cross]; //option A: spinDB
0176     int yspin = yspinpat[sphenix_cross]; //option A: spinDB
0177     //int bspin = bpat[sphenix_cross]; //option B: preset patterns
0178     //int yspin = ypat[sphenix_cross]; //option B: preset patterns
0179 
0180     /*
0181     float ny_offset = 0.85; float nx_offset = 0.275;
0182     float sy_offset = 0.0; float sx_offset = 0.0;
0183 
0184     n_y = n_y - ny_offset;
0185     n_x = n_x - nx_offset;
0186     s_y = s_y - sy_offset;
0187     s_x = s_x - sx_offset;
0188     */
0189 
0190     // phi = -pi/2 to the left of polarized proton going direction (PHENIX convention - https://journals.aps.org/prd/pdf/10.1103/PhysRevD.88.032006)
0191     // SMD goes from -x to +x from left to right when looking at detector (from perspective of Blue-North, Yellow-South)
0192     float n_phi = atan2(n_y, -n_x) - TMath::Pi() / 2;
0193     if (n_phi < -TMath::Pi()) {n_phi += 2 * TMath::Pi();} 
0194     else if (n_phi > TMath::Pi()) {n_phi -= 2 * TMath::Pi();}
0195 
0196     float s_phi = atan2(s_y, -s_x) - TMath::Pi() / 2;
0197     if (s_phi < -TMath::Pi()) {s_phi += 2 * TMath::Pi();} 
0198     else if (s_phi > TMath::Pi()) {s_phi -= 2 * TMath::Pi();}
0199     
0200 
0201     if (zdcN1_adc > ZDC1CUT && zdcN2_adc > ZDC2CUT && veto_NF < VETOCUT && veto_NB < VETOCUT && showerCutN == 1)
0202     {
0203       if (sqrt(n_x*n_x + n_y*n_y) > RMINCUT && sqrt(n_x*n_x+n_y*n_y) < RMAXCUT)
0204       {
0205     if (bspin == 1)
0206     {
0207       nxy_hits_up[0]->Fill(n_x,n_y);
0208       Nup[0]->Fill(n_phi);
0209     }
0210     else if (bspin == -1)
0211     {
0212       nxy_hits_down[0]->Fill(n_x,n_y);
0213       Ndown[0]->Fill(n_phi);
0214     }
0215         
0216     if (yspin == 1)
0217     {
0218       nxy_hits_up[1]->Fill(n_x,n_y);
0219       Nup[1]->Fill(n_phi);
0220     }
0221     else if (yspin == -1)
0222     {
0223       nxy_hits_down[1]->Fill(n_x,n_y);
0224       Ndown[1]->Fill(n_phi);
0225     }
0226 
0227       }
0228     }
0229 
0230 
0231     if (zdcS1_adc > ZDC1CUT && zdcS2_adc > ZDC2CUT  && veto_SF < VETOCUT && veto_SB < VETOCUT && showerCutS == 1)
0232     {
0233       if (sqrt(s_x*s_x + s_y*s_y) > RMINCUT && sqrt(s_x*s_x+s_y*s_y) < RMAXCUT)
0234       {
0235         if (bspin == 1)
0236     {
0237       sxy_hits_up[0]->Fill(s_x,s_y);
0238       Sup[0]->Fill(s_phi);
0239     }
0240     else if (bspin == -1)
0241     {
0242       sxy_hits_down[0]->Fill(s_x,s_y);
0243       Sdown[0]->Fill(s_phi);
0244     }
0245 
0246     if (yspin == 1)
0247     {
0248       sxy_hits_up[1]->Fill(s_x,s_y);
0249       Sup[1]->Fill(s_phi);
0250     }
0251     else if (yspin == -1)
0252     {
0253       sxy_hits_down[1]->Fill(s_x,s_y);
0254       Sdown[1]->Fill(s_phi);
0255     }
0256       }
0257     }
0258   }
0259   //=====================================================================================//
0260 
0261   
0262   
0263 
0264   // get sqrt asymmetry
0265   for (int i = 0; i < nasymbins; i++)
0266   {
0267     float phi = i*(2*TMath::Pi()/nasymbins) - (TMath::Pi() - TMath::Pi()/nasymbins);
0268 
0269     int phibin = i; // N_Left
0270     int phibin2 = (phibin + (int)(nasymbins/2.)) % nasymbins; // N_Right
0271     
0272 
0273     //== 1: Blue North, 2: Yellow North, 3: Blue South, 4: Yellow South
0274     float sqrtrawasym[4] = {0};
0275     float sqrtrawasymerr[4] = {0};
0276 
0277     for (int j = 0; j < 2; j++)
0278     {
0279 
0280       sqrtrawasym[j%2] = sqrtasym(Nup[j%2]->GetBinContent(phibin+1),Nup[j%2]->GetBinContent(phibin2+1),
0281                   Ndown[j%2]->GetBinContent(phibin+1),Ndown[j%2]->GetBinContent(phibin2+1));
0282 
0283       sqrtrawasymerr[j%2] = sqrtasymerr(Nup[j%2]->GetBinContent(phibin+1),Nup[j%2]->GetBinContent(phibin2+1),
0284                     Ndown[j%2]->GetBinContent(phibin+1),Ndown[j%2]->GetBinContent(phibin2+1),
0285                     Nup[j%2]->GetBinError(phibin+1),Nup[j%2]->GetBinError(phibin2+1),
0286                     Ndown[j%2]->GetBinError(phibin+1),Ndown[j%2]->GetBinError(phibin2+1));
0287 
0288       sqrtrawasym[j%2+2] = sqrtasym(Sup[j%2]->GetBinContent(phibin+1),Sup[j%2]->GetBinContent(phibin2+1),
0289                     Sdown[j%2]->GetBinContent(phibin+1),Sdown[j%2]->GetBinContent(phibin2+1));
0290 
0291       sqrtrawasymerr[j%2+2] = sqrtasymerr(Sup[j%2]->GetBinContent(phibin+1),Sup[j%2]->GetBinContent(phibin2+1),
0292                       Sdown[j%2]->GetBinContent(phibin+1),Sdown[j%2]->GetBinContent(phibin2+1),
0293                       Sup[j%2]->GetBinError(phibin+1),Sup[j%2]->GetBinError(phibin2+1),
0294                       Sdown[j%2]->GetBinError(phibin+1),Sdown[j%2]->GetBinError(phibin2+1));    
0295       
0296       gSqrtAsym[j%2]->SetPoint(i,phi,sqrtrawasym[j%2]);
0297       gSqrtAsym[j%2]->SetPointError(i,0,sqrtrawasymerr[j%2]);
0298 
0299       gSqrtAsym[j%2+2]->SetPoint(i,phi,sqrtrawasym[j%2+2]);
0300       gSqrtAsym[j%2+2]->SetPointError(i,0,sqrtrawasymerr[j%2+2]);
0301     }
0302 
0303   }
0304 
0305 
0306   
0307   TCanvas *ncall = new TCanvas("ncall","ncall",1200,600);
0308   //DrawCanvas(runnumber,nxy_hits_up[0],nxy_hits_down[0],gSqrtAsym[0],gSqrtAsym[1]);
0309   DrawAsym(runnumber,gSqrtAsym[0],gSqrtAsym[1]);
0310 
0311   TCanvas *scall = new TCanvas("scall","scall",1200,600);
0312   //DrawCanvas(runnumber,sxy_hits_up[1],sxy_hits_down[1],gSqrtAsym[2],gSqrtAsym[3]);
0313   DrawAsym(runnumber,gSqrtAsym[2],gSqrtAsym[3]);
0314 
0315   std::string outname = Form("asym_run%d.root", runnumber);
0316   TFile outfile(outname.c_str(), "RECREATE");
0317   for (int i = 0; i < 4; ++i) { gSqrtAsym[i]->Write(); }
0318   for (int i = 0; i < 2; ++i)
0319   {
0320     nxy_hits_up[i]->Write(); nxy_hits_down[i]->Write();
0321     sxy_hits_up[i]->Write(); sxy_hits_down[i]->Write();
0322     Nup[i]->Write(); Ndown[i]->Write(); Sup[i]->Write(); Sdown[i]->Write();
0323   }
0324   ncall->Write(); scall->Write();
0325   outfile.Close();
0326 
0327 }
0328 
0329 
0330 
0331 float sqrtasym(float N1, float N2, float N3, float N4){
0332   float asym = (sqrt(N1*N4)-sqrt(N3*N2))/(sqrt(N1*N4)+sqrt(N3*N2));
0333   return asym;
0334 }
0335 
0336 
0337 
0338 float sqrtasymerr(float N1, float N2, float N3, float N4, float E1, float E2, float E3, float E4){
0339   float D1 = sqrt(N2*N3*N4/N1)*E1;
0340   float D2 = sqrt(N1*N3*N4/N2)*E2;
0341   float D3 = sqrt(N1*N2*N4/N3)*E3;
0342   float D4 = sqrt(N1*N2*N3/N4)*E4;
0343   float asym_err = (1/pow(sqrt(N1*N4)+sqrt(N3*N2),2))*sqrt(pow(D1,2)+pow(D2,2)+pow(D3,2)+pow(D4,2));
0344   return asym_err;
0345 }
0346 
0347 
0348 void FitAsym(TGraphErrors *gAsym, double &eps, double &epserr, double &chi2)
0349 {
0350   double fitLow = -TMath::Pi()/2;
0351   double fitHigh = TMath::Pi()/2;
0352 
0353   TF1 *sin1 = new TF1("sin1","-[0] * TMath::Sin(x - [1])",fitLow,fitHigh);
0354   sin1->SetParLimits(0,0,0.2);
0355   sin1->SetParLimits(1,0,0.05);
0356 
0357   gAsym->GetYaxis()->SetRangeUser(-0.05,0.05);
0358   gAsym->GetXaxis()->SetRangeUser(fitLow,fitHigh);
0359   gAsym->Fit("sin1","MR");
0360   eps = sin1->GetParameter(0);
0361   epserr = sin1->GetParError(0);
0362   chi2 = sin1->GetChisquare();
0363 }
0364 
0365 
0366 void DrawAsym(int runnumber, TGraphErrors *gSqrtAsymBlue, TGraphErrors *gSqrtAsymYellow)
0367 {
0368   TPad *pad1 = new TPad("pad1", "Pad1", 0.0, 0.0, 0.5, 1.0);
0369   pad1->SetRightMargin(0.1);
0370   pad1->SetTopMargin(0.1); 
0371   pad1->SetLeftMargin(0.15); 
0372   pad1->SetBottomMargin(0.15); 
0373   TPad *pad2 = new TPad("pad2", "Pad2", 0.5, 0.0, 1.0, 1.0);
0374   pad2->SetRightMargin(0.1);
0375   pad2->SetTopMargin(0.1); 
0376   pad2->SetLeftMargin(0.15); 
0377   pad2->SetBottomMargin(0.15);
0378 
0379   pad1->Draw();
0380   pad2->Draw();
0381 
0382   // only fit in a length pi interval
0383   double fitLow = -TMath::Pi()/2;
0384   double fitHigh = TMath::Pi()/2;
0385 
0386 
0387   pad1->cd();
0388 
0389   TF1 *sin1 = new TF1("sin1","-[0] * TMath::Sin(x - [1])",fitLow,fitHigh);
0390   sin1->SetParLimits(0,0,0.2);
0391   //sin1->SetParLimits(1,0,0.2);
0392 
0393   gSqrtAsymBlue->GetXaxis()->SetTitle("#phi, rad");
0394   gSqrtAsymBlue->GetXaxis()->SetTitleSize(0.045);
0395   gSqrtAsymBlue->GetXaxis()->SetTitleOffset(0.85);
0396   gSqrtAsymBlue->GetYaxis()->SetTitle("#epsilon(#phi)");
0397   gSqrtAsymBlue->GetYaxis()->SetTitleSize(0.045);
0398   gSqrtAsymBlue->GetYaxis()->SetTitleOffset(1.35);
0399   gSqrtAsymBlue->GetYaxis()->SetRangeUser(-0.05,0.05);
0400   gSqrtAsymBlue->GetXaxis()->SetRangeUser(fitLow,fitHigh);
0401   gSqrtAsymBlue->SetMarkerStyle(21);
0402   gSqrtAsymBlue->Fit("sin1","MR");
0403   gSqrtAsymBlue->Draw("ap");
0404   
0405   TLine *line1 = new TLine(fitLow, 0, fitHigh, 0);
0406   line1->SetLineStyle(2);
0407   line1->SetLineWidth(2);
0408   line1->SetLineColor(kBlack);
0409   line1->Draw();
0410 
0411   TLatex *latex = new TLatex();
0412   latex->SetTextSize(0.04);
0413   latex->SetTextFont(42);
0414   latex->DrawLatexNDC(0.6, 0.8, "#bf{sPHENIX} internal");
0415   latex->DrawLatexNDC(0.6, 0.75, Form("Run %d", runnumber));
0416 
0417   TLatex *latexFit = new TLatex();
0418   latexFit->SetTextSize(0.03);
0419   latexFit->SetTextAlign(13);
0420 
0421   
0422   latexFit->DrawLatexNDC(0.2, 0.875, "|#epsilon_{raw}|sin(#phi - #phi_{0}) fit:");
0423   latexFit->DrawLatexNDC(0.2, 0.825, Form("|#epsilon_{raw}| = %.2f%% #pm %.2f%%", sin1->GetParameter(0)*100, sin1->GetParError(0)*100));
0424   latexFit->DrawLatexNDC(0.2, 0.775, Form("#phi_{0} + #pi = %.2f #pm %.2f", sin1->GetParameter(1), sin1->GetParError(1)));
0425   latexFit->DrawLatexNDC(0.2, 0.715, Form("#chi^{2}/NDF = %.2f/%d", sin1->GetChisquare(), sin1->GetNDF()));
0426 
0427 
0428   pad2->cd();
0429   TF1 *sin2 = new TF1("sin2","-[0] * TMath::Sin(x - [1])",fitLow,fitHigh);
0430   sin2->SetParLimits(0,0,0.2);
0431   //sin2->SetParLimits(1,0,0.2);
0432 
0433   gSqrtAsymYellow->GetXaxis()->SetTitle("#phi, rad");
0434   gSqrtAsymYellow->GetXaxis()->SetTitleSize(0.045);
0435   gSqrtAsymYellow->GetXaxis()->SetTitleOffset(0.85);
0436   gSqrtAsymYellow->GetYaxis()->SetTitle("#epsilon(#phi)");
0437   gSqrtAsymYellow->GetYaxis()->SetTitleSize(0.045);
0438   gSqrtAsymYellow->GetYaxis()->SetTitleOffset(1.35);
0439   gSqrtAsymYellow->GetYaxis()->SetRangeUser(-0.05,0.05);
0440   gSqrtAsymYellow->GetXaxis()->SetRangeUser(fitLow,fitHigh);
0441   gSqrtAsymYellow->SetMarkerStyle(21);
0442   gSqrtAsymYellow->Fit("sin2","MR");
0443   gSqrtAsymYellow->Draw("ap");
0444   
0445   
0446   TLine *line2 = new TLine(fitLow, 0, fitHigh, 0);
0447   line2->SetLineStyle(2);
0448   line2->SetLineWidth(2);
0449   line2->SetLineColor(kBlack);
0450   line2->Draw();
0451 
0452 
0453   TLatex *latexFit2 = new TLatex();
0454   latexFit2->SetTextSize(0.03);
0455   latexFit2->SetTextAlign(13);
0456 
0457   
0458   latexFit2->DrawLatexNDC(0.2, 0.875, "|#epsilon_{raw}|sin(#phi - #phi_{0}) fit:");
0459   latexFit2->DrawLatexNDC(0.2, 0.825, Form("|#epsilon_{raw}| = %.2f%% #pm %.2f%%", sin2->GetParameter(0)*100, sin2->GetParError(0)*100));
0460   latexFit2->DrawLatexNDC(0.2, 0.775, Form("#phi_{0} + #pi = %.2f #pm %.2f", sin2->GetParameter(1), sin2->GetParError(1)));
0461   latexFit2->DrawLatexNDC(0.2, 0.715, Form("#chi^{2}/NDF = %.2f/%d", sin2->GetChisquare(), sin2->GetNDF()));
0462 
0463 
0464 }
0465 
0466 
0467 void DrawCanvas(int runnumber, TH2 *hits_up, TH2 *hits_down, TGraphErrors *gSqrtAsymBlue, TGraphErrors *gSqrtAsymYellow)
0468 {
0469 
0470 // Create custom pads
0471   TPad *pad1 = new TPad("pad1", "Pad1", 0.0, 0.5, 0.5, 1.0);
0472   TPad *pad2 = new TPad("pad2", "Pad2", 0.5, 0.5, 1.0, 1.0);
0473   TPad *pad3 = new TPad("pad3", "Pad3", 0.0, 0.0, 0.5, 0.5); 
0474   TPad *pad4 = new TPad("pad4", "Pad4", 0.5, 0.0, 1.0, 0.5);
0475 
0476   // Draw pads on the canvas
0477   pad1->Draw();
0478   pad2->Draw();
0479   pad3->Draw();
0480   pad4->Draw();
0481 
0482   
0483 
0484   pad1->cd();
0485   hits_up->SetTitle("p^{#uparrow} + p");
0486   hits_up->GetXaxis()->SetTitle("x (cm)");
0487   hits_up->GetYaxis()->SetTitle("y (cm)");
0488   hits_up->Draw("colz");
0489   pad2->cd();
0490   hits_down->SetTitle("p^{#downarrow} + p");
0491   hits_down->Draw("colz");
0492   hits_down->GetXaxis()->SetTitle("x (cm)");
0493   hits_down->GetYaxis()->SetTitle("y (cm)");
0494   
0495 
0496   // only fit in a length pi interval
0497   double fitLow = -TMath::Pi()/2;
0498   double fitHigh = TMath::Pi()/2;
0499 
0500   pad3->cd();
0501   TF1 *sin1 = new TF1("sin1","-[0] * TMath::Sin(x - [1])",fitLow,fitHigh);
0502   sin1->SetParLimits(0,0,0.1);
0503 
0504   gSqrtAsymBlue->GetXaxis()->SetTitle("#phi, rad");
0505   gSqrtAsymBlue->GetXaxis()->SetTitleSize(0.045);
0506   gSqrtAsymBlue->GetXaxis()->SetTitleOffset(0.85);
0507   gSqrtAsymBlue->GetYaxis()->SetTitle("#epsilon(#phi)");
0508   gSqrtAsymBlue->GetYaxis()->SetTitleSize(0.045);
0509   gSqrtAsymBlue->GetYaxis()->SetTitleOffset(0.85);
0510   gSqrtAsymBlue->GetYaxis()->SetRangeUser(-0.05,0.05);
0511   gSqrtAsymBlue->GetXaxis()->SetRangeUser(fitLow,fitHigh);
0512   gSqrtAsymBlue->SetMarkerStyle(21);
0513   gSqrtAsymBlue->Fit("sin1","MR");
0514   gSqrtAsymBlue->Draw("ap");
0515   
0516   TLine *line = new TLine(fitLow, 0, fitHigh, 0); 
0517   line->SetLineStyle(2);
0518   line->SetLineWidth(2);
0519   line->SetLineColor(kBlack);
0520   line->Draw();
0521   line->Draw();
0522 
0523   TLatex *latex = new TLatex();
0524   latex->SetTextSize(0.04);
0525   latex->SetTextFont(42);
0526   latex->DrawLatexNDC(0.65, 0.8, "#bf{sPHENIX} internal");
0527   latex->DrawLatexNDC(0.65, 0.75, Form("Run %d", runnumber));
0528 
0529   TLatex *latexFit = new TLatex();
0530   latexFit->SetTextSize(0.04);
0531   latexFit->SetTextAlign(13); 
0532 
0533   
0534   latexFit->DrawLatexNDC(0.2, 0.85, "|#epsilon_{raw}|sin(#phi - #phi_{0}) fit:");
0535   latexFit->DrawLatexNDC(0.2, 0.8, Form("|#epsilon_{raw}| = %.2f%% #pm %.2f%%", sin1->GetParameter(0)*100, sin1->GetParError(0)*100));
0536   latexFit->DrawLatexNDC(0.2, 0.75, Form("#phi_{0} + #pi = %.2f #pm %.2f", sin1->GetParameter(1), sin1->GetParError(1)));
0537   latexFit->DrawLatexNDC(0.2, 0.7, Form("#chi^{2}/NDF = %.2f/%d", sin1->GetChisquare(), sin1->GetNDF()));
0538 
0539 
0540   pad4->cd();
0541 
0542   TF1 *sin2 = new TF1("sin2","-[0] * TMath::Sin(x - [1])",fitLow,fitHigh);
0543   sin2->SetParLimits(0,0,0.1);
0544 
0545 
0546   
0547   gSqrtAsymYellow->GetXaxis()->SetTitle("#phi, rad");
0548   gSqrtAsymYellow->GetXaxis()->SetTitleSize(0.045);
0549   gSqrtAsymYellow->GetXaxis()->SetTitleOffset(0.85);
0550   gSqrtAsymYellow->GetYaxis()->SetTitle("#epsilon(#phi)");
0551   gSqrtAsymYellow->GetYaxis()->SetTitleSize(0.045);
0552   gSqrtAsymYellow->GetYaxis()->SetTitleOffset(0.85);
0553   gSqrtAsymYellow->GetYaxis()->SetRangeUser(-0.05,0.05);
0554   gSqrtAsymYellow->GetXaxis()->SetRangeUser(fitLow,fitHigh);
0555   gSqrtAsymYellow->SetMarkerStyle(21);
0556   gSqrtAsymYellow->Fit("sin2","MR");
0557   gSqrtAsymYellow->Draw("ap");
0558   
0559   
0560   line->SetLineStyle(2);
0561   line->SetLineWidth(2);
0562   line->SetLineColor(kBlack);
0563   line->Draw();
0564   line->Draw();
0565 
0566 
0567   TLatex *latexFit2 = new TLatex();
0568   latexFit2->SetTextSize(0.04);
0569   latexFit2->SetTextAlign(13);  
0570 
0571   // Draw each line
0572   latexFit2->DrawLatexNDC(0.2, 0.85, "|#epsilon_{raw}|sin(#phi - #phi_{0}) fit:");
0573   latexFit2->DrawLatexNDC(0.2, 0.8, Form("|#epsilon_{raw}| = %.2f%% #pm %.2f%%", sin2->GetParameter(0)*100, sin2->GetParError(0)*100));
0574   latexFit2->DrawLatexNDC(0.2, 0.75, Form("#phi_{0} + #pi = %.2f #pm %.2f", sin2->GetParameter(1), sin2->GetParError(1)));
0575   latexFit2->DrawLatexNDC(0.2, 0.7, Form("#chi^{2}/NDF = %.2f/%d", sin2->GetChisquare(), sin2->GetNDF()));
0576 
0577 
0578 
0579 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!