analysis/ZDC-CALIB/zdc_calib.C

0001 //get zdc weigthts per side for 3 modules
0002 //updated for sPHENIX ZDC
0003 //Ejiro Umaka, Peter Steinberg
0004
0005
0006 #define Nmodules 3
0007 float E_n;
0008
0009 float no_booster_weight[3] = {1,1,1};
0010 std::vector<float> skip_zdc_mod = {1,1,1};
0011
0012 //change per distribution
0013 float uncorr_lo = 133.138;
0014 float uncorr_hi = 224.234;
0015
0016
0017 bool calculateWeights(float energy,
0018               std::vector< std::vector<float> >& in_data,
0019               TVectorD& out_weights,
0020               TMatrixD& out_matrix
0021               )
0022 {
0023
0024   int n_zdc = in_data.size();
0025   int dim = n_zdc+1;
0026
0027   std::cout << "Calculate weights for " << n_zdc << " modules!" << std::endl;
0028
0029   TVectorD mom1(dim);
0030
0031   size_t nevt = in_data.at(0).size();
0032   for (size_t i = 0;i<nevt;i++)
0033     {
0034       for (size_t ir=0;ir<n_zdc;ir++)
0035     {
0036       mom1(ir) += in_data.at(ir).at(i)/double(nevt);
0037       for (size_t ic=0;ic<n_zdc;ic++)
0038         {
0039           out_matrix(ir,ic) += 2*in_data.at(ir).at(i)*in_data.at(ic).at(i)/double(nevt);
0040         }
0041     }
0042     }
0043
0044
0045   for (size_t ir=0;ir<n_zdc;ir++)
0046     {
0047       out_matrix(ir,n_zdc)=mom1(ir);
0048       out_matrix(n_zdc,ir)=mom1(ir);
0049     }
0050
0051   TVectorD solution(dim);
0052   solution(n_zdc) = energy;
0053
0054   TDecompLU lu(out_matrix);
0055   lu.Decompose();
0056   bool status = lu.Solve(solution);
0057
0058   if (!status)
0059    {
0060      std::cout << "Fail!" << endl;
0061     }
0062   else
0063     {
0064       std::cout << "Success!" << endl;
0065     }
0066
0067   solution.Print();
0068
0069   out_weights = solution;
0070
0071   return status;
0072 }
0073
0074
0075 //put input file and set constraint energy to 100 GeV
0076 void zdc_calib(TString filename = "", float e = 100)
0077
0078 {
0079
0080    TString s_zdc = "";
0081    for (int iz=0;iz<Nmodules;iz++)
0082     {
0083       s_zdc += TString::Itoa(skip_zdc_mod[iz],10);
0084     }
0085
0086   TFile fout("output"+filename+"_"+s_zdc+"_fout.root","RECREATE");
0087   TH1F h_e_uncorr("h_e_uncorr","Uncorr. energy;E_{uncorr}",3000,0,6000);
0088   TH2F h_e_uncorr2D("h_e_uncorr2D","Uncorr. energy;E_{ZDC} [GeV]};E_{EM} [GeV];",100,0,300,100,0,300);
0089   TH1F h_e_corr("h_e_corr","corr. energy;E_{corr}",3000,0,6000);
0090   TH2F h_e_corr2D("h_e_corr2D","corr. energy;E_{ZDC} [GeV];E_{EM} [GeV]",100,0,600,100,0,600);
0091
0092   int n_zdc = std::accumulate(skip_zdc_mod.begin(),skip_zdc_mod.end(),0);
0093   std::cout << "Working with " << n_zdc << " modules!" << std::endl;
0094
0095   TMatrixD* zdcMatrix = new TMatrixD(n_zdc+1,n_zdc+1);
0096   TVectorD* zdcWeights = new TVectorD(n_zdc+1);
0097
0098
0099   std::vector< std::vector<float> > data;
0100   data.resize(n_zdc);
0101   std::vector< std::vector<float> > all_data;
0102   all_data.resize(n_zdc);
0103
0104   ifstream ifs(filename);
0105
0106   int event;
0107   float x,y;
0108   std::vector<float> z(Nmodules);
0109   int Nevents;
0110
0111
0112
0113 float thissum = 0.;
0114   while (!ifs.eof())
0115     {
0116
0117       ifs >> z[0] >> z[1] >> z[2];
0118
0119       float sum = 0;
0120       float sum_had = 0;
0121       bool fail_event = false;
0122
0123       for (int iz=0;iz<Nmodules;iz++)
0124     {
0125       z[iz] = (z[iz]<4000?z[iz]:0)/no_booster_weight[iz];
0126
0127       if (!skip_zdc_mod[iz] && z[iz]>10) fail_event = true;
0128
0129       sum += z[iz]*skip_zdc_mod[iz];
0130       if (iz>0) sum_had += z[iz]*skip_zdc_mod[iz];
0131     }
0132       if (fail_event) continue;
0133
0134       int ind = 0;
0135       if (sum>uncorr_lo&&sum<uncorr_hi)
0136     {
0137       for (int iz=0;iz<Nmodules;iz++)
0138         {
0139           if (skip_zdc_mod[iz])
0140         {
0141           data.at(ind).push_back(z[iz]);
0142           ind++;
0143         }
0144         }
0145     }
0146
0147       ind = 0;
0148       for (int iz=0;iz<Nmodules;iz++)
0149     {
0150       if (skip_zdc_mod[iz])
0151         {
0152           all_data.at(ind).push_back(z[iz]);
0153           ind++;
0154         }
0155     }
0156
0157       h_e_uncorr.Fill(sum);
0158       h_e_uncorr2D.Fill(sum_had,sum-sum_had);
0159     }
0160
0161   calculateWeights(e, data, *zdcWeights, *zdcMatrix);
0162
0163   int Nevt = all_data.at(0).size();
0164   std::cout << Nevt  << " events used!" << std::endl;
0165   for (int i = 0;i<Nevt;i++)
0166     {
0167
0168       float e_corr = 0;
0169       float e_corr_had = 0;
0170       int ind = 0;
0171       for (int iz=0;iz<Nmodules;iz++)
0172     {
0173       if (skip_zdc_mod[iz])
0174         {
0175           e_corr += all_data.at(ind).at(i) * (*zdcWeights)(ind);
0176           if (iz>0) e_corr_had += all_data.at(ind).at(i) * (*zdcWeights)(ind);
0177           ind++;
0178         }
0179     }
0180
0181       h_e_corr.Fill(e_corr);
0182       h_e_corr2D.Fill(e_corr_had,e_corr-e_corr_had);
0183
0184     }
0185
0186   fout.cd();
0187   zdcWeights->Write("zdcTBWeights");
0188   fout.Write();
0189   fout.Close();
0190
0191 }
0192