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, int clt)
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 nocorN = 0;
0180   float nblair = 0;
0181   float ttseg = 0;
0182   int fillnum = get_fillnum(rn);
0183   TFile* ingraph = TFile::Open("mbd_to_col_map.root");
0184   TGraph* mbdcolmap = (TGraph*)ingraph->Get("mbdtogr1");
0185   mbdcolmap->SetBit(TGraph::kIsSortedX);
0186   for(int i=1; i<nseg+1; ++i)
0187     {
0188       //cout << "start" << endl;
0189       //cout << "test1" << endl;
0190       TFile* file = TFile::Open(("/sphenix/user/jocl/projects/trigcount_files/"+to_string(rn)+"/triggercounter_5z_"+to_string(rn)+"_"+to_string(i)+(clt?"_clt":"")+".root").c_str());
0191
0192       //cout << "test2" << endl;
0193       //cout << "gettree" << endl;
0194       TTree* tree = (TTree*)file->Get("_tree");
0195       //cout << "test3" << endl;
0196       //cout << "gottree" << endl;
0197       long long unsigned int segloraw[64] = {0};
0198       long long unsigned int seghiraw[64] = {0};
0199       long long unsigned int segloscaled[64] = {0};
0200       long long unsigned int seghiscaled[64] = {0};
0201       long long unsigned int seglolive[64] = {0};
0202       long long unsigned int seghilive[64] = {0};
0203       long long unsigned int startBCO = 0;
0204       long long unsigned int endBCO = 0;
0205       int badflag = 0;
0206
0207       long long unsigned int segevt = 0;
0208       //cout << "setbranches" << endl;
0209       //cout << "test4" << endl;
0210       tree->SetBranchAddress("badFlag",&badflag);
0211       tree->SetBranchAddress("startScal",segloscaled);
0212       tree->SetBranchAddress("endScal",seghiscaled);
0213       tree->SetBranchAddress("startRaw",segloraw);
0214       tree->SetBranchAddress("endRaw",seghiraw);
0215       tree->SetBranchAddress("startLive",seglolive);
0216       tree->SetBranchAddress("endLive",seghilive);
0217       tree->SetBranchAddress("nevt",&segevt);
0218       tree->SetBranchAddress("trigCounts",trigcounts);
0219       tree->SetBranchAddress("eMBDlive",embdlive);
0220       tree->SetBranchAddress("startBCO",&startBCO);
0221       tree->SetBranchAddress("endBCO",&endBCO);
0222       //cout << "test5" << endl;
0223       //cout << "getentry" << endl;
0224       tree->GetEntry(0);
0225       //cout << "test6" << endl;
0226       //cout << "gotentry" << endl;
0227       float tSeg = (endBCO-startBCO)/(9.4e6);
0228       ttseg += tSeg;
0229       if(tSeg == 0)
0230     {
0231       cout << "rn " << rn << " seg " << i << " tseg = 0" << endl;
0232       continue;
0233     }
0234       //cout << endBCO << " " << startBCO << endl;
0235       //cout << badflag << endl;
0236       nevt += segevt;
0237
0238       //cout << "test7" << endl;
0239       int nBunch = -1;
0240       for(int i=0; i<719; ++i)
0241         {
0242           if(fillnum == fills[i][0])
0243             {
0244               nBunch = fills[i][1];
0245               break;
0246             }
0247         }
0248
0249       if(fillnum == -999 || nBunch == -1)
0250     {
0251       nBunch = 111;
0252     }
0253
0254       if(nBunch != 111)
0255     {
0256       cout << "NBunch/rn: " << nBunch << " " << rn << endl;
0257     }
0258
0259       if((rn >= 48863 && rn <= 48867) || (rn <= 51508 && rn >= 51506) || rn==52444) nBunch = 111;
0260       if(i==1)
0261     {
0262       for(int j=0; j<64; ++j)
0263         {
0264           firstrawscaler[j] = segloraw[j];
0265           firstlivescaler[j] = seglolive[j];
0266           firstscaledscaler[j] = segloscaled[j];
0267         }
0268     }
0269       for(int j=0; j<64; ++j)
0270     {
0271       if(seghilive[j] != 0)
0272         {
0273           lastrawscaler[j] = seghiraw[j];
0274           lastlivescaler[j] = seghilive[j];
0275           lastscaledscaler[j] = seghiscaled[j];
0276         }
0277     }
0278       for(int j=0; j<5; ++j)
0279     {
0280       totembdlive[j] += embdlive[j];
0281     }
0282       float avgTot[64] = {0};
0283       float rB = nBunch*78.2e3; //nBunch obtained from CAD tarball provided by Kin Yip
0284       float rM = (seghiraw[10]-segloraw[10])/tSeg; //tSeg from BCO difference between start and end of segment * beam clock (9.4 MHz).
0285       rM = mbdcolmap->Eval(rM);
0286       if(rM > rB)
0287     {
0288       cout << "rM > rB!! " << rB << " " << rM << " " << tSeg << endl;
0289     }
0290       if(rM > 0.82*rB)
0291     {
0292       cout << "rM too high!" << rB << " " << rM << " " << tSeg << endl;
0293     }
0294       float pSum = 0;
0295       for(int k=1; k<10; ++k)
0296     {
0297       pSum += pow(-log(1-rM/rB),k) / factorial(k-1); //sum probability * number of collisions for rate
0298     }
0299       //cout << endl;
0300       float ncorrseg = 0;
0301       for(int j=0; j<64; ++j)
0302     {
0303       for(int k=0; k<6; ++k)
0304         {
0305           tottrigcounts[k][j] += trigcounts[k][j]; //number of triggers actually in segment for three zvtx selections, all trigger bits
0306         }
0307
0308       sumgoodscaled[j] += (seghiscaled[j] - segloscaled[j]); //scaled counts as taken from difference of scaled scaler at start and end of segment
0309       sumgoodraw[j] += (seghiraw[j] - segloraw[j]);
0310       sumgoodlive[j] += (seghilive[j] - seglolive[j]); //use regular live counts for all triggers except 10
0311       if(j==10)
0312         {
0313           nmbdc += tSeg * (rB-rM) * pSum; //N_{coll} for trigger 10
0314           ncorrseg = tSeg * (rB-rM) * pSum;
0315
0316           nblair += (seghiraw[j]-segloraw[j])*(1-log(1-rM/rB));
0317         }
0318     }
0319       //cout << i << " " << sumgoodscaled[10] << " " << sumgoodscaled[18] << endl;
0320       /*
0321       if(i%10 == 0)
0322     {
0323       cout << "Segment " << i << " " << tSeg << " " << (rB - rM) << " " << pSum << " " << seghiraw[10] - segloraw[10] << endl;
0324     }
0325       */
0326       cout << "Run " << rn << " Segment " << i << " time: " << tSeg << " beam rate: " << rB << " bunches: " << nBunch <<" True collision rate: " << rM << " N_{corrected}: " << ncorrseg << " N_{uncorrected}: " << rM*tSeg << " ???: " << (seghiraw[10]-segloraw[10])*42/25.2 << endl;
0327       //cout << i << endl;
0328       nocorN += tSeg*rM;
0329       file->Close();
0330       //cout << "closed" << endl;
0331     }
0332   cout << "uncorrected / corrected total collisions " << nocorN << " " << nmbdc << endl;
0333
0334   //int totalgood = zhist->Integral();
0335   for(int i=0; i<64; ++i)
0336     {
0337       avgPS[i] = sumgoodscaled[i]==0?0:((double)sumgoodlive[i])/sumgoodscaled[i]; //calculate average prescale for each trigger
0338       totlive[i] = lastlivescaler[i] - firstlivescaler[i];
0339       totscaled[i] = lastscaledscaler[i] - firstscaledscaler[i];
0340       totraw[i] = lastrawscaler[i] - firstrawscaler[i];
0341     }
0342
0343
0344   TFile* outfile = TFile::Open(("triggeroutput_nblair_"+to_string(rn)+(clt?"_clt":"")+".root").c_str(),"RECREATE");
0345   outfile->cd();
0346   TTree* outt = new TTree("outt","A persevering date tree");
0347
0348   outt->Branch("totlive",totlive,"totlive[64]/l");
0349   outt->Branch("totscaled",totscaled,"totscaled[64]/l");
0350   outt->Branch("totraw",totraw,"totraw[64]/l");
0351   outt->Branch("sumgoodscaled",sumgoodscaled,"sumgoodscaled[64]/l");
0352   outt->Branch("sumgoodraw",sumgoodraw,"sumgoodraw[64]/l");
0353   outt->Branch("sumgoodlive",sumgoodlive,"sumgoodlive[64]/l");
0354   outt->Branch("avgPS",avgPS,"avgPS[64]/D");
0355   outt->Branch("tottrigcounts",tottrigcounts,"tottrigcounts[6][64]/l");
0356   outt->Branch("totembdlive",totembdlive,"totembdlive[5]/D");
0357   outt->Branch("nevt",&nevt,"nevt/l");
0358   outt->Branch("nmbdc",&nmbdc,"nmbdc/F");
0359   outt->Branch("nocorN",&nocorN,"nocorN/F");
0360   outt->Branch("nblair",&nblair,"nblair/F");
0361   outt->Branch("ttseg",&ttseg,"ttseg/F");
0362   outt->Fill();
0363   outt->Write();
0364   //zhist->Write();
0365   outfile->Close();
0366   return 0;
0367 }