LuminosityCounterGoodRuns/run/analyze.C

0001 static const float radius_EM = 93.5;
0002 static const float minz_EM = -130.23;
0003 static const float maxz_EM = 130.23;
0004
0005 static const float radius_IH = 127.503;
0006 static const float minz_IH = -170.299;
0007 static const float maxz_IH = 170.299;
0008
0009 static const float radius_OH = 225.87;
0010 static const float minz_OH = -301.683;
0011 static const float maxz_OH = 301.683;
0012
0013 float get_emcal_mineta_zcorrected(float zvertex) {
0014   float z = minz_EM - zvertex;
0015   float eta_zcorrected = asinh(z / (float)radius_EM);
0016   return eta_zcorrected;
0017 }
0018
0019 float get_emcal_maxeta_zcorrected(float zvertex) {
0020   float z = maxz_EM - zvertex;
0021   float eta_zcorrected = asinh(z / (float)radius_EM);
0022   return eta_zcorrected;
0023 }
0024
0025 float get_ihcal_mineta_zcorrected(float zvertex) {
0026   float z = minz_IH - zvertex;
0027   float eta_zcorrected = asinh(z / (float)radius_IH);
0028   return eta_zcorrected;
0029 }
0030
0031 float get_ihcal_maxeta_zcorrected(float zvertex) {
0032   float z = maxz_IH - zvertex;
0033   float eta_zcorrected = asinh(z / (float)radius_IH);
0034   return eta_zcorrected;
0035 }
0036
0037 float get_ohcal_mineta_zcorrected(float zvertex) {
0038   float z = minz_OH - zvertex;
0039   float eta_zcorrected = asinh(z / (float)radius_OH);
0040   return eta_zcorrected;
0041 }
0042
0043 float get_ohcal_maxeta_zcorrected(float zvertex) {
0044   float z = maxz_OH - zvertex;
0045   float eta_zcorrected = asinh(z / (float)radius_OH);
0046   return eta_zcorrected;
0047 }
0048
0049 bool check_bad_jet_eta(float jet_eta, float zertex, float jet_radius) {
0050   float emcal_mineta = get_emcal_mineta_zcorrected(zertex);
0051   float emcal_maxeta = get_emcal_maxeta_zcorrected(zertex);
0052   float ihcal_mineta = get_ihcal_mineta_zcorrected(zertex);
0053   float ihcal_maxeta = get_ihcal_maxeta_zcorrected(zertex);
0054   float ohcal_mineta = get_ohcal_mineta_zcorrected(zertex);
0055   float ohcal_maxeta = get_ohcal_maxeta_zcorrected(zertex);
0056   float minlimit = emcal_mineta;
0057   if (ihcal_mineta > minlimit) minlimit = ihcal_mineta;
0058   if (ohcal_mineta > minlimit) minlimit = ohcal_mineta;
0059   float maxlimit = emcal_maxeta;
0060   if (ihcal_maxeta < maxlimit) maxlimit = ihcal_maxeta;
0061   if (ohcal_maxeta < maxlimit) maxlimit = ohcal_maxeta;
0062   minlimit += jet_radius;
0063   maxlimit -= jet_radius;
0064   return jet_eta < minlimit || jet_eta > maxlimit;
0065 }
0066
0067 int get_scaler(int runnumber, int bit)
0068 {
0069   TSQLServer *db = TSQLServer::Connect("pgsql://sphnxdaqdbreplica:5432/daq","","");
0070
0071   if (db)
0072     {
0073       printf("Server info: %s\n", db->ServerInfo());
0074     }
0075   else
0076     {
0077       printf("bad\n");
0078     }
0079
0080   TSQLRow *row;
0081   TSQLResult *res;
0082   TString cmd = "";
0083   char sql[1000];
0084
0085   sprintf(sql, "select scaled from gl1_scalers where runnumber = %d and index = %d;", runnumber, bit);
0086   res = db->Query(sql);
0087   row = res->Next();
0088   try
0089     {
0090       int retval = stoi(row->GetField(0));
0091       delete row;
0092       delete res;
0093       delete db;
0094       return retval;
0095     }
0096   catch(...)
0097     {
0098       delete row;
0099       delete res;
0100       delete db;
0101       return 0;
0102     }
0103 }
0104 int get_fillnum(int runnumber)
0105 {
0106
0107
0108   TSQLServer *db = TSQLServer::Connect("pgsql://sphnxdbreplica:5432/spinDB","phnxrc","");
0109   if(!db)
0110     {
0111       cout << "FAILED TO CONNECT TO DB!" << endl;
0112       throw exception();
0113     }
0114   TSQLRow *row;
0115   TSQLResult *res;
0116   char sql[1000];
0117
0118
0119
0120   sprintf(sql, "select fillnumber, runnumber from spin where runnumber = %d;", runnumber);
0121   res = db->Query(sql);
0122   if(!res) return -1;
0123   row = res->Next();
0124   if(!row) return -1;
0125
0126   string fillnum(row->GetField(0));
0127   return stoi(fillnum);
0128 }
0129
0130 unsigned int factorial(int n)
0131 {
0132   if(n>10)
0133     {
0134       cout << "Don't do that." << endl;
0135       throw exception();
0136     }
0137   if(n==0) return 1;
0138   if(n==1) return 1;
0139   return n*factorial(n-1);
0140 }
0141
0142 int analyze(int rn, int nseg)
0143 {
0144
0145   ifstream fillnumfile("/sphenix/user/jocl/projects/analysis/LuminosityCounterGoodRuns/run/fillnumbers.txt");
0146   int fills[719][2];
0147
0148   for(int i = 0; i<719; ++i)
0149     {
0150       fillnumfile >> fills[i][0] >> fills[i][1];
0151     }
0152
0153   long long unsigned int firstlivescaler[64] = {0};
0154   long long unsigned int lastlivescaler[64] = {0};
0155   long long unsigned int sumgoodlive[64] = {0};
0156   long long unsigned int totlive[64] = {0};
0157
0158   long long unsigned int firstscaledscaler[64] = {0};
0159   long long unsigned int lastscaledscaler[64] = {0};
0160   long long unsigned int sumgoodscaled[64] = {0};
0161   long long unsigned int totscaled[64] = {0};
0162
0163   long long unsigned int firstrawscaler[64] = {0};
0164   long long unsigned int lastrawscaler[64] = {0};
0165   long long unsigned int sumgoodraw[64] = {0};
0166   long long unsigned int totraw[64] = {0};
0167
0168   long long unsigned int trigcounts[6][64] = {0};
0169   double embdlive[5] = {0};
0170
0171   long long unsigned int nevt = 0;
0172
0173   double avgPS[64] = {0};
0174   long long unsigned int tottrigcounts[6][64] = {0};
0175   double totembdlive[5] = {0};
0176   //TH1D* zhist;
0177   bool gothist = false;
0178   float nmbdc = 0;
0179   float nblair = 0;
0180   float ttseg = 0;
0181   int fillnum = get_fillnum(rn);
0182   for(int i=1; i<nseg+1; ++i)
0183     {
0184       //cout << "start" << endl;
0185       //cout << "test1" << endl;
0186       TFile* file = TFile::Open(("/sphenix/user/jocl/projects/trigcount_files/"+to_string(rn)+"/triggercounter_5z_"+to_string(rn)+"_"+to_string(i)+".root").c_str());
0187       //cout << "test2" << endl;
0188       //cout << "gettree" << endl;
0189       TTree* tree = (TTree*)file->Get("_tree");
0190       //cout << "test3" << endl;
0191       //cout << "gottree" << endl;
0192       long long unsigned int segloraw[64] = {0};
0193       long long unsigned int seghiraw[64] = {0};
0194       long long unsigned int segloscaled[64] = {0};
0195       long long unsigned int seghiscaled[64] = {0};
0196       long long unsigned int seglolive[64] = {0};
0197       long long unsigned int seghilive[64] = {0};
0198       long long unsigned int startBCO = 0;
0199       long long unsigned int endBCO = 0;
0200       int badflag = 0;
0201
0202       long long unsigned int segevt = 0;
0203       //cout << "setbranches" << endl;
0204       //cout << "test4" << endl;
0205       tree->SetBranchAddress("badFlag",&badflag);
0206       tree->SetBranchAddress("startScal",segloscaled);
0207       tree->SetBranchAddress("endScal",seghiscaled);
0208       tree->SetBranchAddress("startRaw",segloraw);
0209       tree->SetBranchAddress("endRaw",seghiraw);
0210       tree->SetBranchAddress("startLive",seglolive);
0211       tree->SetBranchAddress("endLive",seghilive);
0212       tree->SetBranchAddress("nevt",&segevt);
0213       tree->SetBranchAddress("trigCounts",trigcounts);
0214       tree->SetBranchAddress("eMBDlive",embdlive);
0215       tree->SetBranchAddress("startBCO",&startBCO);
0216       tree->SetBranchAddress("endBCO",&endBCO);
0217       //cout << "test5" << endl;
0218       //cout << "getentry" << endl;
0219       tree->GetEntry(0);
0220       //cout << "test6" << endl;
0221       //cout << "gotentry" << endl;
0222       float tSeg = (endBCO-startBCO)/(9.4e6);
0223       ttseg += tSeg;
0224       if(tSeg == 0)
0225     {
0226       cout << "rn " << rn << " seg " << i << " tseg = 0" << endl;
0227       continue;
0228     }
0229       //cout << endBCO << " " << startBCO << endl;
0230       //cout << badflag << endl;
0231       nevt += segevt;
0232
0233       //cout << "test7" << endl;
0234       int nBunch = -1;
0235       for(int i=0; i<719; ++i)
0236         {
0237           if(fillnum == fills[i][0])
0238             {
0239               nBunch = fills[i][1];
0240               break;
0241             }
0242         }
0243
0244       if(fillnum == -999 || nBunch == -1)
0245     {
0246       nBunch = 111;
0247     }
0248
0249       if(nBunch != 111)
0250     {
0251       cout << "NBunch/rn: " << nBunch << " " << rn << endl;
0252     }
0253
0254       if((rn >= 48863 && rn <= 48867) || (rn <= 51508 && rn >= 51506)) nBunch = 111;
0255       if(i==1)
0256     {
0257       for(int j=0; j<64; ++j)
0258         {
0259           firstrawscaler[j] = segloraw[j];
0260           firstlivescaler[j] = seglolive[j];
0261           firstscaledscaler[j] = segloscaled[j];
0262         }
0263     }
0264       for(int j=0; j<64; ++j)
0265     {
0266       if(seghilive[j] != 0)
0267         {
0268           lastrawscaler[j] = seghiraw[j];
0269           lastlivescaler[j] = seghilive[j];
0270           lastscaledscaler[j] = seghiscaled[j];
0271         }
0272     }
0273       for(int j=0; j<5; ++j)
0274     {
0275       totembdlive[j] += embdlive[j];
0276     }
0277       float avgTot[64] = {0};
0278       float rB = nBunch*78.2e3; //nBunch obtained from CAD tarball provided by Kin Yip
0279       float rM = (seghiraw[10]-segloraw[10])/tSeg; //tSeg from BCO difference between start and end of segment * beam clock (9.4 MHz)
0280       if(rM > rB)
0281     {
0282       cout << "rM > rB!! " << rB << " " << rM << " " << tSeg << endl;
0283     }
0284       if(rM > 0.82*rB)
0285     {
0286       cout << "rM too high!" << rB << " " << rM << " " << tSeg << endl;
0287     }
0288       float pSum = 0;
0289       for(int k=1; k<10; ++k)
0290     {
0291       pSum += pow(-log(1-rM/rB),k) / factorial(k-1); //sum probability * number of collisions for rate
0292     }
0293       //cout << endl;
0294       for(int j=0; j<64; ++j)
0295     {
0296       for(int k=0; k<6; ++k)
0297         {
0298           tottrigcounts[k][j] += trigcounts[k][j]; //number of triggers actually in segment for three zvtx selections, all trigger bits
0299         }
0300
0301       sumgoodscaled[j] += (seghiscaled[j] - segloscaled[j]); //scaled counts as taken from difference of scaled scaler at start and end of segment
0302       sumgoodraw[j] += (seghiraw[j] - segloraw[j]);
0303       sumgoodlive[j] += (seghilive[j] - seglolive[j]); //use regular live counts for all triggers except 10
0304       if(j==10)
0305         {
0306           nmbdc += tSeg * (rB-rM) * pSum; //N_{coll} for trigger 10
0307           nblair += (seghiraw[j]-segloraw[j])*(1-log(1-rM/rB));
0308         }
0309     }
0310       //cout << i << " " << sumgoodscaled[10] << " " << sumgoodscaled[18] << endl;
0311       if(i%10 == 0)
0312     {
0313       cout << "Segment " << i << " " << tSeg << " " << (rB - rM) << " " << pSum << " " << seghiraw[10] - segloraw[10] << endl;
0314     }
0315       //cout << "Run " << rn << " Segment " << i << " time: " << tSeg << " beam rate: " << rB << " bunches: " << nBunch <<" MBD live rate: " << rM << " Sum(n*p(n)): " << pSum << " lambda: " << -log(1-rM/rB) << endl
0316       //cout << i << endl;
0317       file->Close();
0318       //cout << "closed" << endl;
0319     }
0320
0321   cout << "sumgoodraw / nmbdc " << sumgoodraw[10] << " " << nmbdc << endl;
0322
0323   //int totalgood = zhist->Integral();
0324   for(int i=0; i<64; ++i)
0325     {
0326       avgPS[i] = sumgoodscaled[i]==0?0:((double)sumgoodlive[i])/sumgoodscaled[i]; //calculate average prescale for each trigger
0327       totlive[i] = lastlivescaler[i] - firstlivescaler[i];
0328       totscaled[i] = lastscaledscaler[i] - firstscaledscaler[i];
0329       totraw[i] = lastrawscaler[i] - firstrawscaler[i];
0330     }
0331
0332
0333   TFile* outfile = TFile::Open(("triggeroutput_nblair_"+to_string(rn)+".root").c_str(),"RECREATE");
0334   outfile->cd();
0335   TTree* outt = new TTree("outt","A persevering date tree");
0336
0337   outt->Branch("totlive",totlive,"totlive[64]/l");
0338   outt->Branch("totscaled",totscaled,"totscaled[64]/l");
0339   outt->Branch("totraw",totraw,"totraw[64]/l");
0340   outt->Branch("sumgoodscaled",sumgoodscaled,"sumgoodscaled[64]/l");
0341   outt->Branch("sumgoodraw",sumgoodraw,"sumgoodraw[64]/l");
0342   outt->Branch("sumgoodlive",sumgoodlive,"sumgoodlive[64]/l");
0343   outt->Branch("avgPS",avgPS,"avgPS[64]/D");
0344   outt->Branch("tottrigcounts",tottrigcounts,"tottrigcounts[6][64]/l");
0345   outt->Branch("totembdlive",totembdlive,"totembdlive[5]/D");
0346   outt->Branch("nevt",&nevt,"nevt/l");
0347   outt->Branch("nmbdc",&nmbdc,"nmbdc/F");
0348   outt->Branch("nblair",&nblair,"nblair/F");
0349   outt->Branch("ttseg",&ttseg,"ttseg/F");
0350   outt->Fill();
0351   outt->Write();
0352   //zhist->Write();
0353   outfile->Close();
0354   return 0;
0355 }