C++ TTree::SetEntryList方法代码示例

double NewChi2Func(const double *xx ){
 const Double_t scale = xx[0];
 fileInMC->cd();
 TTree* MyTreeMC = (TTree*) fileInMC->Get(treeNameMC.c_str());
 outFile->cd();
 TString nameDATA = Form("hDATA_%d_%d_%.5f",Data_or_MC,nIter,ScaleTrue);
 TH1F* hDATA = (TH1F*) outFile->Get(nameDATA);
 
 TString NameMC = Form("hMC_Chi2_%.5f",scale);
 TH1F* hMC;
 if (!gROOT->FindObject(NameMC.Data())){
  hMC = new TH1F(NameMC,NameMC,numBINS,minBINS,maxBINS);
  hMC->Reset();
  
  MyTreeMC->SetEntryList(0); 
  MyTreeMC->Draw(">> myListMCTot",(AdditionalCut + Form(" && (ET * (1+(%f)))>%f",scale,minET)).Data(),"entrylist");
  TEntryList *mylistMCTot = (TEntryList*)gDirectory->Get("myListMCTot");
  MyTreeMC->SetEntryList(mylistMCTot);
  
  TString DrawMC = Form("(%s * (1+(%f)))>>%s",variableName.c_str(),scale,NameMC.Data());
  MyTreeMC->Draw(DrawMC);
//   MyTreeMC->Draw(DrawMC,(AdditionalCut+Form("&& (ET * (1+(%f)))>%f",scale,minET)).Data());
  hMC->Sumw2();
  hMC->Scale(hDATA->GetEffectiveEntries()/hMC->GetEffectiveEntries());
  outFile->cd();
  hMC->Write();
 } 
 else {
//   std::cerr << " NewChi2 old " << NameMC.Data() << std::endl;
  hMC = (TH1F*) outFile->Get(NameMC.Data());
 }
 
 double result = myChi2(hMC,hDATA);
 return result;
}

double Chi2F(const double *xx ){
 const Double_t scale = xx[0];
 fileInMC->cd();
 TTree* MyTreeMC = (TTree*) fileInMC->Get(treeNameMC.c_str());
 outFile->cd();
 TString nameDATA = Form("hDATA_%d_%d_%.5f",Data_or_MC,nIter,ScaleTrue);
 TH1F* hDATA = (TH1F*) outFile->Get(nameDATA);
 TString NameMC = Form("hMC_Chi2_%.5f",scale);
 
 //std::cerr << " NameMC = "  << NameMC.Data() << " => " << scale << std::endl;
 
 TH1F* hMC;
 if (!gROOT->FindObject(NameMC.Data())){
  hMC = new TH1F(NameMC,NameMC,numBINS,minBINS,maxBINS);
  hMC->Reset();

//  std::cerr << " AdditionalCut.Data() = " << AdditionalCut.Data() << std::endl;
  
  MyTreeMC->SetEntryList(0); 
  MyTreeMC->Draw(">> myListMCTot",(AdditionalCut + Form(" && (ET * (1+(%f)))>%f",scale,minET)).Data(),"entrylist");
  TEntryList *mylistMCTot = (TEntryList*)gDirectory->Get("myListMCTot");
  MyTreeMC->SetEntryList(mylistMCTot);
  
  TString DrawMC = Form("(%s * (1+(%f)))>>%s",variableName.c_str(),scale,NameMC.Data());
//  std::cerr << " DrawMC = " << DrawMC.Data() << std::endl;
  MyTreeMC->Draw(DrawMC);
//   MyTreeMC->Draw(DrawMC,(AdditionalCut+Form("&& (ET * (1+(%f)))>%f",scale,minET)).Data());
  hMC->Sumw2();
  hMC->Scale(hDATA->GetEffectiveEntries()/hMC->GetEffectiveEntries());
  outFile->cd();
  hMC->Write();
 }
 else {
//   std::cerr << " KM old "  << NameMC.Data() << std::endl;
  hMC = (TH1F*) outFile->Get(NameMC.Data());
 }
 
 outFile->cd();
//  hDATA.Write();
//  hMC.Write();
 
 double result = hMC->KolmogorovTest(hDATA,"M");
// double result = - hMC->KolmogorovTest(hDATA,"X");
//  double result = hMC->Chi2Test(&hDATA,"CHI2/NDF");
//=========> E' QUESTO! ==> double result = hMC->Chi2Test(hDATA,"CHI2");
 //double result = - hMC.Chi2Test(&hDATA,""); ///==== http://root.cern.ch/root/html/TH1.html#TH1:Chi2Test
 return result;
}

int main(int argc, char** argv)
{ 
 TDRStyle();
 
 gStyle->SetPadTopMargin(0.2);
 gStyle->SetPadBottomMargin(0.2);
 gStyle->SetPadLeftMargin(0.07);
 gStyle->SetPadRightMargin(0.23);
 gStyle->cd(); 
 
 
 std::cout << " " << std::endl;
 std::cout << " " << std::endl;
 std::cout << " " << std::endl; 
 std::cout << "     ___|          | _)  |         ___|                            |                   " << std::endl;
 std::cout << "   \\___ \\   __ \\   |  |  __|     \\___ \\    _` |  __ `__ \\   __ \\   |   _ \\     " << std::endl; 
 std::cout << "         |  |   |  |  |  |             |  (   |  |   |   |  |   |  |   __/             " << std::endl;
 std::cout << "   _____/   .__/  _| _| \\__|     _____/  \\__,_| _|  _|  _|  .__/  _| \\___|          " << std::endl;
 std::cout << "           _|                                              _|                          " << std::endl;
 std::cout << " " << std::endl;
 std::cout << " " << std::endl;
 std::cout << "    _ \\                      _ \\                                   " << std::endl;
 std::cout << "   |   |  |   |  __ \\       |   |   _` |  __ \\    _` |   _ \\      " << std::endl;
 std::cout << "   __ <   |   |  |   |      __ <   (   |  |   |  (   |   __/         " << std::endl;
 std::cout << "  _| \\_\\ \\__,_| _|  _|     _| \\_\\ \\__,_| _|  _| \\__, | \\___| " << std::endl;
 std::cout << "                                                |___/                " << std::endl;
 std::cout << " " << std::endl;
 std::cout << " " << std::endl;
 
 
 //Check if all nedeed arguments to parse are there                                                                                                                               
 if(argc != 2)
 {
  std::cerr << ">>>>> analysis.cpp::usage: " << argv[0] << " configFileName" << std::endl ;
  return 1;
 }


 // Parse the config file                                                                                                                                                          
 parseConfigFile (argv[1]) ;
 
 std::string treeName  = gConfigParser -> readStringOption("Input::treeName");
 std::string inputFileName = gConfigParser -> readStringOption("Input::inputFileName");

 std::vector<int> runRangesMin = gConfigParser -> readIntListOption("Options::runRangesMin");
 std::vector<int> runRangesMax = gConfigParser -> readIntListOption("Options::runRangesMax");
 
 TTree *treeJetLepVect; 
 
 std::string CutFile = gConfigParser -> readStringOption("Selections::CutFile"); 
 std::vector<std::string> vCut; ///====> only the first cut is used!
 std::cout << " nCuts   = " << ReadFileCut(CutFile, vCut) << std::endl;
 
 ///==== output file ====
 std::string outputDirectory = gConfigParser -> readStringOption("Output::outputDirectory");
 std::string OutFileName    = gConfigParser -> readStringOption("Output::outFileName");
 
 ///==== debug flag (begin) ==== 
 bool  debug = false; 
 try {
  debug = gConfigParser -> readBoolOption("Input::debug");
 }
 catch (char const* exceptionString){
  std::cerr << " exception = " << exceptionString << std::endl;
 }
 std::cout << ">>>>> input::debug  " << debug  << std::endl;  
 ///==== debug flag (end) ==== 
 
 
 
 
 ///==== program ====
 
 TFile* f = new TFile(inputFileName.c_str(), "READ");    
 treeJetLepVect = (TTree*) f->Get(treeName.c_str());
  
 int totNumEvents = treeJetLepVect->GetEntries(vCut.at(0).c_str());
 int numOutputFiles = runRangesMax.size();
 std::cout << " totNumEvents   = " << totNumEvents   << std::endl;
 std::cout << " numOutputFiles = " << numOutputFiles << std::endl;
 
 treeJetLepVect->SetEntryList(0); 
 treeJetLepVect->Draw(">> myList",vCut.at(0).c_str(),"entrylist");
 TEntryList *myList = (TEntryList*)gDirectory->Get("myList");
 treeJetLepVect->SetEntryList(myList); 
 
 std::cout << " tot = " << myList->GetN() << " = " << totNumEvents << " =? " << treeJetLepVect -> GetEntries () << std::endl;
 
 TH1F* MyHistoCounterMC = (TH1F*) f->Get("AllEvents/totalEvents");
 
 for (int iOutFile = 0; iOutFile < numOutputFiles; iOutFile++) {
  std::cout << " iOutFile = " << iOutFile << " : " << numOutputFiles << std::endl;
  TString outputRootFileName = Form("%s/%s_%d.root",outputDirectory.c_str(),OutFileName.c_str(), iOutFile); 
  TFile outputRootFile ( outputRootFileName.Data(), "RECREATE") ;
  outputRootFile.cd () ;
  outputRootFile.mkdir ("ntupleEcalAlignment") ;
  outputRootFile.cd ("ntupleEcalAlignment") ;  
  
  TTree* cloneTree = treeJetLepVect -> CloneTree (0) ;
  
//.........这里部分代码省略.........

//.........这里部分代码省略.........
 MyTreeMC->SetBranchAddress("MT",&MT);
 MyTreeMC->SetBranchAddress("EoP",&EoP);
 MyTreeMC->SetBranchAddress("eta",&eta);
 
 MyTreeMC->SetBranchAddress("E5x5",&E5x5);
 MyTreeMC->SetBranchAddress("p",&p);
 MyTreeMC->SetBranchAddress("eleES",&eleES);
 MyTreeMC->SetBranchAddress("eleFBrem",&eleFBrem);
 
 numEntriesMC = MyTreeMC->GetEntries();
 
 ///==== prepare minuit ====

 fitMin->SetRange(MinScanRange,MaxScanRange); 
 
 double step[1] = {0.001};
 double variable[1] = {0.0};
 minuit->SetLimitedVariable(0,"Scale" , variable[0]  , step[0] , MinScan  , MaxScan );
 

///===========================
///==== DATA Scale search ====
 ScaleTrue = -1000; ///==== default
 Data_or_MC = 1; ///=== 1 = Data;  0 = MC;
 numEvents = MyTreeDATA->GetEntries(); //==== number of events in Data sample
 outFile->cd();
 vET_data.clear();
 nIter = 1000000000; ///==== less than 1000000000 iterations at the end !!!
 TString nameDATA = Form("hDATA_%d_%d_%.5f",Data_or_MC,nIter,ScaleTrue);
 TH1F hDATA(nameDATA,nameDATA,numBINS,minBINS,maxBINS);
 
 MyTreeDATA->Draw(">> myList",(AdditionalCut + Form(" && ET > %f",minET)).Data(),"entrylist");
 TEntryList *mylist = (TEntryList*)gDirectory->Get("myList");
 MyTreeDATA->SetEntryList(mylist);
 
 MyTreeDATA->Draw(Form("%s >> %s",variableName.c_str(),nameDATA.Data()));
 ConvertStdVectDouble(vET_data,MyTreeDATA->GetV1(),mylist->GetN());
 
 hDATA.Write();
  
 std::cerr << "... I'm minimizing ... DATA analysis" << std::endl;
 std::cerr << ">>>>>>> numEvents = " << numEvents << " => " << vET_data.size() << " selected (=" << mylist->GetN() << ")" << std::endl;
 numSelectedData = vET_data.size();
 
 
 ///===== Chi2 ====
 std::cerr << " === Chi2 === " << std::endl;
 minuit->SetFunction(functorChi2);

 TGraph * grChi2 = new TGraph(iNoSteps);
 minuit->Scan(iPar_NoBG,iNoSteps,grChi2->GetX(),grChi2->GetY(),MinScan,MaxScan);

// TGraph * grChi2 = new TGraph();
// for (int iStep = 0; iStep < iNoSteps; iStep++){
//  double x = MinScan + (MaxScan - MinScan) / iNoSteps * (iStep+0.5);
//  double y = Chi2F(&x);
//  grChi2->SetPoint(iStep+1,x,y);
// }
 grChi2->Draw("AL");
 outFile->cd();
 minuit->PrintResults();
 outFile->cd();
 grChi2->SetTitle("grChi2");
 grChi2->Write();
 const double *outParametersTemp = minuit->X();
 const double *errParametersTemp = minuit->Errors();

///**** NewChi2 ****
void doMC_NewChi2(){
 TTree* MyTreeMC = (TTree*) fileInMC->Get(treeNameMC.c_str());
 for (nIter = 0; nIter<maxIter; nIter++){
  if (!(nIter%1)) std::cerr << ">>> nIter = " << nIter << " : " << maxIter << std::endl;
  vET_data.clear();
  outFile->cd();
  TString nameDATA = Form("hDATA_%d_%d_%.5f",Data_or_MC,nIter,ScaleTrue);
  TH1F hDATA(nameDATA,nameDATA,numBINS,minBINS,maxBINS);
  
  MyTreeMC->Draw(">> myListMC",(AdditionalCut+Form("&& (ET * (1+(%f)))>%f",ScaleTrue,minET)).Data(),"entrylist");
  TEntryList *myListMC = (TEntryList*)gDirectory->Get("myListMC");

  MyTreeMC->SetEntryList(0); 
  TEntryList *listMCHere = new TEntryList("listMCHere","listMCHere");    
  for (int iEvt = 0; iEvt < numSelectedData; iEvt ++){
   listMCHere->Enter(myListMC->GetEntry(gRandom->Uniform(0,myListMC->GetN())));
  }
 
 MyTreeMC->SetEntryList(listMCHere);
 MyTreeMC->Draw(Form("(1+%f) * %s >> %s",ScaleTrue,variableName.c_str(),nameDATA.Data()));
  
 ConvertStdVectDouble(vET_data,MyTreeMC->GetV1(),numSelectedData);
   
  ///==== newChi2 ====
  minuit->SetFunction(functorNewChi2);
  TGraph * grNewChi2 = new TGraph(iNoSteps);
  minuit->Scan(iPar_NoBG,iNoSteps,grNewChi2->GetX(),grNewChi2->GetY(),MinScan,MaxScan);
  grNewChi2->Draw("AL");
  outFile->cd();
  minuit->PrintResults();
  const double *outParametersNewChi2 = minuit->X();
  const double *errParametersNewChi2 = minuit->Errors();
 
  double minNewChi2 = grNewChi2->Eval(outParametersNewChi2[0]);
  ///==== end newChi2 ====
 
  ///==== Save the whole shape of LL/Chi2 ====
  for (unsigned int ii=0; ii < iNoSteps; ii++){
   double X_ii = (MaxScan - MinScan) / iNoSteps * ii + MinScan;
   
   Alpha   = X_ii;
   Chi2    = 0;
   LL      = 0;
   NewChi2 = grNewChi2->Eval(X_ii);
   
   myTreeChi2->Fill();
  }
  

  ///===== Look for minima =====
  double a;
  double b;
  double c;
  
  double errX_low = -9999;
  double errX_up = 9999;
  int err_low = 0;
  int err_up = 0;
  for (unsigned int ii=0; ii < iNoSteps; ii++){
   double X_ii = (MaxScan - MinScan) / iNoSteps * ii + MinScan;
   double here = grNewChi2->Eval(X_ii);
   if (err_low == 0){
    if (here < (minNewChi2 + DELTA_CHI2)){
     errX_low = X_ii;
     err_low = 1;
    }
   }
   else if (err_up == 0 && here > (minNewChi2 + DELTA_CHI2) && X_ii > outParametersNewChi2[0]){
    errX_up = X_ii; 
    err_up = 1;
   }
  }
  
  AlphaMean = outParametersNewChi2[0];
  AlphaMinus = errX_low;
  AlphaPlus = errX_up;
  
  
  grNewChi2->Fit("fitMin","RMQ");
  c = fitMin->GetParameter(0);
  b = fitMin->GetParameter(1);
  a = fitMin->GetParameter(2);
  AlphaMean_Fit  = -b / (2*a);
  AlphaMinus_Fit = (-b + 2 * sqrt(a)) / (2*a);  ///==== delta Chi2 = 1 
  AlphaPlus_Fit  = (-b - 2 * sqrt(a)) / (2*a);  ///==== delta Chi2 = 1   
  
  myTreeNewChi2_Result->Fill();
  
 // delete listMCHere;
 }
 
}

///**** LL ****
void doMC_LL(){
 TTree* MyTreeMC = (TTree*) fileInMC->Get(treeNameMC.c_str());

 for (nIter = 0; nIter<maxIter; nIter++){
  if (!(nIter%1)) std::cerr << ">>> nIter = " << nIter << " : " << maxIter << std::endl;
  vET_data.clear();
  outFile->cd();
  TString nameDATA = Form("hDATA_%d_%d_%.5f",Data_or_MC,nIter,ScaleTrue);
  TH1F hDATA(nameDATA,nameDATA,numBINS,minBINS,maxBINS);
  
  MyTreeMC->Draw(">> myListMC",(AdditionalCut+Form("&& (ET * (1+(%f)))>%f",ScaleTrue,minET)).Data(),"entrylist");
  TEntryList *myListMC = (TEntryList*)gDirectory->Get("myListMC");
  
  MyTreeMC->SetEntryList(0); 
  TEntryList *listMCHere = new TEntryList("listMCHere","listMCHere");    
  for (int iEvt = 0; iEvt < numSelectedData; iEvt ++){
   listMCHere->Enter(myListMC->GetEntry(gRandom->Uniform(0,myListMC->GetN())));
  }
  
  MyTreeMC->SetEntryList(listMCHere);
  MyTreeMC->Draw(Form("(1+%f) * %s >> %s",ScaleTrue,variableName.c_str(),nameDATA.Data()));
  
  ConvertStdVectDouble(vET_data,MyTreeMC->GetV1(),numSelectedData);
  
  ///==== likelihood ====
 std::cerr << " === LL === " << std::endl;
 std::cerr << " === pseudo vET_data.size() = " << vET_data.size() << std::endl;
 
 
  minuit->SetFunction(functorLL); 
  TGraph * grLL_temp = new TGraph(iNoSteps);
  minuit->Scan(iPar_NoBG,iNoSteps,grLL_temp->GetX(),grLL_temp->GetY(),MinScan,MaxScan);
  TGraph * grLL = new TGraph();
  int nPointLL = 0;
 for (unsigned int iStep = 0; iStep < iNoSteps; iStep++){
   double x = MinScan + (MaxScan - MinScan) / iNoSteps * (iStep+0.5);
   double y = LLFunc(&x);
   if (y != numberDATA * numEvents) {
    grLL->SetPoint(nPointLL,x,y);
    nPointLL++;
   }
  }
  grLL->Draw("AL");
  outFile->cd();
  minuit->PrintResults();
  const double *outParametersTemp2 = minuit->X();
  const double *errParametersTemp2 = minuit->Errors();
  
  double *outParametersLL = new double;
  double *errParametersLL = new double;
  outParametersLL[0] = outParametersTemp2[0];
  errParametersLL[0] = errParametersTemp2[0];
  
  std::cerr << " nPointLL = " << nPointLL << std::endl;
  
  double minLL = grLL->Eval(outParametersLL[0]);
  ///==== end likelihood ====
  
    ///==== Save the whole shape of LL/Chi2 ====
  for (unsigned int ii=0; ii < iNoSteps; ii++){
   double X_ii = (MaxScan - MinScan) / iNoSteps * ii + MinScan;
   
   Alpha   = X_ii;
   Chi2    = 0;
   LL      = grLL->Eval(X_ii);
   NewChi2 = 0;
   
   myTreeChi2->Fill();
  }

  ///===== Look for minima =====
  double a;
  double b;
  double c;
  
  
  double errX_low = -9999;
  double errX_up = 9999;
  int err_low = 0;
  int err_up = 0;
  for (unsigned int ii=0; ii < iNoSteps; ii++){
   double X_ii = (MaxScan - MinScan) / iNoSteps * ii + MinScan;
   double here = grLL->Eval(X_ii);
   if (err_low == 0){
    if (here < (minLL + DELTA_LL)){
     errX_low = X_ii;
     err_low = 1;
    }
   }
   else if (err_up == 0 && here > (minLL + DELTA_LL) && X_ii > outParametersLL[0]){
    errX_up = X_ii; 
    err_up = 1;
   }
  }
  
  AlphaMean = outParametersLL[0];
  AlphaMinus = errX_low;
  AlphaPlus = errX_up;
 
//.........这里部分代码省略.........

double LLFunc(const double *xx ){
 const Double_t scale = xx[0];

 ///==== Likelihood ==== 
 fileInMC->cd();
 TTree* MyTreeMC = (TTree*) fileInMC->Get(treeNameMC.c_str());
 TString NameMCGr = Form("grMC_%.5f",scale);
 TString NameMC = Form("hMC_%.5f",scale);
 TH1F* hMC;
// TGraphErrors* tempGrLL;
 if (!gROOT->FindObject(NameMC.Data())){
  hMC = new TH1F(NameMC,NameMC,NBINTemplate,MinTemplate,MaxTemplate);
  hMC->Reset();
  
  MyTreeMC->SetEntryList(0); 
  MyTreeMC->Draw(">> myListMCTot",(AdditionalCut + Form(" && (ET * (1+(%f)))>%f",scale,minET)).Data(),"entrylist");
  TEntryList *mylistMCTot = (TEntryList*)gDirectory->Get("myListMCTot");
  MyTreeMC->SetEntryList(mylistMCTot);
  
  TString DrawMC = Form("(%s * (1+(%f)))>>%s",variableName.c_str(),scale,NameMC.Data());
  std::cerr << " DrawMC = " << DrawMC.Data() << std::endl;
//   std::cerr << " LL CUT = " << (AdditionalCut+Form("&& (ET * (1+(%f)))>%f",scale,minET)).Data() << std::endl;
// MyTreeMC->Draw(DrawMC,(AdditionalCut+Form("&& (ET * (1+(%f)))>%f",scale,minET)).Data());
  MyTreeMC->Draw(DrawMC);
  hMC->Scale(10./numEvents);
  outFile->cd();
  hMC->Write();
//  tempGrLL = new TGraphErrors(buildGEfromH (*hMC));
//  tempGrLL->SetName(NameMCGr);
//  tempGrLL->Write();
 }
 else {
  std::cerr << " LL old " << NameMC.Data() << std::endl;
  hMC = (TH1F*) outFile->Get(NameMC.Data());
//  tempGrLL = (TGraphErrors*) outFile->Get(NameMCGr.Data());
 }
 numberDATA = vET_data.size();
 double result = 1.;
 
 std::cerr << " numberDATA = " << numberDATA << " hMC->GetEntries() = " << hMC->GetEntries() << std::endl;
 for (int iEvt = 0; iEvt < numberDATA; iEvt ++){
  double ET = vET_data.at(iEvt);
  int bin = ( ET - MinTemplate ) / Delta;
  if (bin > 0 && bin <= NBINTemplate){
  std::cerr << " ok here " << std::endl;
//   result *= (tempGrLL->Eval(hMC->GetBinCenter(bin)));
//    std::cerr << " hMC->GetBinContent(" << bin << ") = " << hMC->GetBinContent(bin) << " result = " << result << std::endl;
   if (hMC->GetBinContent(bin) != 0) {
    result -= log(hMC->GetBinContent(bin));
   }
   else {
    result = numberDATA * numEvents;
    std::cerr << " >>>>>>>>>>>>>>>> it's ZERO !!!! " << std::endl; //====== > too verbose!
    return result;
   }
//    if (hMC->GetBinContent(bin) == 0) {
//     std::cerr << " result = " << result << " hMC.GetBinContent(" << bin << ":" << NBINTemplate << ") = " << hMC.GetBinContent(bin) << " scale = " << scale << std::endl;
//    }
  }
 }
//  outFile->cd();
//  hMC.Write();
 
//  if (result != 0) result = -log(result);
//  else {
//   result = numberDATA * numEvents;
//   std::cerr << " it's ZERO !!!! " << std::endl;
//  }
 
 ///==== end Likelihood ====    

 return result;
}

void plotL1Digis(Long64_t entryNum=1)
{
  TH1::SetDefaultSumw2();
  
  TString inputFile = "L1UpgradeAnalyzer.root";
  TString inputFileunp = "L1UnpackedUnpacker.root";
  
  TFile *inFile = TFile::Open(inputFile);
  TFile *inFileunp = TFile::Open(inputFileunp);
  
  TTree *emulatorResults = (TTree*)inFile->Get("L1UpgradeAnalyzer/L1UpgradeTree");
  TTree *unpackerResults = (TTree*)inFileunp->Get("UnpackerResults/L1UpgradeTree");
  
  //TCut entryCut = Form("Entry$ == %i", entryNum);
  //TCut entryCutUnpacker = Form("Entry$ == %i",entryNum+179);

  TEntryList* elist = new TEntryList(unpackerResults);
  
  for ( int bx(0); bx<100; bx++) {
	if ( !unpackerResults->Draw("Entry$",Form("FEDBXID == %i", bx), "nodraw" ) ) continue;
  	
  	elist->Enter(unpackerResults->GetV1()[0], unpackerResults);
  }
  
  unpackerResults->SetEntryList(elist);
  
  elist->Print();

  const int nHISTS = 30;
  TString labels[nHISTS] = {"region_et", "region_eta", "region_phi",
			    "egcand_rank", "egcand_eta", "egcand_phi",
			    "central_jet_hwPt", "central_jet_hwEta", "central_jet_hwPhi",
			    "forward_jet_hwPt", "forward_jet_hwEta", "forward_jet_hwPhi",
			    "ETT", "HTT",
			    "MET_Rank", "MET_Phi", "MHT_Rank", "MHT_Phi",
			    "iso_egamma_hwPt", "iso_egamma_hwEta", "iso_egamma_hwPhi",
			    "noniso_egamma_hwPt", "noniso_egamma_hwEta", "noniso_egamma_hwPhi",
			    "tau_hwPt", "tau_hwEta", "tau_hwPhi",
			    "isotau_hwPt", "isotau_hwEta", "isotau_hwPhi"};

  TString projectionnames[nHISTS] = {"legacyregion_et", "legacyregion_gctEta", "legacyregion_gctPhi",
				     "legacyemcand_rank", "legacyemcand_regionEta", "legacyemcand_regionPhi",
				     "jet_hwPt", "jet_hwEta", "jet_hwPhi",
				     "jet_hwPt", "jet_hwEta", "jet_hwPhi",
				     "etsum_hwPt", "etsum_hwPt",
				     "etsum_hwPt","etsum_hwPhi","etsum_hwPt","etsum_hwPhi",
				     "egamma_hwPt", "egamma_hwEta", "egamma_hwPhi",
				     "egamma_hwPt", "egamma_hwEta", "egamma_hwPhi",
			         "tau_hwPt", "tau_hwEta", "tau_hwPhi",
			         "isotau_hwPt", "isotau_hwEta", "isotau_hwPhi"};

  TCut projectioncuts[nHISTS] = {"", "", "",
				 "", "", "",
				 "(jet_hwQual&0x2)!=0x2&&jet_bx==0","(jet_hwQual&0x2)!=0x2&&jet_bx==0","(jet_hwQual&0x2)!=0x2&&jet_bx==0",
				 "(jet_hwQual&0x2)==0x2&&jet_bx==0","(jet_hwQual&0x2)==0x2&&jet_bx==0","(jet_hwQual&0x2)==0x2&&jet_bx==0",
				 "etsum_type==0&&etsum_bx==0","etsum_type==1&&etsum_bx==0",
				 "etsum_type==2&&etsum_bx==0","etsum_type==2&&etsum_bx==0","etsum_type==3&&etsum_bx==0","etsum_type==3&&etsum_bx==0",
				 "egamma_hwIso==1&&egamma_bx==0", "egamma_hwIso==1&&egamma_bx==0", "egamma_hwIso==1&&egamma_bx==0",
				 "egamma_hwIso==0&&egamma_bx==0", "egamma_hwIso==0&&egamma_bx==0", "egamma_hwIso==0&&egamma_bx==0",
				 "tau_bx==0", "tau_bx==0", "tau_bx==0",
				 "isotau_bx==0", "isotau_bx==0", "isotau_bx==0"};
  Int_t minBin[nHISTS] = {0, 0, 0,
			  0, 0, 0,
			  0, 0, 0,
			  0, 0, 0,
			  0, 0,
			  0, 0, 0, 0,
  			  0, 0, 0,
			  0, 0, 0,
			  0, 0, 0,
			  0, 0, 0};
  Int_t maxBin[nHISTS] = {40, 22, 25,
			  64, 22, 25,
			  64,22,25,
			  64,22,25,
			  600, 300,
			  200, 70, 200, 30,
			  64, 22, 25,
			  64, 22, 25,
			  64, 22, 25,
			  64, 25, 25};

  TH1I *hists[nHISTS][2];
  TH1D *divs[nHISTS];
  TCanvas *c[nHISTS];

  for(int i = 0; i < nHISTS; ++i)
  {
    hists[i][0] = new TH1I(labels[i], ";"+labels[i], maxBin[i]-minBin[i], minBin[i], maxBin[i]);
    hists[i][1] = (TH1I*)hists[i][0]->Clone(labels[i]+"unpacked");
    divs[i] = new TH1D(labels[i]+"div", ";"+labels[i], maxBin[i]-minBin[i], minBin[i], maxBin[i]);

    //emulatorResults->Project(hists[i][0]->GetName(), projectionnames[i], projectioncuts[i]&&entryCut);
    //unpackerResults->Project(hists[i][1]->GetName(), projectionnames[i], projectioncuts[i]&&entryCutUnpacker);

    emulatorResults->Project(hists[i][0]->GetName(), projectionnames[i], projectioncuts[i]);
    unpackerResults->Project(hists[i][1]->GetName(), projectionnames[i], projectioncuts[i]);

    divs[i]->Divide(hists[i][1], hists[i][0]);
    for (int m=1;m<=hists[i][1]->GetNbinsX();m++){
//.........这里部分代码省略.........

int main() {
	Int_t nbins = 800, count = 0;
    Double_t integral;
	TFile *input = new TFile("FilterRMSComparison.root");
    TFile *templatefile = new TFile("/home/marko/Desktop/H4Analysis/ntuples/Templates_APDs.old.root");
	TTree *MyTree = (TTree*) input->Get("RMS");
    TCanvas *can1 = new TCanvas("can1", "canvas", 1200,900);
    TCanvas *can2 = new TCanvas("can2", "canvas", 1200,1200);
    TCanvas *can3 = new TCanvas("can3", "canvas", 1200,1200);
    can1->Divide(1,3);
    can2->Divide(1,2);
    can3->Divide(1,2);
    MyTree->SetEntryList(0);
    TString listcut = "abs(unfilteredbslope)<6 && unfilteredampfit>700";
    MyTree->Draw(">>myList", listcut, "entrylist");
    TEntryList *myList = (TEntryList*) gDirectory->Get("myList");
    MyTree->SetEntryList(myList);
    Int_t nevents = myList->GetN();
	MyTree->Draw("unfilteredevent", "abs(unfilteredbslope)<6 && unfilteredampfit>700", "goff");
    Double_t *vTemp = MyTree->GetV1();
    Int_t *vEvent = new Int_t[nevents];
    for (int iEntry = 0; iEntry<nevents; iEntry++) {
        vEvent[iEntry] = vTemp[iEntry];
    }
	TString plot, plot2, cut;
	char name[50];
    TH1F *histoave = new TH1F("histoave","Wave Pulse Average", nbins, -40, 120);
    TH1F *histoavefft = new TH1F("histoavefft","Wave Pulse Average FFT", nbins, 0, 5);
    TH1F *histoaveph = new TH1F("histoaveph","Wave Pulse Average Phase", nbins, 0, 800);
    TH1F *originaltemplate = (TH1F*) templatefile->Get("APD2_E50_G50_prof");
    originaltemplate->Rebin(16);
    TH1F *templatehisto = new TH1F("templatehisto", "Template = Green, Average = Red", nbins, -40, 120); 
    TH1F *difference = new TH1F("difference","Template vs. Average Difference", nbins, -40, 120);
    TH1F *percentdifference = new TH1F("percentdifference","Template vs. Average Percent Difference", nbins, -40, 120);
    for (Int_t i=0;i<nbins;i++) {
        templatehisto->SetBinContent(i+1, originaltemplate->GetBinContent(i+1));
    }
    templatehisto->SetLineColor(kGreen+3);
    templatehisto->SetLineWidth(4);
    can2->cd(2);
    histoave->GetYaxis()->SetRangeUser(-.120,1.2);
    //templatehisto->Draw();
    histoave->SetStats(0);
    histoave->Draw();
	for (Int_t i=0;i<nevents;i++) {
    //for (Int_t i=0;i<10;i++) {
        if (vEvent[i]==513) continue; //event for which electronics die out for a bit halfway through
		TString histoname = "TempHisto_";
        histoname += i;
        TString histoname2 = "TempHisto2_";
        histoname2 += i;
        TH2F* TempHisto = new TH2F (histoname, "Temp Histo", nbins, -40, 120, 1000, -120, 800); //nanoseconds
        TH2F* TempHisto2 = new TH2F (histoname2, "Temp Histo", nbins, -40, 120, 1000, -120, 800); //nanoseconds
        TString h1name = "h1001_";
        h1name += i;
        TString h1name2 = "h1002_";
        h1name2 += i;
        TString h1name2fft = "h1002fft_";
        h1name2 += i;
        TString h1name2ph = "h1002ph_";
        h1name2 += i;
        TString h1name3 = "h1003_";
        h1name3 += i;
        TString h1name4 = "h1004_";
        h1name4 += i;
        TH1F *h1001 = new TH1F(h1name,"Red = Unfiltered, Blue = Filtered", nbins, -40, 120);
        TH1F *h1002 = new TH1F(h1name2,"h1002", nbins, -40, 120);
        TH1F *h1002fft = new TH1F(h1name2fft,"h1002fft", nbins, 0, 5);
        TH1F *h1002ph = new TH1F(h1name2ph,"h1002ph", nbins, 0, 800);
        TH1F *h1003 = new TH1F(h1name3,"Filtered WF - Unfiltered WF", nbins, -40, 120);
        TH1F *h1004 = new TH1F(h1name4,"Percent Change in Filtered WF - Unfiltered WF", nbins, -40, 120);
        plot = "unfilteredwfval:(unfilteredwftime-unfilteredtimeref)>>";
        plot += histoname;
        plot2 = "filteredwfval:(filteredwftime-filteredtimeref)>>";
        plot2 += histoname2;
        cut = "abs(unfilteredbslope)<6 && unfilteredampfit>700 && unfilteredevent==";
        cut += vEvent[i];
        MyTree->Draw(plot, cut, "goff");
        TempHisto = (TH2F*) gDirectory->Get(histoname);
        h1001 = transform2Dto1D(TempHisto);
        MyTree->Draw(plot2, cut, "goff");
        TempHisto2 = (TH2F*) gDirectory->Get(histoname2);
        h1002 = transform2Dto1D(TempHisto2);
        if (h1002->GetBinCenter(h1002->GetMaximumBin()) < 30) continue;
        sprintf(name, "Good Events/Event%d", vEvent[i]);
        strcat(name, ".png");
        h1001->SetLineColor(kRed);
        h1002->SetLineColor(kBlue);
        h1001->SetStats(0);
        h1002->SetStats(0);
        //can1->cd(1);
        can1->cd();
        h1001->GetXaxis()->SetTitle("Time (ns)");
        h1001->GetYaxis()->SetTitle("Amplitude");
        h1001->Draw();
        h1002->Draw("same");
        //for (Int_t i=0;i<nbins;i++) {
        //    h1003->SetBinContent(i+1, (h1002->GetBinContent(i+1))-(h1001->GetBinContent(i+1)));
        //    if (h1001->GetBinContent(i+1) != 0) h1004->SetBinContent(i+1, ((h1002->GetBinContent(i+1))-(h1001->GetBinContent(i+1)))/(h1001->GetBinContent(i+1)));
        //    else h1004->SetBinContent(i+1, 0);
//.........这里部分代码省略.........