///////////////////////////////////////////////////////////////////////////// // // \author Jason C. Webb (jwebb@iucf.indiana.edu) // ///////////////////////////////////////////////////////////////////////////// // // A root script to verify that everything is fine (or at least // not obviously out of whack) with the EEmc towers in each // run. Basically, we test the adc spectra from each run in each // tower against a set of "reference" histograms (most likely a // "good" run. This reference file should be specified in the // call to checkTowers, otherwise it selects one AT RANDOM from // the current directory (actually, it selects the first file in // the directory, but you see the point of setting the file.) // // This code is slow to execute. You have been warned. // // Note that this code is not really specific to the problem of finding // problems with the towers during the experiment. One could generalize // this script to compare any sets of histograms in any files (so long // as the directory structure remains the same between them). // // PREPROCESSOR DIRECTIVES ////////////////////////////////////////////////// // // If defined, output will be a comma separated list of 4 bits, each // corresponding to the order of the test (chi^2, ped shift, width // and NULL adc values). // //#define __MACHINE_READABLE__ // ///////////////////////////////////////////////////////////////////////////// // // Usage: // // First, generate pedestal histogram files as described here: // // http://www.star.bnl.gov/eemc/how_to/doTwPed.html // // (Just use the accumTwPed.C script, but make sure that the trigger // conditions are reasonable.) // // Next, cd to the directory containing the output ".hist.root" files // and execute this script, either interactively -- // // .x checkRuns(); // to run with defaults // // .x checkRuns( "file.hist.root", // reference file (defaults to 1st file) // 2.0, // cut on chi^2 / ndf // 1.0, // max pedestal shift // 0.3, // max increase in rms // 100 ) // max bin to consider // // or in batch mode // // $ root -q -b checkRuns.C // // $ root -q -b 'checkRuns.C("file.hist.root",...)' // // The result will be a table with run number (e.g. filename) in the first // column, followed by a comma separated list of tower errors in that // run (the first row lists the towers). Note that this file output is // designed to be read in by a spreadsheet or parsed by a scripting // language. // ///////////////////////////////////////////////////////////////////////////// // // Revision History // // $Log$ // void checkRuns( //TString *ref = NULL, Char_t *cref = NULL, //"R4123132.ezTree-ped.hist.root", Double_t chi2cut = 2.0, Double_t maxShift = 1.1, Double_t maxWidthChange = 0.3, Int_t maxBin = 70 ) { TString *ref = NULL; if ( cref ) ref = new TString(cref); ////////////////////////////////////////////////////////////////// // // Print out a legend for the user // #ifndef __MACHINE_READABLE__ std::cout << "LEGEND:" << endl; std::cout << "X - chi^2 between run and reference > " << chi2cut << endl; std::cout << "P - pedestal shift by > " << maxShift << endl; std::cout << "W - change in width by > " << maxWidthChange * 100. << "%" << endl; std::cout << "0 - mean < channel 2" << endl << endl; std::cout << "maximum ADC considererd = " << maxBin << endl << endl; #else std::cout << "LEGEND:" << endl; std::cout << "1000 - chi^2 between run and reference > " << chi2cut << endl; std::cout << "0100 - pedestal shift by > " << maxShift << endl; std::cout << "0010 - change in width by > " << maxWidthChange * 100. << "%" << endl; std::cout << "0001 - mean < channel 2" << endl << endl; std::cout << "maximum ADC considererd = " << maxBin << endl << endl; #endif std::cout << endl; ////////////////////////////////////////////////////////////////// // // Get a list of the system files in the current directory // and loop over them // TList *sysFiles = TSystemDirectory("pwd",".").GetListOfFiles(); TIter next(sysFiles); TSystemFile *file; TFile *refFile = NULL; TList *refHistKeys; TString refName = ""; TList *refHistos = new TList(); // List of histograms TH1F hmask; Int_t hmask_set = 0; while ( file = (TSystemFile *)next() ) { // Abort if this is not a ".hist.root" file if ( !TString(file -> GetName()).EndsWith(".hist.root") ) continue; ////////////////////////////////////////////////////////////////// // // Code which only executes for the reference file, or the first // file in the directory if the reference is not defined. // // If the reference file isn't defined, use the first ".hist.root" // file encountered. if ( !ref ) { ref = new TString(file -> GetName()); // std::cout << "reference file: " << ref -> Data() << endl; } // If the reference file has yet to be connected as a TFile, do so // and load in the histograms from that file if ( !refFile ) { refFile = new TFile(ref->Data()); refFile -> cd("pedTw"); refHistKeys = gDirectory -> GetListOfKeys(); refName = TString(refFile -> GetName()); refName.Resize( refName.First(".") ); TKey *k; TIter nextkey(refHistKeys); // Read in the histograms and stuff them into a TList, so // we only need to read these histograms in one time. Int_t idx = 0; std::cout << " run file / tower, "; while ( k = (TKey *)nextkey() ) { TH1F *h = (TH1F *)k -> ReadObj() ; // Create a histogram to "mask off" bins 100 and higher // via histogram multiplication. if ( !hmask_set ) { hmask = (TH1F)h -> Clone("hmask"); hmask -> Reset(); for ( int ii = 0;ii Fill(ii); hmask_set++; } // zero out bins above adc = 100 h -> Multiply(&hmask); refHistos -> Add ( h ); // And output a comma separated list of the tower "names" std::cout << h -> GetName() << ", "; } std::cout << "END" << endl; } // If the file being considered is the reference file, punt if ( TString(file -> GetName()) == *ref ) continue; // Connect the new root file and prepare to loop over // histograms. TFile *testFile = new TFile ( file -> GetName() ); testFile -> cd("pedTw"); TList *testHistKeys = gDirectory -> GetListOfKeys(); // Output the filename as the first entry in the row (NO endl here) std::cout << file -> GetName() << "," << std::flush; // Loop over the reference histograms, and perform various tests // of the shape and agreement with the reference histograms TH1F *refHist; TIter nexthist(refHistos); Int_t idx = 0; while ( refHist = (TH1F *)nexthist() ) { // Read in the test histogram. WARNING: this assumes that // the lists for the reference and test histograms are // filled in a sensible way. IF THIS IS NOT THE CASE, // you need to use the (slower) FindObject( refHist -> GetName() ) // method. TKey *testKey = (TKey *)testHistKeys -> At(idx++); TH1F *testHist = (TH1F *)testKey -> ReadObj(); // Zero out bins above adc = 100 testHist -> Multiply(&hmask); /////////////////////////////////////////////////////// // // Here we test the compatability of the two histograms // using a chi^2 test // Double_t myChi2 = 0.0; Double_t myNdf = 0.0; Double_t rr = refHist -> Integral(); Double_t tt = testHist -> Integral(); Double_t aa = sqrt( rr / tt ); for ( int ii = 0; ii < maxBin; ii++ ) { Double_t yyRef = refHist -> GetBinContent(ii); Double_t yyTest = testHist -> GetBinContent(ii); if ( yyRef + yyTest != 0.0 ) { myChi2 += ( yyRef / aa - yyTest * aa )**2 / ( yyRef + yyTest ); myNdf++; } } if ( myNdf != 0.0 ) { myChi2 /= myNdf; } #ifndef __MACHINE_READABLE__ if ( myChi2 > chi2cut ) std::cout << "X" << std::flush; #else std::cout << ( myChi2 > chi2cut ) ? 0 : 1 << std::flush; #endif /////////////////////////////////////////////////////// // // Next we check for a shift in the mean // refHist ->Fit("gaus","ORQ","",0.,maxBin); TF1 * refFit=refHist->GetFunction("gaus"); testHist ->Fit("gaus","ORQ","",0.,maxBin); TF1 * testFit=testHist->GetFunction("gaus"); Double_t mmRef = refFit->GetParameter(1); //refHist -> GetMean(); Double_t mmTest = testFit->GetParameter(1); //testHist -> GetMean(); //printf("%p %p \n",refFit, testFit); //printf("%s %f %f\n",refHist ->GetName(), refFit->GetParameter(1),); #ifndef __MACHINE_READABLE__ if ( fabs( mmRef - mmTest ) > maxShift ) std::cout << "P" << std::flush; #else std::cout << ( fabs( mmRef - mmTest ) > maxShift ) ? 0 : 1 << std::flush; #endif /////////////////////////////////////////////////////// // // Check for a change in the widths of the distributions // Double_t rmsRef = refFit->GetParameter(2); //refHist -> GetRMS(); Double_t rmsTest =testFit->GetParameter(2); // testHist -> GetRMS(); #ifndef __MACHINE_READABLE__ if ( fabs ( rmsRef - rmsTest ) > maxWidthChange * rmsRef ) std::cout << "W" << std::flush; #else std::cout << ( fabs ( rmsRef - rmsTest ) > maxWidthChange * rmsRef ) ? 0 : 1 << std::flush; #endif /////////////////////////////////////////////////////// // // Check for a lot of NULL readings // #ifndef __MACHINE_READABLE__ if ( mmRef < 2.0 ) std::cout << "0" << std::flush; #else std::cout << ( mmRef < 2.0 ) ? 0 : 1 << std::flush; #endif // Comma-separated list of errors in each tower std::cout << "," << std::flush; // Cleanup for loop over towers in a given run. Here we // delete any objects created with the "new" command // which we don't wan't to keep around... // (nothing to delete... good). } // Loop over towers in a given run std::cout << "END" << endl; std::cout << std::flush; // Cleanup, delete pointers we no longer need, etc... delete testFile; } // Loop over all files in the pwd of the filesystem }