Robert Michaels, rom@jlab.org, Jefferson
Lab Hall A, May 2003
//-------------------------------------------------------- // tscalevt_main.C // // Test of the ROC10/11 scalers which are read in the datastream // every synch event (typ. every 100 events). // // R. Michaels, Jan 2002 //-------------------------------------------------------- #include#include #include "THaScaler.h" #include "THaCodaFile.h" #include "THaEvData.h" #ifndef __CINT__ #include "TROOT.h" #include "TFile.h" #include "TH1.h" #include "TH2.h" #include "TProfile.h" #include "TNtuple.h" #include "TRandom.h" #endif #define DEBUG 0 #define SKIPEVENT 0 // How many scaler events to skip // (This reduces effect of clock granularity) #define MYROC 11 // Needed for the SKIPEVENT trick #define BCM_CUT1 3000 // cut on BCM (x1 gain) to require beam on. #define BCM_CUT3 10000 // cut on BCM (x3 gain) to require beam on. #define BCM_CUT10 30000 // cut on BCM (x10 gain) to require beam on. #define NBCM 6 // number of BCM signals using namespace std; int main(int argc, char* argv[]) { int i,iev,evnum,trig,iskip,status; Int_t clkp, clkm, lastclkp, lastclkm; Float_t sum,asy; TString filename = "run.dat"; char bank[100] = "EvLeft"; // Event stream, Left HRS cout << "Analyzing CODA file " << filename << endl; cout << "Bank = " << bank << endl << flush; evnum = 10000000; if (argc > 1) evnum = atoi(argv[1]); THaScaler *scaler = new THaScaler(bank); if (scaler->Init("1-9-2002") == -1) { cout << "Error initializing scaler " << endl; return 1; } THaCodaFile *coda = new THaCodaFile(filename); THaEvData evdata; // Pedestals. Left, Right Arms. u1,u3,u10,d1,d3,d10 Float_t bcmpedL[NBCM] = { 188.2, 146.2, 271.6, 37.8, 94.2, 260.2 }; Float_t bcmpedR[NBCM] = { 53.0, 41.8, 104.1, 0., 101.6, 254.6 }; Float_t bcmped[NBCM]; if (strcmp(bank,"EvLeft") == 0) { cout << "Using Left Arm BCM pedestals" << endl; for (i = 0; i < NBCM; i++) { bcmped[i] = bcmpedL[i]; } } else { cout << "Using Right Arm BCM pedestals" << endl; for (i = 0; i < NBCM; i++) { bcmped[i] = bcmpedR[i]; } } // Initialize root and output TROOT scalana("scalroot","Hall A scaler analysis"); TFile hfile("scaler.root","RECREATE","Scaler data in Hall A"); // Define the ntuple here // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 char rawrates[]="time:u1:u3:u10:d1:d3:d10:t1:t2:t3:t4:t5:clkp:clkm:tacc:"; // 15 16 17 18 19 20 21 22 23 24 25 26 char asymmetries[]="au1:au3:au10:ad1:ad3:ad10:at1:at2:at3:at4:at5:aclk"; int nlen = strlen(rawrates) + strlen(asymmetries); char *string_ntup = new char[nlen+1]; strcpy(string_ntup,rawrates); strcat(string_ntup,asymmetries); TNtuple *ntup = new TNtuple("ascal","Scaler Rates",string_ntup); Float_t farray_ntup[27]; // Note, dimension is same as size of string_ntup (i.e. nlen+1) status = 0; iev = 0; iskip = 0; lastclkp = 0; lastclkm = 0; while (status == 0) { status = coda->codaRead(); if (status != 0) { cout << "coda status nonzero. assume EOF" << endl; goto quit; } evdata.LoadEvent(coda->getEvBuffer()); // Dirty trick to average over larger time intervals (depending on // SKIPEVENT) to reduce the fluctuations due to clock. if (evdata.GetRocLength(MYROC) > 16) iskip++; if (SKIPEVENT != 0 && iskip < SKIPEVENT) continue; iskip = 0; scaler->LoadData(evdata); // Not every trigger has new scaler data, so skip if not new. if ( !scaler->IsRenewed() ) continue; if (iev++ > evnum) goto quit; // Fill the ntuple here // Note, we must average the two helicities here to get non-helicity rates Double_t time = (scaler->GetPulser(1,"clock") + scaler->GetPulser(-1,"clock"))/1024; farray_ntup[0] = time; farray_ntup[1] = 0.5*(scaler->GetBcmRate(1,"bcm_u1") + scaler->GetBcmRate(-1,"bcm_u1")); farray_ntup[2] = 0.5*(scaler->GetBcmRate(1,"bcm_u3") + scaler->GetBcmRate(-1,"bcm_u3")); farray_ntup[3] = 0.5*(scaler->GetBcmRate(1,"bcm_u10") + scaler->GetBcmRate(-1,"bcm_u10")); farray_ntup[4] = 0.5*(scaler->GetBcmRate(1,"bcm_d1") + scaler->GetBcmRate(-1,"bcm_d1")); farray_ntup[5] = 0.5*(scaler->GetBcmRate(1,"bcm_d3") + scaler->GetBcmRate(-1,"bcm_d3")); farray_ntup[6] = 0.5*(scaler->GetBcmRate(1,"bcm_d10") + scaler->GetBcmRate(-1,"bcm_d10")); for (trig = 1; trig <= 5; trig++) { farray_ntup[6+trig] = 0.5*(scaler->GetTrigRate(1,trig) + scaler->GetTrigRate(-1,trig)); } clkp = scaler->GetNormData(1,"clock"); clkm = scaler->GetNormData(-1,"clock"); farray_ntup[12] = clkp - lastclkp; farray_ntup[13] = clkm - lastclkm; lastclkp = clkp; lastclkm = clkm; farray_ntup[14] = 0.5*(scaler->GetNormRate(1,"TS-accept") + scaler->GetNormRate(-1,"TS-accept")); string bcms[] = {"bcm_u1", "bcm_u3", "bcm_u10", "bcm_d1", "bcm_d3", "bcm_d10"}; for (int ibcm = 0; ibcm < 6; ibcm++ ) { sum = scaler->GetBcmRate(1,bcms[ibcm].c_str()) - bcmped[ibcm] + scaler->GetBcmRate(-1,bcms[ibcm].c_str()) - bcmped[ibcm]; asy = -999; if (sum != 0) { asy = (scaler->GetBcmRate(1,bcms[ibcm].c_str()) - scaler->GetBcmRate(-1,bcms[ibcm].c_str())) / sum; } farray_ntup[15+ibcm] = asy; } for (trig = 1; trig <= 5; trig++) { asy = -999; if (scaler->GetBcmRate(1,"bcm_u3") > BCM_CUT3 && scaler->GetBcmRate(-1,"bcm_u3") > BCM_CUT3) { sum = scaler->GetTrigRate(1,trig)/scaler->GetBcmRate(1,"bcm_u3") + scaler->GetTrigRate(-1,trig)/scaler->GetBcmRate(-1,"bcm_u3"); if (sum != 0) { asy = (scaler->GetTrigRate(1,trig)/scaler->GetBcmRate(1,"bcm_u3") - scaler->GetTrigRate(-1,trig)/scaler->GetBcmRate(-1,"bcm_u3")) / sum; } } farray_ntup[20+trig] = asy; } sum = scaler->GetPulser(1,"clock") + scaler->GetPulser(-1,"clock"); asy = -999; if (sum != 0) { asy = (scaler->GetPulser(1,"clock") - scaler->GetPulser(-1,"clock")) / sum; } farray_ntup[26] = asy; if (DEBUG) { cout << "event type " << evdata.GetEvType() << " event number " << evdata.GetEvNum() << endl; scaler->Print(); if (farray_ntup[12] < 1000) cout << "Low clock" << endl; cout << " clock(+) " << scaler->GetNormRate(1,"clock"); cout << " clock(-) " << scaler->GetNormRate(-1,"clock") << endl; if (farray_ntup[16]>0.005||farray_ntup[16]<-0.005) cout << "Big asy" << endl; } ntup->Fill(farray_ntup); } quit: hfile.Write(); hfile.Close(); exit(0); }