Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:15:28 

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