Int_t read_parms (Int_t j, TString filename, Float_t *yparms, Float_t *yparms_err ) {    

  vector <TString> sdme;
  vector <double> parms;
  vector <double> parms_err;


    TString infile = filename+".fit2";   // file with parameters
    
    cout << "Opening parameters file: " << infile.Data() << endl;
  
    // now read and print fitted values
    ifstream parameters;
    parameters.open (infile.Data());
    if (!parameters) {
        cout << "ERROR: Failed to open data file= " << infile.Data() << endl;
        return 1;   // failed to open data file
    }

    TString CanvasName = filename;
    // cout << "j=" << j << " CanvasName=" << CanvasName << endl;
    c[j] = new TCanvas(CanvasName,CanvasName,200,10,1000,700);
    
    TString line;
    while (line.ReadLine(parameters)){
        
        TObjArray *tokens = line.Tokenize("\t");
        Int_t ntokens = tokens->GetEntries();
        
        // cout << " ntokens=" << ntokens << " line=" << line.Data() << endl;
	Int_t jmax=ntokens/3;
        for (Int_t j=0; j<jmax; j++){
	  sdme.push_back( (((TObjString*)tokens->At(3*j))->GetString()) );
          parms.push_back( (((TObjString*)tokens->At(3*j+1))->GetString()).Atof() );
          parms_err.push_back( (((TObjString*)tokens->At(3*j+2))->GetString()).Atof());
        }
        
    }   // end loop over lines
    
        TString title = filename;
        TLatex *t1 = new TLatex(0.2,0.85,title);
        t1->SetNDC();
        t1->SetTextSize(0.04);
        t1->Draw();
        
        for (Int_t j=0; j<sdme.size()-2; j++) {     // -2 to eliminate Sigma and P
	  cout << sdme[j] << "=" << parms[j] << " err=" << parms_err[j] << endl;
        
            TString sdmename;
            sdmename = sdme[j];
            title.Form("%s = \t%.3f #pm %.3f\n",sdmename.Data(),parms[j],parms_err[j]);
	    // cout << "title=" << title << endl;
        	TLatex *t1 = new TLatex(0.2,0.75 - 0.05*j,title);
        	t1->SetNDC();
        	t1->SetTextSize(0.04);
        	t1->Draw();
        }
        

	// cout << " parms.size()=" << parms.size() << endl;
      for (jj=0; jj<parms.size()-2; jj++) {
	yparms[jj] = parms[jj];
	yparms_err[jj] = parms_err[jj];
      }

    
    parameters.close();

    return 0;
    
}

int convertToInt(const TObjArray& run_ls_event, int idx)
{
  TObjString* entry = dynamic_cast<TObjString*>(run_ls_event.At(idx));
  assert(entry);
  return atoi(entry->GetString().Data());
}

int64_t gr_engine (const char* name, int rx, int ry, int rdx, int rdy)
{
 double x[64000],y[64000],dy[64000],dx[64000],  bu;
 int i,j;
 FILE * pFile;
 int MAXLINES=64000;
 // char mystring[1500];// one line
 TString oneline, title=name, token;

  printf("i... gr_engine: going to open filename=%s\n", name );

  pFile=fopen( name ,"r" ); 
  if (pFile==NULL) {
    printf("cannot open %s,STOPping\n", name ); 
    return 0;
  } // error
  //  printf("file opened\n%s","");
  //.......................readout HERE.......i
  i=0;
  int lastlen;// remove spaces
  while ((i<MAXLINES)&&( feof(pFile)==0) ){
    //    printf("  reading line %d...\n", i );
    if ( oneline.Gets (pFile, kTRUE) ){//chop true ... continue if not eof
      //      printf("     string==%s\n",  oneline.Data() ); 

    //purify .................. start
    do {
      lastlen=oneline.Length();
      if (oneline.Index(" ")==0){oneline=oneline(1,oneline.Length()-1);}
      lastlen=oneline.Length();
      if (oneline.Index(" ")==lastlen){oneline=oneline(0,oneline.Length()-1);}
       oneline.ReplaceAll("\t"," ");
       oneline.ReplaceAll("  "," ");
       oneline.ReplaceAll("  "," ");
       //       printf(" ---    string==<%s> (%d)\n",  oneline.Data(), lastlen ); 
    }while( lastlen!=oneline.Length());
    //    printf("     string==<%s> (%d)\n",  oneline.Data(), lastlen ); 
    //purify .................. stop
    //printf("S=%s\n", oneline.Data() );
   //........ parse oneline
    if (  (oneline.Index("#")==0) || 
	  (oneline.Index("@")==0) ||
	  (oneline.Index("END")==0) 
 ){ // starts with # - COMMENT HERE
      //       printf( "COMM:%s\n", oneline.Data() );
      if (title.Length()==0){ title=oneline( 1,oneline.Length()-1  );}
     //     title=oneline( oneline.Length(),oneline.Length()-1  );
  }else{// DATA HERE
      TObjArray *tar; 
      if (oneline.Length()>0){ // 1digit also exist...
       if (oneline.Index("#")>0){oneline.Remove(oneline.Index("#") );}// #
       if (oneline.Index("@")>0){oneline.Remove(oneline.Index("@") );}// #
       if (oneline.Index("END")>0){oneline.Remove(oneline.Index("END") );}// #
       //       oneline.ReplaceAll(oneline.Index("\t"),1," ");
       oneline.ReplaceAll("\t"," ");
      x[i]=0; y[i] =0;dx[i]=0;dy[i]=0; j=0;  //  go through the columns 
      // printf( "NotCM:%s\n", oneline.Data() );

      /*
       * r?? contains a column to use for x,dx,y,dy:
       */
      tar= oneline.Tokenize(" ");
      //      printf("entries==%d\n", tar->GetEntries() );
      while( (j<=rx && rx>=0) || (j<=ry && ry>=0) || (j<=rdy && rdy>=0)|| (j<=rdx && rdx>=0)){  
	if (j<tar->GetEntries()){
         token= ((TObjString*)(tar->At(j)))->GetString();
	 //	 printf(" %3d %3d  token <%s>\n", i,j,  token.Data() );
	 bu= token.Atof(); 
	 if (i<5){ 	 printf("%3d.  token <%s> = %f\n", i, token.Data(), bu ); }
         if (rx==j)  { x[i] =bu;}
         if (ry==j)  { y[i] =bu;}
         if (rdx==j) { dx[i] =bu;}
         if (rdy==j) { dy[i] =bu;} 
	}//j<entries
         j++;
      }//small while 

     //tar->Delete();  // abandon tar 
      i++;// skip to a next datum if point
     }// if oneline  lenght>1 go thru all
     //back the loop
   }//DATAH ERE end
    //     printf("%s","data line ended\n");


    }//if oneline.Gets false
   
   //........ parse oneline
  }//while
  //.......................readout HERE end....i
  //.......................close, create the TGraphErrors
  fclose( pFile);
  // TITLE HERE----------------
  //  title=name; DONE earlier
  title.ReplaceAll(".","_"); 
  title.ReplaceAll(" ","_"); 
  printf("i... %d elements read. title= /%s/\n", i-1, title.Data() );
 TGraphErrors *g=new TGraphErrors(i,x,y,dx,dy);
 g->SetMarkerStyle(22);
 g->SetTitle( title.Data()  );
//.........这里部分代码省略.........

//.........这里部分代码省略.........
		if (ev->edge(0,k)) continue;
		if (ev->range(0,k)<=0) continue;
		if (ev->range(0,k)>160) continue;
		if (ev->maxpixel(0,k)>250) continue;
		if (ev->x(0,k)<24) continue;
		if (ev->y(0,k)<24) continue;
		if (ev->x(0,k)>1000) continue;
		if (ev->y(0,k)>1000) continue;
		if (ev->maxpixel(0,k)/ev->E(0,k)>0.25) continue;
    	      	    
		xx = ev->x(0,k)-512;
        yy = ev->y(0,k)-512;      
        //if (yy >760-512&&yy<780-512) continue;//what is this?
        r2 = (xx*xx+yy*yy)* 16.66/1024*16.66/1024; 
        
        passTrig=0;
        passTrack++;
         
		if(passTrack==1)
		{
		E1 = ev->E(0,k)/lightCalib;
		if (useRadial) E1 *= pow(1+ r2/(16.123*16.123),3.5857/2 ); 
		Etrack = E1;
		phi = ev->phi(0,k);
        phi = atan2(sin(phi),cos(phi)) * 180./ 3.1416;   
	  	TrackX = xx;
		TrackY = yy;
		range_ccd=ev->range(0,k)*mmperpixel;
		eh1->Fill(Etrack);
		}
		else
		continue;
             
		TObjArray* arr = ev->waveform_vectors();
		FastWfVector* mesh = (FastWfVector*) arr->At(0);
		CspWfVector* anode = (CspWfVector*) arr->At(1);
		CspWfVector* veto = (CspWfVector*) arr->At(2);
	     	  
	    //loop over triggers
	    ntrig=mesh->size();
	    
		for (int nt = 0; nt < mesh->size(); nt++){	  
		diff=1e9;
		anodeRMS=mesh->at(nt).getRMS();
		meshRMS=anode->at(nt).getRMS();
	    if (mesh->at(nt).getRMS() >0.004){cutnum=1; continue;}//for neutron calib run: >0.002
	    if (anode->at(nt).getRMS()>0.0006){cutnum=2; continue;}//for neutron calib run: >0.003
	    if (fabs(mesh->at(nt,0).getFastRise0())>0.5){cutnum=3; continue;}
	    if (fabs(mesh->at(nt,0).getFastRise10())>0.5){cutnum=4; continue;}
	    if (fabs(mesh->at(nt,0).getFastRise25())>0.5){cutnum=5; continue;}
	    if (fabs(mesh->at(nt,0).getFastRise50())>0.5){cutnum=6; continue;}
	    if (fabs(mesh->at(nt,0).getFastRise75())>0.5){cutnum=7; continue;}  
	    if (fabs(mesh->at(nt,0).getFastRise90())>0.5){cutnum=8; continue;}
	    if (fabs(mesh->at(nt,0).getRise0())>0.5){cutnum=9; continue;}
	    if (fabs(mesh->at(nt,0).getRise10())>0.5){cutnum=10; continue;}
	    if (fabs(mesh->at(nt,0).getRise25())>0.5){cutnum=11; continue;}
	    if (fabs(mesh->at(nt,0).getRise50())>0.5){cutnum=12; continue;}
	    if (fabs(mesh->at(nt,0).getRise75())>0.5){cutnum=13; continue;}
	    if (fabs(mesh->at(nt,0).getRise90())>0.5){cutnum=14; continue;}
	    if (fabs(mesh->at(nt,0).getSlowFall0())>0.5){cutnum=15; continue;}
	    if (fabs(mesh->at(nt,0).getSlowFall10())>0.5){cutnum=16; continue;}
	    if (fabs(mesh->at(nt,0).getSlowFall25())>0.5){cutnum=17; continue;}
	    if (fabs(mesh->at(nt,0).getSlowFall50())>0.5){cutnum=18; continue;}
	    if (fabs(mesh->at(nt,0).getSlowFall75())>0.5){cutnum=19; continue;}
	    if (fabs(mesh->at(nt,0).getSlowFall90())>0.5){cutnum=20; continue;}
	    if (fabs(mesh->at(nt,0).getFall0())>0.5){cutnum=21; continue;}

void ContourList(){
 
   const Double_t PI = TMath::Pi(); 
    
   TCanvas* c = new TCanvas("c","Contour List",0,0,600,600);
   c->SetRightMargin(0.15);
   c->SetTopMargin(0.15);
 
   Int_t i, j, TotalConts;
    
   Int_t nZsamples   = 80;
   Int_t nPhiSamples = 80;

   Double_t HofZwavelength = 4.0;       // 4 meters 
   Double_t dZ             =  HofZwavelength/(Double_t)(nZsamples - 1);
   Double_t dPhi           = 2*PI/(Double_t)(nPhiSamples - 1);
 
   TArrayD z(nZsamples);
   TArrayD HofZ(nZsamples);
   TArrayD phi(nPhiSamples);
   TArrayD FofPhi(nPhiSamples);

   
   // Discretized Z and Phi Values
   for ( i = 0; i < nZsamples; i++) {
      z[i] = (i)*dZ - HofZwavelength/2.0;
      HofZ[i] = SawTooth(z[i], HofZwavelength);
   }

   for(Int_t i=0; i < nPhiSamples; i++){
      phi[i] = (i)*dPhi;      
      FofPhi[i] = sin(phi[i]); 
   }

   // Create Histogram
   TH2D *HistStreamFn = new TH2D("HstreamFn", 
   "#splitline{Histogram with negative and positive contents. Six contours are defined.}{It is plotted with options CONT LIST to retrieve the contours points in TGraphs}", 
   nZsamples, z[0], z[nZsamples-1], nPhiSamples, phi[0], phi[nPhiSamples-1]);
     
   // Load Histogram Data
   for (Int_t i = 0; i < nZsamples; i++) {
      for(Int_t j = 0; j < nPhiSamples; j++){
         HistStreamFn->SetBinContent(i,j, HofZ[i]*FofPhi[j]);
      } 
   }

   gStyle->SetPalette(1); 
   gStyle->SetOptStat(0);
   gStyle->SetTitleW(0.99);
   gStyle->SetTitleH(0.08);
   
   Double_t contours[6];
   contours[0] = -0.7;
   contours[1] = -0.5;
   contours[2] = -0.1;
   contours[3] =  0.1;
   contours[4] =  0.4;
   contours[5] =  0.8;
     
   HistStreamFn->SetContour(6, contours); 
 
   // Draw contours as filled regions, and Save points
   HistStreamFn->Draw("CONT Z LIST");
   c->Update(); // Needed to force the plotting and retrieve the contours in TGraphs
    
   // Get Contours
   TObjArray *conts = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
   TList* contLevel = NULL;
   TGraph* curv     = NULL;    
   Int_t nGraphs    = 0;
   Int_t TotalConts = 0;
    
   
   if (conts == NULL){
      printf("*** No Contours Were Extracted!\n");
      TotalConts = 0;
      return;
   } else {
      TotalConts = conts->GetSize();
   }

   printf("TotalConts = %d\n", TotalConts);
        
   for(i = 0; i < TotalConts; i++){       
      contLevel = (TList*)conts->At(i);
      printf("Contour %d has %d Graphs\n", i, contLevel->GetSize());
      nGraphs += contLevel->GetSize();
   }
    
   nGraphs = 0;
    
   TCanvas* c1 = new TCanvas("c1","Contour List",610,0,600,600);
   c1->SetTopMargin(0.15);
   TH2F *hr = new TH2F("hr",
   "#splitline{Negative contours are returned first (highest to lowest). Positive contours are returned from}{lowest to highest. On this plot Negative contours are drawn in red and positive contours in blue.}",
   2, -2, 2, 2, 0, 6.5);

   hr->Draw();
   Double_t x0, y0, z0;
   TLatex l;
//.........这里部分代码省略.........

void fitter(
            const string workDirName="Test", // Working directory
            // Select the type of datasets to fit
            bool fitData     = true,         // Fits Data if true, otherwise fits MC
            bool fitPbPb     = true,         // Fits PbPb datasets
            bool fitPP       = true,         // Fits PP datasets
            // Select the type of object to fit
            bool incJpsi     = true,         // Includes Jpsi model
            bool incPsi2S    = true,         // Includes Psi(2S) model
            bool incBkg      = true,         // Includes Background model
            // Select the fitting options
            bool cutCtau     = false,        // Apply prompt ctau cuts
            bool doSimulFit  = false,        // Do simultaneous fit
            bool wantPureSMC = false,        // Flag to indicate if we want to fit pure signal MC
            int  numCores    = 2,            // Number of cores used for fitting
            // Select the drawing options
            bool setLogScale = true,         // Draw plot with log scale
            bool incSS       = false,        // Include Same Sign data
            bool zoomPsi     = false,        // Zoom Psi(2S) peak on extra pad
            int  nBins       = 54            // Number of bins used for plotting
            ) 
{
  // -------------------------------------------------------------------------------
  // STEP 0: INITIALIZE THE FITTER WORK ENVIROMENT
  // The work enviroment is divided as follows:
  /*
    main |-> Macros: Contain all the macros
         |-> Input   |-> <WorkDir> : Contain Input File, Bin and Parameter List for a given work directory (e.g. 20160201)
	 |-> Output  |-> <WorkDir> : Contain Output Plots and Results for a given work directory (e.g. 20160201)
	 |-> DataSet : Contain all the datasets (MC and Data)
  */

  if (!checkSettings(fitData, fitPbPb, fitPP, incJpsi, incPsi2S, incBkg, cutCtau, doSimulFit, wantPureSMC, setLogScale, zoomPsi, incSS, numCores, nBins)) { return; }

  map<string,string> DIR;
  if(!iniWorkEnv(DIR, workDirName)){ return; }
 
  // -------------------------------------------------------------------------------
  // STEP 1: CREATE/LOAD THE ROODATASETS
  /*
    Input : List of TTrees with format:  TAG <tab> FILE_NAME
    Output: Collection of RooDataSets splitted by tag name, including OS and SS dimuons.
  */
  
  const string InputTrees = DIR["input"] + "InputTrees.txt";
  map<string, vector<string> > InputFileCollection;
  if(!getInputFileNames(InputTrees, InputFileCollection)){ return; }
  
  TObjArray* aDSTAG = new TObjArray(); // Array to store the different tags in the list of trees
  aDSTAG->SetOwner(true);
  map<string, RooWorkspace> Workspace;

  for(map<string, vector<string> >::iterator FileCollection=InputFileCollection.begin(); FileCollection!=InputFileCollection.end(); ++FileCollection) {
    // Get the file tag which has the following format: DSTAG_COLL , i.e. DATA_PP 
    string FILETAG = FileCollection->first;  
    string DSTAG   = FILETAG;
    if (FILETAG.size()) {
      DSTAG.erase(DSTAG.find("_"));
    } else {
      cout << "[ERROR] FILETAG is empty!" << endl;
    }
    // Extract the filenames
    vector<string> InputFileNames = FileCollection->second; 
    string         OutputFileName;
    // If we have data, check if the user wants to fit data
    if ( (FILETAG.find("DATA")!=std::string::npos) && fitData==true ) {
      if ( (FILETAG.find("PP")!=std::string::npos)   && !fitPP   ) continue; // If we find PP, check if the user wants PP
      if ( (FILETAG.find("PbPb")!=std::string::npos) && !fitPbPb ) continue; // If we find PbPb, check if the user wants PbPb
      OutputFileName = DIR["dataset"] + "DATASET_" + FILETAG + ".root";
      if(!tree2DataSet(Workspace[DSTAG], InputFileNames, FILETAG, OutputFileName)){ return; }
      if (!aDSTAG->FindObject(DSTAG.c_str())) aDSTAG->Add(new TObjString(DSTAG.c_str()));
    }
    // If we find MC, check if the user wants to fit MC
    if ( (FILETAG.find("MC")!=std::string::npos) && fitData==false ) {
      if ( (FILETAG.find("PP")!=std::string::npos)    && !fitPP    ) continue; // If we find PP, check if the user wants PP
      if ( (FILETAG.find("PbPb")!=std::string::npos)  && !fitPbPb  ) continue; // If we find PbPb, check if the user wants PbPb
      if ( (FILETAG.find("JPSI")!=std::string::npos)  && !incJpsi  ) continue; // If we find Jpsi MC, check if the user wants to include Jpsi
      if ( (FILETAG.find("PSI2S")!=std::string::npos) && !incPsi2S ) continue; // If we find Psi2S MC, check if the user wants to include Psi2S
      OutputFileName = DIR["dataset"] + "DATASET_" + FILETAG + ".root";
      if(!tree2DataSet(Workspace[DSTAG], InputFileNames, FILETAG, OutputFileName)){ return; }
      if (!aDSTAG->FindObject(DSTAG.c_str())) aDSTAG->Add(new TObjString(DSTAG.c_str()));
      if (wantPureSMC)
      {
        OutputFileName = DIR["dataset"] + "DATASET_" + FILETAG + "_PureS" + ".root";
        if(!tree2DataSet(Workspace[Form("%s_PureS",DSTAG.c_str())], InputFileNames, FILETAG, OutputFileName)){ return; }
      }
    } 
  }
  if (Workspace.size()==0) {
    cout << "[ERROR] No onia tree files were found matching the user's input settings!" << endl; return;
  }

  // -------------------------------------------------------------------------------
  // STEP 2: LOAD THE INITIAL PARAMETERS
  /*
    Input : List of initial parameters with format PT <tab> RAP <tab> CEN <tab> iniPar ... 
    Output: two vectors with one entry per kinematic bin filled with the cuts and initial parameters
  */
  
  string InputFile;
//.........这里部分代码省略.........

//_____________________________
void muIDCutsOptim::Loop()
{
  //  In a ROOT session, you can do:
  //      root> .L muIDCutsOptim.C
  //      root> muIDCutsOptim t
  //      root> t.GetEntry(12); // Fill t data members with entry number 12
  //      root> t.Show();       // Show values of entry 12
  //      root> t.Show(16);     // Read and show values of entry 16
  //      root> t.Loop();       // Loop on all entries
  //
  
  //  This is the loop skeleton where:
  //    jentry is the global entry number in the chain
  //    ientry is the entry number in the current Tree
  //  Note that the argument to GetEntry must be:
  //    jentry for TChain::GetEntry
  //    ientry for TTree::GetEntry and TBranch::GetEntry
  //

  if (fChain == 0) return;
  
  
  //===== Definition of file type, particle and if we want the high purity cut included in the soft muon ID cuts
  const char* fileType = "MC";
  const char* particle = "JPsi";
  fIspp = kFALSE;
  fIncludeHighPurity = kFALSE;
  //=====
  
  
  Long64_t nentries(0);
  if ( !strcmp(fileType,"MC") ) nentries = fChain->GetEntries();
  else nentries = 1500000;
  
  std::cout << "# Events = " << nentries << std::endl;
  
  
  //===== Definition of some cuts
  Double_t leMinvSig = 3.0; // Signal range
  Double_t ueMinvSig = 3.2;
  
  Double_t leMinvBkg1 = 2.7; // Bkg 1 range (sideband)
  Double_t ueMinvBkg1 = 2.9;
  
  Double_t leMinvBkg2 = 3.3; // Bkg 1 range (sideband)
  Double_t ueMinvBkg2 = 3.5;

  Double_t Ptmin = 3.0; // Pt cut (6.5) (3.0)
  Double_t Ptmax = 6.5; // (12.0) (6.5)
  
  Double_t Ymin = 2.0; // Y cut (0.0) (2.0)
  Double_t Ymax = 2.4; // (2.4) (2.4)
  //=====
  
  
  //===== Definition of MuID variables and binning
  const int nvar = 19;
  const char* varname[19] = {"isGoodMuon", "highPurity", "TrkMuArb", "TMOneStaTight", "nPixValHits",
    "nMuValHits", "nTrkHits", "normChi2_inner", "normChi2_global", "nPixWMea",
    "nTrkWMea", "StationsMatched", "dxy", "dxyErr", "dz",
    "dzErr", "ptErr_inner", "ptErr_global","VtxProb"};
  const int nbins[19] = {2, 2, 2, 2, 11,
    56, 36, 100, 100, 6,
    19, 7, 100, 100, 100,
    100, 100, 100, 200};
  const double lowedge[19] = {0, 0, 0, 0, 0,
    0, 0, 0, 0, 0,
    0, 0, 0, 0, 0,
    0, 0, 0, 0};
  const double hiedge[19] = {2, 2, 2, 2, 11,
    56, 36, 10, 50, 6,
    19, 7, 2, 0.25, 20,
    5, 0.15, 0.15, 1.};
  //=====
  
  
  SelectBranches(nvar,varname); // Activates only desired branches (General and MuID ones in this case)
  
  
  //====== Create histograms and arrays to store them ========
  TObjArray* aDistr = new TObjArray();
  TObjArray* aSig = new TObjArray();
  TObjArray* aBkg = new TObjArray();
  TObjArray* aSigC = new TObjArray();
  TObjArray* aBkgC = new TObjArray();
  TObjArray* aGlob = new TObjArray();
  TObjArray* aMinvCuts = new TObjArray();
  TObjArray* aSigSoftMuCutsButOne = new TObjArray();
  TObjArray* aBkgSoftMuCutsButOne = new TObjArray();
  TObjArray* aSigSoftMuCuts = new TObjArray();
  TObjArray* aBkgSoftMuCuts = new TObjArray();
  
  map<TString, TH1F*> hists_sig, hists_bkg, hists_sig_distr, hists_bkg_distr, hists_distr, hists_sig_distr_CutTest, hists_bkg_distr_CutTest, hists_sig_SofMu, hists_bkg_SofMu;
  
  TString sHPName("");
  if ( fIncludeHighPurity ) sHPName += "highPurity"; // This is for the histograms creation
  
  TObjArray* dummyA(0x0);
  TString partType("");
//.........这里部分代码省略.........

//
// *** Configuration script for phi->KK analysis with 2010 runs ***
// 
// A configuration script for RSN package needs to define the followings:
//
// (1) decay tree of each resonance to be studied, which is needed to select
//     true pairs and to assign the right mass to all candidate daughters
// (2) cuts at all levels: single daughters, tracks, events
// (3) output objects: histograms or trees
//
Bool_t RsnConfigPhiTPC
(
   AliRsnAnalysisTask *task,
   Bool_t              isMC,
   Bool_t              isMix,
   Bool_t              useCentrality,
   AliRsnCutSet       *eventCuts
)
{
   if (!task) ::Error("RsnConfigPhiTPC", "NULL task");
   
   // we define here a suffix to differentiate names of different setups for the same resonance
   // and we define also the name of the list of tracks we want to select for the analysis
   // (if will fail if no lists with this name were added to the RsnInputHandler)
   const char *suffix     = "tpc";
   const char *listName   = "kaonTPC";
   Bool_t      useCharged =  kTRUE;
   Int_t       listID     = -1;
   
   // find the index of the corresponding list in the RsnInputHandler
   AliAnalysisManager        *mgr   = AliAnalysisManager::GetAnalysisManager();
   AliMultiInputEventHandler *multi = dynamic_cast<AliMultiInputEventHandler*>(mgr->GetInputEventHandler());
   if (multi) {
      TObjArray          *array = multi->InputEventHandlers();
      AliRsnInputHandler *rsn   = (AliRsnInputHandler*)array->FindObject("rsnInputHandler");
      if (rsn) {
         AliRsnDaughterSelector *sel = rsn->GetSelector();
         listID = sel->GetID(listName, useCharged);
      }
   }
   if (listID >= 0)
      ::Info ("RsnConfigPhiTPC.C", "Required list '%s' stays in position %d", listName, listID);
   else {
      ::Error("RsnConfigPhiTPC.C", "Required list '%s' absent in handler!", listName);
      return kFALSE;
   }
   
   // ----------------------------------------------------------------------------------------------
   // -- DEFINITIONS -------------------------------------------------------------------------------
   // ----------------------------------------------------------------------------------------------
   
   // PAIR DEFINITIONS:
   // this contains the definition of particle species and charge for both daughters of a resonance,
   // which are used for the following purposes:
   // --> species is used to assign the mass to the daughter (e.g. for building invariant mass)
   // --> charge is used to select what tracks to use when doing the computation loops
   // When a user wants to compute a like-sign background, he must define also a pair definition
   // for each like-sign: in case of charged track decays, we need one for ++ and one for --
   // Last two arguments are necessary only in some cases (but it is not bad to well initialize them):
   // --> PDG code of resonance, which is used for selecting true pairs, when needed
   // --> nominal resonance mass, which is used for computing quantities like Y or Mt
   AliRsnPairDef *phi_kaonP_kaonM = new AliRsnPairDef(AliRsnDaughter::kKaon, '+', AliRsnDaughter::kKaon, '-', 333, 1.019455);
   AliRsnPairDef *phi_kaonP_kaonP = new AliRsnPairDef(AliRsnDaughter::kKaon, '+', AliRsnDaughter::kKaon, '+', 333, 1.019455);
   AliRsnPairDef *phi_kaonM_kaonM = new AliRsnPairDef(AliRsnDaughter::kKaon, '-', AliRsnDaughter::kKaon, '-', 333, 1.019455);

   // PAIR LOOPS:
   // these are the objects which drive the computations and fill the output histograms
   // each one requires to be initialized with an AliRsnPairDef object, which provided masses,
   // last argument tells if the pair is for mixing or not (this can be also set afterwards, anyway)
   const Int_t     nPairs = 5;
   Bool_t          addPair[nPairs] = {1, 1, 1, 1, 1};
   AliRsnLoopPair *phiLoop[nPairs];
   phiLoop[0] = new AliRsnLoopPair(Form("%s_phi_kaonP_kaonM"     , suffix), phi_kaonP_kaonM, kFALSE);
   phiLoop[1] = new AliRsnLoopPair(Form("%s_phi_kaonP_kaonM_true", suffix), phi_kaonP_kaonM, kFALSE);
   phiLoop[2] = new AliRsnLoopPair(Form("%s_phi_kaonP_kaonM_mix" , suffix), phi_kaonP_kaonM, kTRUE );
   phiLoop[3] = new AliRsnLoopPair(Form("%s_phi_kaonP_kaonP"     , suffix), phi_kaonP_kaonP, kFALSE);
   phiLoop[4] = new AliRsnLoopPair(Form("%s_phi_kaonM_kaonM"     , suffix), phi_kaonM_kaonM, kFALSE);

   // set additional option for true pairs
   // 1) we select only pairs coming from the same mother, which must have the right PDG code (from pairDef)
   // 2) we select only pairs decaying according to the right channel (from pairDef species+charge definitions)
   phiLoop[1]->SetOnlyTrue(kTRUE);
   phiLoop[1]->SetCheckDecay(kTRUE);
   
   // don't add true pairs if not MC
   if (!isMC) addPair[1] = 0;
   addPair[0] = !isMix;
   addPair[1] = !isMix;
   addPair[2] =  isMix;
   addPair[3] = !isMix;
   addPair[4] = !isMix;
   
   // ----------------------------------------------------------------------------------------------
   // -- COMPUTED VALUES & OUTPUTS -----------------------------------------------------------------
   // ----------------------------------------------------------------------------------------------
   
   // All values which should be computed are defined here and passed to the computation objects,
   // since they define all that is computed bye each one, and, in case one output is a histogram
   // they define the binning and range for that value
   //
//.........这里部分代码省略.........

void browseStacks( bool makePictures=false, bool wait=true , bool addHistName = false, Double_t maxYScaleF = 1., 
                  bool logScale = false, bool setMinZero = true) {


  gStyle->SetOptTitle(0);

  bool keep2D=false;

  //fix the hNJet histos
  TList *list = gDirectory->GetList();
  TIterator *iter = list->MakeIterator();
  TObject *obj = 0;
  while(obj = iter->Next()) {
  
    if(TString(obj->GetName()).Contains("hnJet") && obj->InheritsFrom(TH1::Class())) {
      
      int nbins = ((TH1F*)obj)->GetNbinsX();
      float overflow = ((TH1F*)obj)->GetBinContent(nbins+1);
      float lastbinval = ((TH1F*)obj)->GetBinContent(nbins);
      ((TH1F*)obj)->SetBinContent(nbins, overflow+lastbinval);
      ((TH1F*)obj)->GetXaxis()->SetBinLabel(nbins, "#geq4");
    }
  }
  
    
    
    
  // Find out what the names of the existing histograms are
  // The histogram names are XX_YY_ZZ, where XX is the sample,
  // eg, "tt", YY is the actual name, ZZ is the final state, eg, "ee"
  TObjArray* myNames = getMyHistosNames("ttdil","ee",keep2D);
    

  // Now loop over histograms, and make stacks
  TCanvas *c = new TCanvas();
  c->Divide(2,2);
  char* suffix[4];
  suffix[0] = "ee";
  suffix[1] = "mm";
  suffix[2] = "em";
  suffix[3] = "all";
  if (makePictures) c->Print("out/stacks.ps[");
  for (int i=0; i<myNames->GetEntries(); i++) {
     
    for (int sample=0; sample<4; sample++) {
       
       
      hist::stack(Form("st_%s_%s",myNames->At(i)->GetName(),suffix[sample]),
		  Form("%s_%s$",myNames->At(i)->GetName(), suffix[sample]));
      THStack* thisStack = (THStack*) gROOT->FindObjectAny(
							   Form("st_%s_%s", myNames->At(i)->GetName(), suffix[sample]));
       
      thisStack->SetMaximum(thisStack->GetMaximum()*maxYScaleF);
      if(TString(myNames->At(i)->GetName()).Contains("hnJet")) {
	TList* histolist = thisStack->GetHists();
	int hatchcount = 0;
	// 	for(int j = 0; j<histolist->GetSize();j++) {
	// 	  if(TString(histolist->At(j)->GetName()).Contains("tt") ||
	// 	     TString(histolist->At(j)->GetName()).Contains("tautau") ||
	// 	     TString(histolist->At(j)->GetName()).Contains("ww") ) continue;
	// 	  hatch(histolist->At(j)->GetName(), FavoriteHatches[hatchcount]);
	// 	  hatchcount++;
	// 	}
      }
	 
	 
      TLegend* thisLeg = hist::legend(thisStack, "lpf", 0, 0, 0.75, 0.65, 0.99, 0.99);
      c->cd(sample+1);
      if (logScale) gPad->SetLogy(); else gPad->SetLogy(0);
      double stackMax = ((TH1*)thisStack->GetHists()->At(0))->GetMaximum();
      double stackMin = ((TH1*)thisStack->GetHists()->At(0))->GetMinimum();
      thisStack->SetMinimum(stackMin);
      if (setMinZero) thisStack->SetMinimum(0);
      if (logScale && stackMin <=0) thisStack->SetMinimum(1e-2*stackMax);
      if (logScale && stackMax == 0) thisStack->SetMinimum(1e-12); 
      thisStack->Draw("hist");
      string xtitle( ((TH1*)gROOT->FindObjectAny(Form("ttdil_%s_%s", myNames->At(i)->GetName(), suffix[sample])))->GetXaxis()->GetTitle());
      string ytitle( ((TH1*)gROOT->FindObjectAny(Form("ttdil_%s_%s", myNames->At(i)->GetName(), suffix[sample])))->GetYaxis()->GetTitle());
      thisStack->GetXaxis()->SetTitle(xtitle.c_str());
      thisStack->GetYaxis()->SetTitle(ytitle.c_str());
      TString hname = thisStack->GetName();
      if(hname.Contains("hnJet")) {
	thisStack->GetXaxis()->SetLabelSize(0.075);
	thisStack->GetYaxis()->SetLabelSize(0.05);
	thisStack->GetXaxis()->SetTitle("N_{jets}");
      }
      thisLeg->Draw();
	
      TPaveText *pt1 = new TPaveText(0.1, 0.95, 0.4, 0.999, "brNDC");
      pt1->SetName("pt1name");
      pt1->SetBorderSize(0);
      pt1->SetFillStyle(0);
	
      TText *blah;
      if (addHistName) blah = pt1->AddText(hname);
      else blah = pt1->AddText("CMS Preliminary");
      blah->SetTextSize(0.05);
      pt1->Draw();
      c->Modified();
	
//.........这里部分代码省略.........

void GeomAcc::Calculate(TString geomTuple, TString saveName){

  TFile* file = new TFile(geomTuple);
  TTree* tree=(TTree*)file->Get("SelectJpsi/JpsiGen"); 

  double JpsiPT, JpsiRap, CosTheta;
  double MuNegDau, MuPosDau;

  tree->SetBranchAddress("JPsi_PT",&JpsiPT);
  tree->SetBranchAddress("JPsi_RAPIDITY",&JpsiRap);
  tree->SetBranchAddress("JPsi_COSTHETA",&CosTheta);
  tree->SetBranchAddress("MuMinus_THETA",&MuNegDau);
  tree->SetBranchAddress("MuPlus_THETA",&MuPosDau);


  int xbin = 20;
  int ybin = 5;

  TH2D *Num = new TH2D("Num","Num",xbin,0.0,10000.0, ybin,2.0,4.2);
  TH2D *Den = new TH2D("Den","Den",xbin,0.0,10000.0, ybin,2.0,4.2);
  
  TH2D *NumLong = new TH2D("NumLong","NumLong",xbin,0.0,10000.0, ybin,2.0,4.2);
  TH2D *DenLong = new TH2D("DenLong","DenLong",xbin,0.0,10000.0, ybin,2.0,4.2);
 
  TH2D *NumTrans = new TH2D("NumTrans","NumTrans",xbin,0.0,10000.0, ybin,2.0,4.2);
  TH2D *DenTrans = new TH2D("DenTrans","DenTrans",xbin,0.0,10000.0, ybin,2.0,4.2);

  for(int i=0; i<tree->GetEntries(); i++)
  {
    tree->GetEntry(i);

    Den->Fill(JpsiPT,JpsiRap);
    DenLong->Fill(JpsiPT,JpsiRap,(1-CosTheta*CosTheta));
    DenTrans->Fill(JpsiPT,JpsiRap,(1+CosTheta*CosTheta));

    //
    //Calculate Daughters in acceptance criteria
    //
    bool inAcc = Cut(MuPosDau, MuNegDau);

    if(inAcc) {
      Num->Fill(JpsiPT,JpsiRap);
      NumLong->Fill(JpsiPT,JpsiRap,(1-CosTheta*CosTheta));
      NumTrans->Fill(JpsiPT,JpsiRap,(1+CosTheta*CosTheta));
    }
  }

  //
  //Calculate efficiency
  //

  TH2D GeoEff =  (Eff(Num, Den, xbin, ybin, "GeoEff"));
  TH2D GeoLongEff =  (Eff(NumLong, DenLong, xbin, ybin,"GeoLongEff"));
  TH2D GeoTransEff =  (Eff(NumTrans, DenTrans, xbin, ybin,"GeoTransEff"));


  TObjArray* histList = new TObjArray(0);
  histList->Add(Num);
  histList->Add(Den);
  histList->Add(&GeoEff);
  histList->Add(NumLong);
  histList->Add(DenLong);
  histList->Add(&GeoLongEff);
  histList->Add(NumTrans);
  histList->Add(DenTrans);
  histList->Add(&GeoTransEff);
  TFile* saveFile = new TFile("Output/"+saveName, "recreate");
  histList->Write();
  saveFile->Close();

  //
  // Save into final usable format in Efficiency folder
  //

  Format(saveName, "NoPol/");
  Format(saveName, "Long/");
  Format(saveName, "Trans/");

  return;
}

int main(int argc, char* argv[]){
	if (argc != 6){
    		std::cout << "wrong number of arguments: usage ./topDataClass_doAnalysis <file name> <mod #> <row #> <col #> <ch #>" << std::endl;
    		return 0;
  	}

	//define application object
	theApp = new TApplication("App", &argc, argv);
	TString inputFileName = theApp->Argv()[1];
	std::cout << "Input file name "  << inputFileName << std::endl;

	//create target6 interface object
	topDataClass *data = new topDataClass();

	//specify channel of interest
	int inMod = atoi(theApp->Argv()[2]);
	int inRow = atoi(theApp->Argv()[3]);
	int inCol = atoi(theApp->Argv()[4]);
	int inCh = atoi(theApp->Argv()[5]);
	data->setAnalysisChannel( inMod, inRow, inCol, inCh );

	//specify timing marker channel
	data->setTimingMarkerChannel( 0, 0, 1, 2 );

	//specify time window
	data->windowTime = windowTime;

	//open summary tree file
	data->openSummaryTree(inputFileName);

	//create output file
  	TObjArray* strings = inputFileName.Tokenize("/");
  	TObjString* objstring = (TObjString*) strings->At(strings->GetLast());
  	TString inputFileNameBase(objstring->GetString());
	TString outputFileName = "output_topDataClass_doAnalysis_doublePulseSampleDTFit_";
  	outputFileName += inMod;
  	outputFileName += "_";
  	outputFileName += inRow;
  	outputFileName += "_";
  	outputFileName += inCol;
  	outputFileName += "_";
  	outputFileName += inCh;
  	outputFileName += "_";
  	outputFileName += inputFileNameBase;
  	//outputFileName += ".root";
  	std::cout << " outputFileName " << outputFileName << std::endl;
  	outputFile = new TFile( outputFileName , "RECREATE");

	//initialize histograms
	initializeGlobalHistograms();

	//initialize tree branches
	data->setTreeBranches();

	//load pulse info into arrays
	data->selectPulsesForArray();

	//loop over selected events, apply corrections histogram pulse time distributions
	//monitor pulse time vs event #
  	gPulseTimeVsEventNum = new TGraphErrors();
  	for(int entry = 0; entry < data->numUsed; entry++) {
		//skip events not in arrays
  		if( entry >= maxNumEvt )
			continue;

		double pulseTime = data->measurePulseTimeArrayEntry(entry,1);
		double pulseHeight = data->adc_0_A[entry];
		int smpBinNumIn128Array = data->getSmpBinNumIn128Array(entry);
		double smpPos =	data->getSmpPos(entry);

		//apply analysis cuts
		//if( smpBinNumIn128Array >= 127 )
		//	continue;
		//measure pulse time vs event #
		gPulseTimeVsEventNum->SetPoint( gPulseTimeVsEventNum->GetN() , data->eventNum_A[entry],  pulseTime );
		//cut on event #
		if( data->eventNum_A[entry] < 0 )
			continue;

		//histogram selected pulse distributions
		hPulseHeightInitial->Fill( data->adc_0_A[entry] );
		hPulseTimeInitial->Fill( pulseTime );
		hPulseSampleInitial->Fill(data->smp_0_A[entry]);
       		hPulseSmp128Initial->Fill(smpBinNumIn128Array);
		hPulseTimeVsSmp128Initial->Fill(smpBinNumIn128Array, pulseTime );
		hPulseTimeVsSmp128PosInitial->Fill(smpBinNumIn128Array + smpPos, pulseTime );
		hPulseTimeVsHeightInitial->Fill(pulseHeight, pulseTime );
		hPulseTimeVsFTSWInitial->Fill(data->ftsw_A[entry], pulseTime );
   		hFTSWVsSmp128Initial->Fill(smpBinNumIn128Array, data->ftsw_A[entry]);
	}
	
	//loop over selected events, measure time difference between timing marker and selected channel pulses
	gTest = new TGraphErrors();
	data->measurePulseMarkerTimeDifferenceDistribution(hPulseTimeMarkTimeDiffInitial,hPulseTimeMarkTimeDiffVsMarkSmpBinNumInitial );

	//run fit analysis
	doDoublePulseFit(data);

	//write output file
	writeOutputFile();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  }  
  hN= new THnF("exampleFit","exampleFit", ndim, nbins, xmin,xmax);
  //
  // 1.) generate random input points
  //
  for (Int_t ipoint=0; ipoint<TMath::Abs(npoints); ipoint++){
    for (Int_t idim=0; idim<ndim; idim++){
      xyz[idim]=gRandom->Rndm();
    }
    Double_t value=pformula->EvalPar(xyz,0);
    Double_t noise = gRandom->Gaus()*err;
    Double_t noise2 = noise*(1+(gRandom->Rndm()<0.1)*100);  // noise with 10 percent of outliers
    Double_t noiseBreit = gRandom->BreitWigner()*err;
    (*pcstreamIn)<<"testInput"<<
      "val="<<value<<
      "err="<<err<<
      "noise="<<noise<<            // gausian noise 
      "noise2="<<noise2<<          // gausian noise + 10% of outliers
      "noiseBreit="<<noiseBreit;      
    for (Int_t idim=0; idim<ndim; idim++){
      (*pcstreamIn)<<"testInput"<<chxyz[idim]->Data()<<xyz[idim];
    }
    (*pcstreamIn)<<"testInput"<<"\n";
  }
  delete pcstreamIn;
  pcstreamIn  = new TTreeSRedirector("fitNDLocalTestInput.root");
  treeIn= (TTree*)(pcstreamIn->GetFile()->Get("testInput"));
  //   treeIn->Draw("val:xyz0:xyz1>>his(20,0,1,20,0,1)","","profsurf2",10000); // visualization of input
  //
  // 2.) Make local fits 
  //
  if (npoints>0){
    pfitNDIdeal->SetStreamer(pcstreamOutIdeal);
    pfitNDGaus0->SetStreamer(pcstreamOutGaus0);
    pfitNDGaus1->SetStreamer(pcstreamOutGaus1);
    pfitNDGaus2->SetStreamer(pcstreamOutGaus2);
    pfitNDGaus3->SetStreamer(pcstreamOutGaus3);
    pfitNDBreit0->SetStreamer(pcstreamOutBreit0);
    pfitNDBreit1->SetStreamer(pcstreamOutBreit1);
    pfitNDBreit2->SetStreamer(pcstreamOutBreit2);
  }
  //
  pfitNDIdeal->SetHistogram((THn*)(hN->Clone()));
  pfitNDGaus0->SetHistogram((THn*)(hN->Clone()));
  pfitNDGaus1->SetHistogram((THn*)(hN->Clone()));
  pfitNDGaus2->SetHistogram((THn*)(hN->Clone()));
  pfitNDGaus3->SetHistogram((THn*)(hN->Clone()));
  pfitNDGaus2->SetCuts(3,0.8,1);
  pfitNDGaus3->SetCuts(0,0.8,0);

  pfitNDBreit0->SetCuts(3,0.8,1);
  pfitNDBreit1->SetCuts(3,0.8,1);
  pfitNDBreit2->SetCuts(0,0.8,1);
  pfitNDBreit0->SetHistogram((THn*)(hN->Clone()));
  pfitNDBreit1->SetHistogram((THn*)(hN->Clone()));
  pfitNDBreit2->SetHistogram((THn*)(hN->Clone()));
  //
  pfitNDIdeal->MakeFit(treeIn, "val:err", "xyz0:xyz1","Entry$%2==1", "0.02:0.02","2:2",0.0001);
  pfitNDGaus0->MakeFit(treeIn, "val+noise:err", "xyz0:xyz1","Entry$%2==0", "0.02:0.02","2:2",0.0001);  // sample Gaussian1
  pfitNDGaus1->MakeFit(treeIn, "val+noise:err", "xyz0:xyz1","Entry$%2==1", "0.02:0.02","2:2",0.0001);  // sample Gaussian2
  pfitNDGaus2->MakeFit(treeIn, "val+noise2:err", "xyz0:xyz1","Entry$%2==1", "0.02:0.02","2:2",0.0001);  // sample Gaussian2 - with tails robust
  pfitNDGaus3->MakeFit(treeIn, "val+noise2:err", "xyz0:xyz1","Entry$%2==1", "0.02:0.02","2:2",0.0001);  // sample Gaussian2 - with tails non robust
  pfitNDBreit0->MakeFit(treeIn, "val+noiseBreit:err", "xyz0:xyz1","Entry$%2==0", "0.02:0.02","2:2",0.0001);  // sample Breit0
  pfitNDBreit1->MakeFit(treeIn, "val+noiseBreit:err", "xyz0:xyz1","Entry$%2==1", "0.02:0.02","2:2",0.0001);  // sample Breit1
  pfitNDBreit2->MakeFit(treeIn, "val+noiseBreit:err", "xyz0:xyz1","Entry$%2==1", "0.02:0.02","2:2",0.0001);  // sample Breit2 without outlier filtering
  //
  if (npoints<0) return;  // callgrind mode of operation
  pfitNDIdeal->AddVisualCorrection(pfitNDIdeal,1);
  pfitNDGaus0->AddVisualCorrection(pfitNDGaus0,2);
  pfitNDGaus1->AddVisualCorrection(pfitNDGaus1,3);
  pfitNDGaus2->AddVisualCorrection(pfitNDGaus2,7);
  pfitNDGaus3->AddVisualCorrection(pfitNDGaus3,8);
  pfitNDBreit0->AddVisualCorrection(pfitNDBreit0,4);
  pfitNDBreit1->AddVisualCorrection(pfitNDBreit1,5);
  pfitNDBreit2->AddVisualCorrection(pfitNDBreit2,6);
  //

  TObjArray * array = AliNDLocalRegression::GetVisualCorrections();
  for (Int_t i=0; i<array->GetEntries(); i++){
    AliNDLocalRegression * regression = ( AliNDLocalRegression *)array->At(i);
    if (regression==NULL) continue;
    regression->AddVisualCorrection(regression);
    Int_t hashIndex = regression->GetVisualCorrectionIndex();
    treeIn->SetAlias( regression->GetName(), TString::Format("AliNDLocalRegression::GetCorrND(%d,xyz0,xyz1+0)",hashIndex).Data());
  }

  
  delete pcstreamOutGaus0;
  delete pcstreamOutGaus1;
  delete pcstreamOutGaus2;
  delete pcstreamOutGaus3;
  delete pcstreamOutBreit0;
  delete pcstreamOutBreit1;
  delete pcstreamOutBreit2;
  UnitTestGaussNoise(); 
  UnitTestGaussNoisePlusOutliers();
  UnitTestBreitWignerNoise();
  UnitTestStreamer();
  UnitTestContrain();
}

void PlotPubHisto(TObjArray histograms,TEnv *params){
    // This is a modification of the AddHistos macro
    
    // Number of histos to plot:
    Int_t ntot = histograms.GetEntries();
    
    // Check we have what we expect (the order should be: data, qcd, wjets, etc...)
    for(Int_t i = 0; i<ntot; i++){
	if(histograms[i]==0) {
	    cout<<"Error in AddHistos: histogram "<<i<<" is a NULL pointer!"<<endl;
	    return;
	}
	TH1F * hthis = (TH1F*) histograms[i];
	// include the overflow/underflow bins:
	int numbins = hthis->GetNbinsX(); //this is the last bin plotted
	double hicontent = hthis->GetBinContent(numbins);
	double overflow  = hthis->GetBinContent(numbins+1);// this bin con