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