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 "sPhenixStyle.C"
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 smdprofileMaker(const std::string infile = "/sphenix/tg/tg01/coldqcd/jaein213/analysis/ZDCNeutronLocPol/macros/47972/0/zdcneutronlocpol_47972.root", int storenumber = 34782, int runnumber = 47972)
0027 
0028 {
0029   SetsPhenixStyle();
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   TH2D* nxy_hits_south = new TH2D(Form("nxy_hits_all_south"),Form("nxy_hits_all_south"),48,-5.5,5.5,48, -5.92, 5.92);
0051   TH2D* nxy_hits_north = new TH2D(Form("nxy_hits_all_north"),Form("nxy_hits_all_north"),48,-5.5,5.5,48, -5.92, 5.92);
0052  
0053   TH2D* nxy_hits_south_cut = new TH2D(Form("nxy_hits_all_south_cut"),Form("nxy_hits_all_south_cut"),48,-5.5,5.5,48, -5.92, 5.92);
0054   TH2D* nxy_hits_north_cut = new TH2D(Form("nxy_hits_all_north_cut"),Form("nxy_hits_all_north_cut"),48,-5.5,5.5,48, -5.92, 5.92);
0055   TH2D* nxy_hits_south_cut_up = new TH2D(Form("nxy_hits_all_south_cut_up"),Form("nxy_hits_all_south_cut_up"),48,-5.5,5.5,48, -5.92, 5.92);
0056   TH2D* nxy_hits_north_cut_up = new TH2D(Form("nxy_hits_all_north_cut_up"),Form("nxy_hits_all_north_cut_up"),48,-5.5,5.5,48, -5.92, 5.92);
0057   TH2D* nxy_hits_south_up = new TH2D(Form("nxy_hits_all_south_up"),Form("nxy_hits_all_south_up"),48,-5.5,5.5,48, -5.92, 5.92);
0058   TH2D* nxy_hits_north_up = new TH2D(Form("nxy_hits_all_north_up"),Form("nxy_hits_all_north_up"),48,-5.5,5.5,48, -5.92, 5.92);
0059    TH2D* nxy_hits_south_cut_down = new TH2D(Form("nxy_hits_all_south_cut_down"),Form("nxy_hits_all_south_cut_down"),48,-5.5,5.5,48, -5.92, 5.92);
0060   TH2D* nxy_hits_north_cut_down = new TH2D(Form("nxy_hits_all_north_cut_down"),Form("nxy_hits_all_north_cut_down"),48,-5.5,5.5,48, -5.92, 5.92);
0061   TH2D* nxy_hits_south_down = new TH2D(Form("nxy_hits_all_south_down"),Form("nxy_hits_all_south_down"),48,-5.5,5.5,48, -5.92, 5.92);
0062   TH2D* nxy_hits_north_down = new TH2D(Form("nxy_hits_all_north_down"),Form("nxy_hits_all_north_down"),48,-5.5,5.5,48, -5.92, 5.92);
0063   
0064   for (int i = 0; i < 2; i++)
0065   {
0066     gSqrtAsym[i%2] = new TGraphErrors();
0067     gSqrtAsym[i%2]->SetTitle(Form("Square root asymmetry %s North",beam[i].Data()));
0068     gSqrtAsym[i%2+2] = new TGraphErrors();
0069     gSqrtAsym[i%2+2]->SetTitle(Form("Square root asymmetry %s South",beam[i].Data()));
0070 
0071  
0072     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);
0073     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);
0074     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);
0075     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);
0076 
0077     Nup[i] = new TH1D(Form("Nup_%s",beam[i].Data()),Form("Nup_%s",beam[i].Data()),nasymbins,-TMath::Pi(),TMath::Pi()); Nup[i]->Sumw2();
0078     Ndown[i] = new TH1D(Form("Ndown_%s",beam[i].Data()),Form("Ndown_%s",beam[i].Data()),nasymbins,-TMath::Pi(),TMath::Pi()); Ndown[i]->Sumw2();
0079     Sup[i] = new TH1D(Form("Sup_%s",beam[i].Data()),Form("Sup_%s",beam[i].Data()),nasymbins,-TMath::Pi(),TMath::Pi()); Sup[i]->Sumw2();
0080     Sdown[i] = new TH1D(Form("Sdown_%s",beam[i].Data()),Form("Sdown_%s",beam[i].Data()),nasymbins,-TMath::Pi(),TMath::Pi()); Sdown[i]->Sumw2();
0081 
0082   }
0083   //======================================================================//
0084 
0085   //========================== SMD Hit Tree ==============================//
0086   TTree *smdHits = (TTree*)f->Get("smdHits");
0087   int bunchnumber, showerCutN, showerCutS;
0088   float n_x, n_y, s_x, s_y;
0089   float zdcN1_adc, zdcN2_adc, veto_NF, veto_NB;
0090   float zdcS1_adc, zdcS2_adc, veto_SF, veto_SB;
0091   smdHits -> SetBranchAddress ("bunchnumber",       &bunchnumber);
0092   smdHits -> SetBranchAddress ("showerCutN",       &showerCutN);
0093   smdHits -> SetBranchAddress ("showerCutS",       &showerCutS);
0094   smdHits -> SetBranchAddress ("n_x",       &n_x);
0095   smdHits -> SetBranchAddress ("n_y",       &n_y);
0096   smdHits -> SetBranchAddress ("s_x",       &s_x);
0097   smdHits -> SetBranchAddress ("s_y",       &s_y);
0098   smdHits -> SetBranchAddress ("zdcN1_adc",       &zdcN1_adc);
0099   smdHits -> SetBranchAddress ("zdcN2_adc",       &zdcN2_adc);
0100   smdHits -> SetBranchAddress ("zdcS1_adc",       &zdcS1_adc);
0101   smdHits -> SetBranchAddress ("zdcS2_adc",       &zdcS2_adc);
0102   smdHits -> SetBranchAddress ("veto_NF",       &veto_NF);
0103   smdHits -> SetBranchAddress ("veto_NB",       &veto_NB);
0104   smdHits -> SetBranchAddress ("veto_SF",       &veto_SF);
0105   smdHits -> SetBranchAddress ("veto_SB",       &veto_SB);
0106   int nentries = smdHits->GetEntries();
0107   //======================================================================//
0108   
0109   //============================= Option A: Spin DB ================================//
0110   unsigned int qa_level = 0xffff;
0111   SpinDBOutput spin_out("phnxrc");
0112   SpinDBContentv1 spin_cont;
0113   spin_out.StoreDBContent(runnumber,runnumber,qa_level);
0114   SpinDBContent* spin_ptr = &spin_cont;
0115   spin_out.GetDBContentStore(spin_ptr,runnumber);
0116   
0117   int crossingshift = spin_cont.GetCrossingShift();
0118   //crossingshift = 0;
0119   //crossingshift = ixs;
0120   //if (ixs < 0){crossingshift = 120+ixs;}
0121   std::cout << crossingshift << std::endl;
0122   
0123   int bspinpat[120] = {0};
0124   int yspinpat[120] = {0};
0125   for (int i = 0; i < 120; i++)
0126   {
0127     bspinpat[i] = spin_cont.GetSpinPatternBlue(i);
0128     yspinpat[i] = spin_cont.GetSpinPatternYellow(i);
0129     bspinpat[i] *= -1;
0130     yspinpat[i] *= -1;
0131   }
0132   //======================================================================//
0133 
0134   
0135   
0136   //============================= Option B: Get spin pattern manually ================================//
0137   int bpat[120] = {0};
0138   int ypat[120] = {0};
0139   std::string preset_pattern_blue[8];
0140   std::string preset_pattern_yellow[8];
0141   preset_pattern_blue[0] = "+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+*********";
0142   preset_pattern_blue[1] = "-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-*********";
0143   preset_pattern_blue[2] = "+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+*********";
0144   preset_pattern_blue[3] = "-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-*********";
0145   preset_pattern_blue[4] = "++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++-*********";
0146   preset_pattern_blue[5] = "--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--+*********";
0147   preset_pattern_blue[6] = "++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++-*********";
0148   preset_pattern_blue[7] = "--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--+*********";
0149 
0150   
0151   preset_pattern_yellow[0] = "++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++-*********";
0152   preset_pattern_yellow[1] = "++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++-*********";
0153   preset_pattern_yellow[2] = "--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--+*********";
0154   preset_pattern_yellow[3] = "--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--++--+*********";
0155   preset_pattern_yellow[4] = "+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+*********";
0156   preset_pattern_yellow[5] = "+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+--+-++-+-+-+--+-+-+-++-+*********";
0157   preset_pattern_yellow[6] = "-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-*********";
0158   preset_pattern_yellow[7] = "-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-++-+--+-+-+-++-+-+-+--+-*********";
0159 
0160   int spinpatternNo = 3;
0161   if (storenumber == 34485){spinpatternNo = 3;} //pattern names start at '111x111_P1' so index 3  means '111x111_P4' and so on
0162   
0163   for (int i = 0; i < 120; i++)
0164   {
0165     if (preset_pattern_blue[spinpatternNo].at(i) == '+') {bpat[i] = 1;}
0166     else if (preset_pattern_blue[spinpatternNo].at(i) == '-'){bpat[i] = -1;}
0167     else if (preset_pattern_blue[spinpatternNo].at(i) == '*'){bpat[i] = 10;}
0168 
0169     if (preset_pattern_yellow[spinpatternNo].at(i) == '+'){ypat[i] = 1;}
0170     else if (preset_pattern_yellow[spinpatternNo].at(i) == '-'){ypat[i] = -1;}
0171     else if (preset_pattern_yellow[spinpatternNo].at(i) == '*'){ypat[i] = 10;}
0172 
0173   }
0174   //=====================================================================================//
0175 
0176   //======================= process event ==================================//
0177   std::cout << nentries << std::endl;
0178   
0179   int clockOffset = 1; //this is the offset between GL1 and ZDC events
0180 
0181   for (int i = 0; i < nentries - clockOffset; i++)
0182   {
0183     if (i % 1000000 == 0){std::cout << "Entry: " << i << std::endl;}
0184     smdHits->GetEntry(i+clockOffset);
0185    
0186     int sphenix_cross = (bunchnumber + crossingshift) % 120;
0187 
0188     
0189     smdHits->GetEntry(i);
0190 
0191     int bspin = bspinpat[sphenix_cross]; //option A: spinDB
0192     int yspin = yspinpat[sphenix_cross]; //option A: spinDB
0193     //int bspin = bpat[sphenix_cross]; //option B: preset patterns
0194     //int yspin = ypat[sphenix_cross]; //option B: preset patterns
0195 
0196     /*
0197     float ny_offset = 0.85; float nx_offset = 0.275;
0198     float sy_offset = 0.0; float sx_offset = 0.0;
0199 
0200     n_y = n_y - ny_offset;
0201     n_x = n_x - nx_offset;
0202     s_y = s_y - sy_offset;
0203     s_x = s_x - sx_offset;
0204     */
0205 
0206     // phi = -pi/2 to the left of polarized proton going direction (PHENIX convention - https://journals.aps.org/prd/pdf/10.1103/PhysRevD.88.032006)
0207     // SMD goes from -x to +x from left to right when looking at detector (from perspective of Blue-North, Yellow-South)
0208     float n_phi = atan2(n_y, -n_x) - TMath::Pi() / 2;
0209     if (n_phi < -TMath::Pi()) {n_phi += 2 * TMath::Pi();} 
0210     else if (n_phi > TMath::Pi()) {n_phi -= 2 * TMath::Pi();}
0211 
0212     float s_phi = atan2(s_y, -s_x) - TMath::Pi() / 2;
0213     if (s_phi < -TMath::Pi()) {s_phi += 2 * TMath::Pi();} 
0214     else if (s_phi > TMath::Pi()) {s_phi -= 2 * TMath::Pi();}
0215     
0216 
0217     if (zdcN1_adc > ZDC1CUT && zdcN2_adc > ZDC2CUT && veto_NF < VETOCUT && veto_NB < VETOCUT && showerCutN == 1)
0218     {
0219       nxy_hits_north->Fill(n_x, n_y);
0220       if (bspin == 1)
0221       {
0222         nxy_hits_north_up->Fill(n_x, n_y);
0223       }
0224       else if (bspin == -1)
0225       {
0226         nxy_hits_north_down->Fill(n_x, n_y);
0227       }
0228 
0229       if (sqrt(n_x * n_x + n_y * n_y) > RMINCUT && sqrt(n_x * n_x + n_y * n_y) < RMAXCUT)
0230       {
0231         nxy_hits_north_cut->Fill(n_x, n_y);
0232         if (bspin == 1)
0233         {
0234           nxy_hits_north_cut_up->Fill(n_x, n_y);
0235         }
0236         else if (bspin == -1)
0237         {
0238           nxy_hits_north_cut_down->Fill(n_x, n_y);
0239         }
0240 
0241         if (yspin == 1)
0242         {
0243         }
0244         else if (yspin == -1)
0245         {
0246         }
0247       }
0248     }
0249 
0250     if (zdcS1_adc > ZDC1CUT && zdcS2_adc > ZDC2CUT  && veto_SF < VETOCUT && veto_SB < VETOCUT && showerCutS == 1)
0251     {
0252       nxy_hits_south->Fill(s_x, s_y);
0253       if (yspin == 1)
0254       {
0255         nxy_hits_south_up->Fill(s_x, s_y);
0256       }
0257       else if (yspin == -1)
0258       {
0259         nxy_hits_south_down->Fill(s_x, s_y);
0260       }
0261       if (sqrt(s_x * s_x + s_y * s_y) > RMINCUT && sqrt(s_x * s_x + s_y * s_y) < RMAXCUT)
0262       {
0263         nxy_hits_south_cut->Fill(s_x, s_y);
0264         if (yspin == 1)
0265         {
0266           nxy_hits_south_cut_up->Fill(s_x, s_y);
0267         }
0268         else if (yspin == -1)
0269         {
0270           nxy_hits_south_cut_down->Fill(s_x, s_y);
0271         }
0272       }
0273     }
0274   }
0275   //=====================================================================================//
0276 
0277   
0278   
0279 
0280   TFile* ofile = new TFile("Profile.root","RECREATE");
0281   ofile->cd();
0282   nxy_hits_south -> Write();
0283   nxy_hits_south_cut -> Write();
0284   nxy_hits_north -> Write();
0285   nxy_hits_north_cut -> Write();
0286   nxy_hits_south_up -> Write();
0287   nxy_hits_south_cut_up -> Write();
0288   nxy_hits_north_up -> Write();
0289   nxy_hits_north_cut_up -> Write();
0290   nxy_hits_south_down -> Write();
0291   nxy_hits_south_cut_down -> Write();
0292   nxy_hits_north_down -> Write();
0293   nxy_hits_north_cut_down -> Write();
0294 
0295 
0296 
0297   ofile->Close();
0298 }
[ 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!