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

void trigtime(int run, float tmax =160) {
  ostingstream ssrun;
  ssrun << run;
  string srun = srun;
  ssfname << "run" << run << "evt.root";
  string fname = "run" + srun + "evt.root";
  TFile* pfile = new TFile(fname.c_str(), "READ");
  if ( pfile == 0 || !pfile->IsOpen()) {
    cout << "File not found: " << fname << endl;
    return;
  }
  TTree* ptree = dynamic_cast<TTree*>(pfile->Get("DXDisplay/EventTree"));
  if ( ptree == 0 ) {
    cout << "Tree not found." << endl;
    return;
  }
  new TCanvas;
  string stit = "Trigger times; Time [sec]; Trigger (109 is ghost)";
  string sarg = "trigger-(trigger>100)*109+109:tlo-1456885145";
  string sarg = "trigger-(trigger>100)*109+109:tlo-1456885145";
  TH2* ph2 = new TH2F("h2", stit.c_str(), 100, 0, tmax, 8, 107.5, 116.5);
  ph2->SetStats(0);
  ph2->Draw();
  ptree->SetMarkerStyle(2);
  ptree->Draw(sarg.c_str(), "", "same");
}

void drawCutTrees(const char * filter="sinsub_10_3_ad_2")
{

  TFile *f = new TFile (TString::Format("thermalCutTrees_%s.root",filter)); 

  TTree * signal = (TTree*) f->Get("signal_in"); 
  TTree * bg = (TTree*) f->Get("bg_in"); 
  signal->SetMarkerColorAlpha(3,0.05); 
  signal->SetMarkerStyle(1); 

  TCanvas * c = new TCanvas(filter,filter,1000,1000); 
 
  c->Divide(2,2); 

  c->cd(1); 
  gPad->SetLogz(); 
  bg->Draw("coherentHilbertPeak:mapPeak >> h1(100,0,0.5,100,0,300)","","colz"); 
  signal->Draw("coherentHilbertPeak:mapPeak","","psame"); 

  c->cd(2); 
  bg->Draw("coherentHilbertPeak:mapPeak >> h2(100,0,0.5,100,0,300)","dPhiSun > 10","colz"); 
  signal->Draw("coherentHilbertPeak:mapPeak","","psame"); 
  gPad->SetLogz(); 

  c->cd(3); 
  bg->Draw("coherentHilbertPeak:mapPeak >> h3(100,0,0.5,100,0,300)","dPhiNorth > 90","colz"); 
  signal->Draw("coherentHilbertPeak:mapPeak","","psame"); 
  gPad->SetLogz(); 

  c->cd(4); 
  bg->Draw("coherentHilbertPeak:mapPeak >> h4(100,0,0.5,100,0,300)","dPhiNorth > 90 && dPhiSun > 10","colz"); 
  signal->Draw("coherentHilbertPeak:mapPeak","","psame"); 
  gPad->SetLogz(); 

  c = new TCanvas("deconv","deconv",1000,1000); 
 
  c->Divide(2,2); 

  c->cd(1); 
  gPad->SetLogz(); 
  bg->Draw("deconvHilbertPeak:mapPeak >> h1(100,0,0.5,100,0,300)","","colz"); 
  signal->Draw("deconvHilbertPeak:mapPeak","","psame"); 

  c->cd(2); 
  bg->Draw("deconvHilbertPeak:mapPeak >> h2(100,0,0.5,100,0,300)","dPhiSun > 10","colz"); 
  signal->Draw("deconvHilbertPeak:mapPeak","","psame"); 
  gPad->SetLogz(); 

  c->cd(3); 
  bg->Draw("deconvHilbertPeak:mapPeak >> h3(100,0,0.5,100,0,300)","dPhiNorth > 90","colz"); 
  signal->Draw("deconvHilbertPeak:mapPeak","","psame"); 
  gPad->SetLogz(); 

  c->cd(4); 
  bg->Draw("deconvHilbertPeak:mapPeak >> h4(100,0,0.5,100,0,300)","dPhiNorth > 90 && dPhiSun > 10","colz"); 
  signal->Draw("deconvHilbertPeak:mapPeak","","psame"); 
  gPad->SetLogz(); 


}

void evtime() {
  vector<int> runs;
  vector<int> marks;
  vector<int> mcols;
  runs.push_back(13893);
  marks.push_back(2);
  mcols.push_back(2);
  runs.push_back(14009);
  marks.push_back(24);
  mcols.push_back(4);
  runs.push_back(14085);
  marks.push_back(25);
  mcols.push_back(3);
  runs.push_back(14234);
  marks.push_back(5);
  mcols.push_back(3);
  runs.push_back(14434);
  marks.push_back(28);
  mcols.push_back(3);
  string pre = "run";
  string suf = "evt.root";
  double tmax = 160;
  double emax = 180;
  new TCanvas;
  TH1* ph = new TH2F("hevtime", "Event vs. time; Time [sec]; Event",
                     tmax, 0, tmax, emax, 0, emax);
  ph->SetStats(0);
  ph->Draw("axis");
  TLegend* pleg = new TLegend(0.65,0.15,0.85,0.37);
  pleg->SetBorderSize(0);
  for ( unsigned int isam=0; isam<runs.size(); ++isam ) {
    int run = runs[isam];
    int icol = dsindex(run);
    int col = colormap(icol);
    ostringstream ssrun;
    ssrun << run;
    string srun = ssrun.str();
    ostringstream sst0;
    sst0 << t0map(run);
    string fname = pre + srun + suf;
    TFile* pfile = TFile::Open(fname.c_str(), "READ");
    TTree* ptree = dynamic_cast<TTree*>(pfile->Get("DXDisplay/EventTree"));
    ptree->SetMarkerStyle(marks[isam]);
    ptree->SetMarkerColor(col);
cout << run << " " << col << " " << sst0.str() <<  endl;
    if ( ptree == 0 ) {
      cout << "Tree not found" << endl;
      pfile->ls();
      return;
    }
    string sarg = "event:tlo-";
    sarg += sst0.str();
    sarg += ">>hevtime";
    ptree->Draw(sarg.c_str(), "", "same");
    pleg->AddEntry(ptree, srun.c_str(), "p");
  }
  pleg->Draw();
}

void plot(){

  //TH2F *hframe = new TH2F("hframe","Edep_tot vs. track length p=180MeV/C",50,0,50,100,178.0,179.0);
  //TH2F *hframe = new TH2F("hframe","Edep_tot vs. track length p=120MeV/c",50,0,50,100,115.2,116.8);
  //  TH2F *hframe = new TH2F("hframe","Edep_tot vs. track length p=90MeV/c",500,0,50,100,79.5,83.0);
  TH2F *hframe = new TH2F("hframe","Edep_tot vs. track length p=72MeV/c",500,0,50,100,45.5,56.1);
  //  TH2F *hframe = new TH2F("hframe","Edep_tot vs. track length p=65MeV/c",500,0,50,100,0.5,46.1);
  TTree *MyTree = new TTree("MyTree", "MyTree");
  MyTree->ReadFile("Rout_proton72mevOnSolenoid3LIV.txt","track:edep");
  MyTree->SetMarkerStyle(kCircle);
  MyTree->SetMarkerColorAlpha(kRed, 0.35);
  hframe->Draw(); //you can set the axis att via hframe->GetXaxis()..
  MyTree->Draw("edep:track","","same");


  /*
  TTree *MyTree1 = new TTree("MyTree1", "MyTree1");
  MyTree1->ReadFile("Rout_proton120mevOnSolenoid3STD.txt","track1:edep1");
  MyTree1->SetMarkerStyle(kOpenSquare);
  MyTree1->SetMarkerColorAlpha(kBlue, 0.35);
  MyTree1->Draw("edep1:track1","","same");
  */
}

void P_vs_r(void)
{
	gROOT->Reset();
	
	TFile *f = new TFile("hd_root.root");
	TTree *geant = (TTree*)gROOT->FindObject("geant");
	TTree *dana = (TTree*)gROOT->FindObject("dana");

	TCanvas *c1 = new TCanvas("c1");
	c1->SetTicks();
	c1->SetGrid();

	geant->SetMarkerStyle(6);
	geant->SetMarkerSize(0.1);

	geant->Draw("P:sqrt(x*x+y*y)", "P>0.95");
	dana->Draw("P:sqrt(x*x+y*y)", "", "same");
	
	TLatex *lab = new TLatex(400.0, 1.005, "single 1GeV/c proton #theta=5^{o} #phi=180^{o}");
	lab->SetTextSize(0.03);
	lab->Draw();
	
	geant->SetLineColor(geant->GetMarkerColor());
	dana->SetLineColor(dana->GetMarkerColor());
	geant->SetLineWidth(4.0);
	dana->SetLineWidth(4.0);
	
	TLegend *leg = new TLegend(0.5, 0.5, 0.7, 0.7);
	leg->SetFillColor(kWhite);
	leg->AddEntry(geant, "GEANT");
	leg->AddEntry(dana, "DANA");
	leg->Draw();

	c1->SaveAs("P_vs_r.pdf");
	c1->SaveAs("P_vs_r.gif");
}

void kf_plot(){

  TH1F *hframe = new TH1F("hframe","momentum resolution",100,-10.,10.);
  hframe->Draw(); //you can set the axis att via hframe->GetXaxis()..
  gPad->SetLogx(0);
  hframe->GetXaxis()->SetTitle(" #Delta p_T /p_T [%]");
  hframe->GetYaxis()->SetTitle("Number of event");
  TLatex Tl; Tl.SetTextFont(43); Tl.SetTextSize(20);

  double pTresol,pXresol,pYresol;
  int nline =0;
  TTree *MyTree = new TTree("MyTree", "MyTree");
  MyTree->ReadFile("helix_out.txt","pTresol:pXresol:pYresol");
  MyTree->SetMarkerStyle(kCircle);
  MyTree->SetMarkerColorAlpha(kRed, 0.3);
  while(1){
    hframe->Fill(pTresol);
    nline++;
  }
  hframe->Draw("same");

  Tl.DrawLatex(10, 142,   "open: p_T resolution(helix)");

}

void striptree(TString fileName="", Int_t nAPVs=4)
{
  //===========================================================================
  // Change the inputs
  //===========================================================================

  //TString fileName = "Raw_Data_July1_Ped_300V";
  //TString fileName = "Raw_Data_July1_Source_300V";
  //TString fileName = "Raw_Data_July1_Ped_200V";
  //TString fileName = "Raw_Data_Pedestal_300V";
  // TString fileName = "Raw_Data_Source_300V";

  // TString fileName = "Raw_Data_Sept_3_1";
  // TString fileName = "Raw_Data_test_3_7";
  // TString fileName = "Raw_Data_Albox_2011_03_31_1";
  //-- TString fileName = "Raw_Data_04_01_AL_center_1";
  // TString fileName = "Raw_Data_04_01_AL_7820_1";
  //-- TString fileName = "Raw_Data_7817_04_01_1";
  // TString fileName = "Raw_Data_7817_04_01_bkg_1";

  //int nAPVs = 4;

  //TString fileName = "Raw_Data_July1_Source_200V";
  //int nAPVs = 1;

  int nEvents = -1; // -1 means all

  TString dat = fileName;
  TString eps = fileName + "-events.eps";
  TString root = fileName + "-events.root";

  TFile* file = new TFile(root, "RECREATE");

  //===========================================================================
  // Read the data and make a 2D scatter plot
  //===========================================================================

  ifstream in(dat);
  if (!in) {
     cout<< "File " << dat << " not found" <<endl;
     return;
  }

  TString comments = "vertical --> APV0, APV1, ... of event 1, ...";
  TString s;
  while (1) {
    s.ReadLine(in);
    if (s.BeginsWith(comments))
      break;
  }
  gStyle->SetOptStat(10);
  TString title = fileName+"   ADC Counts vs Channel";
  TString title32 = fileName+"   ADC Counts vs Channel for 32 ch";
  TH2D *h2d = new TH2D("h2d", title.Data(), nAPVs*128, 0, nAPVs*128, 200, 0, 200);
  TH2D *h2d32 = new TH2D("h2d32", title.Data(), 16, 0, 16, 200, 0, 200);
  TProfile *profile = new TProfile("profile", title32.Data(), nAPVs*128, 0, nAPVs*128);

  int event, apv, channel, adc;
  // TTree* tree = new TTree("tree", "tree");
  // tree->Branch("event",   &event,   "event/I");
  // tree->Branch("apv",     &apv,     "apv/I");
  // tree->Branch("channel", &channel, "channel/I");
  // tree->Branch("adc",     &adc,     "adc/I");

  Int_t raw[512];
  TTree* etree = new TTree("etree", "etree");
  etree->Branch("raw",   &raw,   "raw[512]/I");    // etree->Draw("raw:Iteration$","Entry$==0")
  etree->SetMarkerStyle(6);
  etree->SetMarkerColor(2);

  int n = 0;
  while (1) {
    in >> s;
    if (!in.good())
      break;
    if (s.BeginsWith("]DATA]")) // last line
      break;
    n++;
    if ((nEvents != -1) && (n > nEvents))
      break;
    event = atoi(s.Data());
    if (event%100 == 0)
      cout << "Processing event " << event << endl;
    for (int i=0; i<2; i++)
      in >> s; // time

    Int_t ichannel = 0;
    for (apv=1; apv<=nAPVs; apv++) {
      in >> s; // header + analog APV data + one tick mark
      int sequenceNumber = 0;
      TString subString = "";
      int length = s.Length();
      for (int j=0; j<length; j++) {
        char c = s[j];
        if (c != ',') {
          subString += c;
        } else {
          if (sequenceNumber >= 12) {
            channel = sequenceNumber - 12;
            adc = atoi(subString.Data());
//.........这里部分代码省略.........

void process()
{
   Int_t raw[512];   // buffer for input signal and bkg trees
   // buffers for output trees
   Int_t sig[512];
   Int_t cmsig[512];
   Int_t cm[16];

   // pedestal
   const char* fbkg_name = "Raw_Data_FZ320P_05_MSSD_2_250V_K237_Pedestal.dat-events.root";
   TFile* fbkg = TFile::Open(fbkg_name);
   if (!fbkg) cout<< "File not found: " << fbkg <<endl<<exitl;

   TTree* tree = (TTree*) fbkg->Get("etree");
   tree->SetBranchAddress("raw", &raw);

   TH2* h2d = (TH2*) fbkg->Get("h2d");
   new TCanvas;
   h2d->Draw();

   TProfile* profile = (TProfile*) fbkg->Get("profile"); 
   //new TCanvas;
   //profile->Draw();

   Int_t pedestal[512];
   for (int i=0; i<512; i++) {
      // pedestal[i] = profile->GetBinContent(i+1) - 0.5;
      // pedestal[i] = profile->GetBinContent(i+1) + 0.5;
      pedestal[i] = profile->GetBinContent(i+1);
   }

   // cout << "\nPedestals for every channel\n" << endl;
   // for (int i=0; i<512; i++) {
   //    cout << pedestal[i] << " ";
   //    if (i>0 && (i+1)%128==0)
   //       cout << endl;
   // }

   cout<< "processing bkg" <<endl;

   // output file with tree
   const char* obfname = "FZ320P_05_MSSD_2-bkg.root";
   TFile* obfile = TFile::Open(obfname, "recreate");

   TTree* btree = new TTree("btree", "btree");
   btree->Branch("sig", &sig, "sig[512]/I");
   btree->Branch("cmsig", &cmsig, "cmsig[512]/I");
   btree->Branch("cm", &cm, "cm[16]/I");
   btree->SetMarkerStyle(6);
   btree->SetMarkerColor(2);

   for (int jentry=0; jentry<tree->GetEntries(); ++jentry)
   {
      tree->GetEvent(jentry);
      // sig
      for (int i=0; i<512; ++i) {
         sig[i] = raw[i] - pedestal[i];
      }
      // calc common mode
      Int_t group32[32];
      Int_t index32[32];
      for (int igroup=0; igroup<16; ++igroup) {
         for (int istrip=0; istrip<32; ++istrip)   // istrip is number inside group of 32
         {
            group32[istrip] = sig[igroup*32 + istrip];
         }
         // sort array group32 in ascending order
         TMath::Sort(32, group32, index32, kFALSE);
         Int_t median = group32[index32[14]];
         cm[igroup] = median;
      }
      // subtract common mode
      for (int istrip=0; istrip<512; ++istrip) {
         Int_t igroup = istrip/32;
         cmsig[istrip] = sig[istrip] - cm[igroup];
      }
      // Fill sig, cmsig, cm
      btree->Fill();
   }
   obfile->Write();

   /////////////////////////////////////////////////
   //
   // signal tree
   //
   /////////////////////////////////////////////////

   cout<< "processing signal" <<endl;
   
   TChain* chain = new TChain("etree");
   chain->Add("Raw_Data_FZ320P_05_MSSD_250V_K237_Position_1.dat-events.root");
   chain->Add("Raw_Data_FZ320P_05_MSSD_250V_K237_Position_2.dat-events.root");
   chain->Add("Raw_Data_FZ320P_05_MSSD_2_250V_K237_Position_3.dat-events.root");
   chain->Add("Raw_Data_FZ320P_05_MSSD_2_250V_K237_Position_4.dat-events.root");
   chain->SetBranchAddress("raw", &raw);

   // output file with tree
   const char* osfname = "FZ320P_05_MSSD_2-signal.root";
   TFile* osfile = TFile::Open(osfname, "recreate");

//.........这里部分代码省略.........

void pulseConvert(const char* ifname)
{
   TFile* ifile = TFile::Open(ifname);
   if (!ifile) {
      cout<< "File not found: " << ifname <<endl;
      return;
   }
   TTree* itree = (TTree*) ifile->Get("pulse");
   if (!itree) {
      cout<< "Error: could not find tree \"pulse\"" <<endl;
      return;
   }

   Int_t event, tc1, tc2;
   Float_t b1_t[1024], b1_c[4096], b2_t[1024], b2_c[4096];

   Float_t* b1_c1 = b1_c;
   Float_t* b1_c2 = b1_c + 1024;
   Float_t* b1_c3 = b1_c + 2048;
   Float_t* b1_c4 = b1_c + 3072;
   Float_t* b2_c1 = b2_c;
   Float_t* b2_c2 = b2_c + 1024;
   Float_t* b2_c3 = b2_c + 2048;
   Float_t* b2_c4 = b2_c + 3072;

   itree->SetBranchAddress("event", &event);
   itree->SetBranchAddress("tc1", &tc1);
   itree->SetBranchAddress("tc2", &tc2);
   itree->SetBranchAddress("b1_t", &b1_t);
   itree->SetBranchAddress("b2_t", &b2_t);
   itree->SetBranchAddress("b1_c", &b1_c);
   itree->SetBranchAddress("b2_c", &b2_c);

   TFile* ofile = TFile::Open(Form("%s.pulse.root",ifname), "recreate");

   // TTree* otree = new TTree("pulse", "old pulse tree");
   TTree* otree = new TTree("p", "new pulse tree");
   otree->SetMarkerStyle(6);
   otree->SetMarkerColor(2);
   otree->SetLineColor(2);

   // otree->Branch("event", &event, "event/I");
   // otree->Branch("tc1", &tc1, "tc1/I");
   // otree->Branch("b1_t",  b1_t, "b1_t[1024]/F");
   // otree->Branch("b1_c1", b1_c1, "b1_c1[1024]/F");
   // otree->Branch("b1_c2", b1_c2, "b1_c2[1024]/F");
   // otree->Branch("b1_c3", b1_c3, "b1_c3[1024]/F");
   // otree->Branch("b1_c4", b1_c4, "b1_c4[1024]/F");
   // otree->Branch("tc2", &tc2, "tc2/I");
   // otree->Branch("b2_t",  b2_t, "b2_t[1024]/F");
   // otree->Branch("b2_c1", b2_c1, "b2_c1[1024]/F");
   // otree->Branch("b2_c2", b2_c2, "b2_c2[1024]/F");
   // otree->Branch("b2_c3", b2_c3, "b2_c3[1024]/F");
   // otree->Branch("b2_c4", b2_c4, "b2_c4[1024]/F");

   otree->Branch("event", &event, "event/I");
   otree->Branch("tc1", &tc1, "tc1/I");
   otree->Branch("t1",  b1_t, "t1[1024]/F");
   otree->Branch("c1", b1_c1, "c1[1024]/F");
   otree->Branch("c2", b1_c2, "c2[1024]/F");
   otree->Branch("c3", b1_c3, "c3[1024]/F");
   otree->Branch("c4", b1_c4, "c4[1024]/F");
   otree->Branch("tc2", &tc2, "tc2/I");
   otree->Branch("t2",  b2_t, "t2[1024]/F");
   otree->Branch("c5", b2_c1, "c5[1024]/F");
   otree->Branch("c6", b2_c2, "c6[1024]/F");
   otree->Branch("c7", b2_c3, "c7[1024]/F");
   otree->Branch("c8", b2_c4, "c8[1024]/F");

   for (Long64_t jentry=0; jentry<itree->GetEntries(); jentry++)
   {
      if (jentry % 1000 == 0) cout<< "jentry = " << jentry <<endl;
      itree->LoadTree(jentry);
      itree->GetEntry(jentry);

      otree->Fill();
   }

   cout<< "Write " << otree->GetEntries() << " events into file " << ofile->GetName() <<endl;
   ofile->Write();
}

void drawtorroidmollers_rings(TString infilename="0", TString outfilenamestem="0")
{
	if (infilename=="0" || outfilenamestem=="0") {
		printf("Usage:\n\t.x drawtorroidmollers_rings_stamps.C(infilename, outfilename)\n\n");
		printf("Where infilename is a .root file \n");
		printf("  and outfilenamestem is the base name of an image file.\n");
		return;
	}  

	gROOT->Reset();
	char inifilename[200];//, outfilenamestem[200], parttype[100];

	Bool_t debug=0;

//	Bool_t doproton=0;
//	Bool_t dopostscript=0;
	Double_t deltaZwindowsize=20; //in mm

	sprintf(inifilename,"%s",infilename.Data());

	const Int_t totplots=6;
	Double_t distance[totplots]={6000,9000,14500,20000,24000,28250};
	TH2F* histos[totplots];

// 	const Int_t totplots=1;
// 	Double_t distance[totplots]={28479};

	char pstitle[200];

	printf("Using data file %s\n",inifilename);
	gROOT->SetStyle("Plain");
	gStyle->SetPalette(1);
	gStyle->SetOptStat(kFALSE);
	TFile *file0 = TFile::Open(inifilename);
	TCanvas *canvas = new TCanvas("canvas","Moller plots",800,800);

	char title[100];
	char weight[200];
	TGraph* scat[totplots];
	TLine line;
	TArc arc;
	arc.SetFillStyle(4000);
	arc.SetNoEdges(kFALSE);
	Double_t pi=3.14159265358979323846;
	Double_t length=600, openingangle=pi/14.;
	Double_t ymin, ymax, xmin, xmax;
	TTree *geant = (TTree*)file0->Get("geant");
	geant->SetMarkerStyle(20);
	geant->SetMarkerSize(0.15);
	for (int i=0; i<totplots; i++) {
		char drawstring[200];
		sprintf(drawstring,"y/1000:x/1000:1000*acos(pz2/sqrt(px2*px2+py2*py2+pz2*pz2))>>htemp%i",i);
		sprintf(weight,"(z>(%.0f)&&z<(%.0f)&&type==0)",distance[i]-deltaZwindowsize,distance[i]+deltaZwindowsize);
		if (debug) printf("%s\n%s\n",drawstring,weight);
		canvas->cd(i+1);
		Int_t num = geant->Draw(drawstring,weight,"colz");
		histos[i]=(TH2F*)gPad->GetPrimitive(Form("htemp%i",i));
		if (num<=0) {
			printf("\n%.0f doesn't seem to be a good choice!\n\n",distance[i]);
		} else {
			if (debug) printf("%i successful events\n",num);
			gPad->SetRightMargin(0.12);
			gPad->SetLeftMargin(0.12);
			gPad->SetTopMargin(0.12);
			gPad->SetBottomMargin(0.12);
			canvas->Update();
			ymin=gPad->GetUymin();
			ymax=gPad->GetUymax();
			xmin=gPad->GetUxmin();
			xmax=gPad->GetUxmax();
			printf("xmin %f, xmax %f, ymin %f, ymax %f, deltax %f, deltay %f\n",xmin,xmax,ymin,ymax,xmax-xmin,ymax-ymin);
		}
		sprintf(title,"z=%.2f m;x  (m);y   (m)", distance[i]/1000.);
		histos[i]->SetTitle(title);
		canvas->Update();
		canvas->Print(Form("%s_%.2fm.png",outfilenamestem.Data(),distance[i]/1000.));
	}	
}

void csv(TString input="tmva.csvoutput.txt", TString par1="par2", TString par2="par3", TString par3="", TString value="eventEffScaled_5") {
  std::cout << "Usage:" << std::endl
            << ".x scripts/csv.C    with default arguments" << std::endl
            << ".x scripts/csv.C(filename, par1, par2, value)" << std::endl
            << std::endl
            << "  Optional arguments:" << std::endl
            << "    filename        path to CSV file" << std::endl
            << "    par1            name of X-parameter branch" << std::endl
            << "    par2            name of Y-parameter branch (if empty, efficiency is drawn as a function of par1)" << std::endl
            << "    value           name of result (efficiency) branch" << std::endl
            << std::endl;

  TTree *tree = new TTree("data", "data");
  tree->ReadFile(input);

  gStyle->SetPalette(1);
  gStyle->SetPadRightMargin(0.14);

  TCanvas *canvas = new TCanvas("csvoutput", "CSV Output", 1200, 900);

  tree->SetMarkerStyle(kFullDotMedium);
  tree->SetMarkerColor(kRed);
  if(par2.Length() > 0) {
    //tree->Draw(Form("%s:%s", par2.Data(), par1.Data()));
    if(par3.Length() > 0)
      tree->Draw(Form("%s:%s:%s:%s", par1.Data(), par2.Data(), par3.Data(), value.Data()), "", "COLZ"); //, "", "Z");
    else
      tree->Draw(Form("%s:%s:%s", par2.Data(), par1.Data(), value.Data()), "", "COLZ"); //, "", "Z");

    TH1 *histo = tree->GetHistogram();
    if(!histo)
      return;

    histo->SetTitle(Form("%s with different classifier parameters", value.Data()));
    histo->GetXaxis()->SetTitle(Form("Classifier parameter %s", par1.Data()));
    histo->GetYaxis()->SetTitle(Form("Classifier parameter %s", par2.Data()));
    if(par3.Length() > 0)
      histo->GetZaxis()->SetTitle(Form("Classifier parameter %s", par3.Data()));
    else
      histo->GetZaxis()->SetTitle("");

    if(par3.Length() == 0) {
      float x = 0;
      float y = 0;
      float val = 0;
      double maxVal = tree->GetMaximum(value);
      double minVal = tree->GetMinimum(value);

      tree->SetBranchAddress(par1, &x);
      tree->SetBranchAddress(par2, &y);
      tree->SetBranchAddress(value, &val);
      TLatex l;
      l.SetTextSize(0.03);
    
      Long64_t nentries = tree->GetEntries();
      for(Long64_t entry=0; entry < nentries; ++entry) {
        tree->GetEntry(entry);
    
        l.SetTextColor(textColor(val, maxVal, minVal));
        l.DrawLatex(x, y, Form("%.3f", val*100));
      }
    }
  }
  else {
    tree->Draw(Form("%s:%s", value.Data(), par1.Data()));
    TH1 *histo = tree->GetHistogram();
    if(!histo) 
      return;
    histo->SetTitle(Form("%s with different classifier parameters", value.Data()));
    histo->GetXaxis()->SetTitle(Form("Classifier parameter %s", par1.Data()));
    histo->GetYaxis()->SetTitle(value);
  }
}

void Loop::myLoop()
{
  // Correction factors:
  double cA[8] = {0.8567,0.8567,0.8700,0.8610,0.8567,0.8550,0.8630,0.8567};
  double cOff[8] = {0.034,0.049,0.049,0.037,0.045,0.049,0.057,0.057};
  
  // Dummy correction factors:
  // double cA[8] = {1,1,1,1,1,1,1,1};
  // double cOff[8] = {0,0,0,0,0,0,0,0};


  TFile *outfile = new TFile("output.root","recreate"); 

  TTree *nt = new TTree("data","ADC readout data");  

  adcCorrectedData myData;
  myData.nch=8;
  nt->Branch("nch",&myData.nch,"nch/I");
  nt->Branch("run",&myData.run,"run/I");
  nt->Branch("evt",&myData.evt,"evt/I");
  nt->Branch("evtSecond",&myData.evtSecond,"evtSecond/I");
  nt->Branch("evtTime",&myData.evtTime,"evtTime/I");
  nt->Branch("pedestal",myData.pedestal,"nch/I");
  nt->SetMarkerStyle(20);

  for (int i=0;i<8;i++){
    nt->Branch(Form("dSize%d",i),&myData.dataSize[i],Form("dSize%d/I",i));
    nt->Branch(Form("data%d",i),myData.data[i],Form("data%d[dSize%d]/F",i,i));
    nt->Branch(Form("time%d",i),myData.time[i],Form("time%d[dSize%d]/F",i,i)); 
  }

  nt->Branch("min",myData.min,"min[8]/F");
  nt->Branch("mintime",myData.mintime,"mintime[8]/F");
  nt->Branch("mintime08",myData.mintime08,"mintime08[8]/F");
  nt->Branch("mintime02",myData.mintime02,"mintime02[8]/F");
  nt->Branch("mintime04",myData.mintime04,"mintime04[8]/F");
  nt->Branch("mintime0802",myData.mintime0802,"mintime0802[8]/F");
  nt->Branch("mintime0804",myData.mintime0804,"mintime0804[8]/F");
  //nt->Branch("minSB",myData.minSB,"minSB[8]/F");
  //nt->Branch("minSBtime",myData.minSBtime,"minSBtime[8]/F");

  char ch;
  unsigned long dSize, dummy, trgTime;
  unsigned long nbr=0, evt=0;

   Long64_t nentries = fChain->GetEntriesFast();
   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries;jentry++) {

      if (jentry % 1000 == 0 ) cout <<jentry<<" / "<<nentries<<" "<<setprecision(2)<<(double)jentry/nentries*100<<"% \r";

      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;

      Float_t min[8];
      Float_t mintime[8];
      Float_t minSB[8];
      Float_t minSBtime[8];
       
      for (int i=0;i<8;i++) {
         myData.min[i]=10e-10;
//         myData.minSB[i]=10e-10;
         myData.mintime[i]=-100;
//         myData.minSBtime[i]=-100;
         myData.dataSize[i]=0;  
      }  

      nbr=0;
      int flag=0;
      int minIt[8];
    
      for (int i=0;i<8;i++) {
         myData.dataSize[i] = 0;
	 myData.pedestal[i] = pedestal[i];
         minIt[i]=-1;
	 
         for (int j=0;j<data[i].size();j++)
         {
            double pulseHeight = data[i][j]-pedestal[i]; //(((double)data[i][j])/255.-0.5 - cOff[i])/cA[i];
            double mytime = time[i][j];

	    myData.data[i][myData.dataSize[i]] = pulseHeight;
	    myData.time[i][myData.dataSize[i]] = mytime;
	    myData.dataSize[i]++;

  	    // fill signal region
  
            if (pulseHeight<myData.min[i]) { 
   	           myData.min[i]=pulseHeight;
	           myData.mintime[i]=mytime;
		   minIt[i]=myData.dataSize[i]-1;
	    }

            /*
            // fill sideband region
            if (pulseHeight<myData.minSB[i]) { 
               if (i!=0) {   
	          if (fabs(mytime-myData.minSBtime[0])>400&&fabs(mytime-myData.minSBtime[0])<500){ 
   	             myData.minSB[i]=pulseHeight;
//.........这里部分代码省略.........

void plot2(){

  TH2F *hframe = new TH2F("hframe","Edep_tot vs. track length for various momenta",50,0,50,1000,0.5,180.0);
  hframe->Draw(); //you can set the axis att via hframe->GetXaxis()..
  gPad->SetLogx(0);
  hframe->GetXaxis()->SetTitle("Track Length [mm]");
  hframe->GetYaxis()->SetTitle("Total Energy [MeV]");
  TLatex Tl; Tl.SetTextFont(43); Tl.SetTextSize(20);
  
  TTree *MyTree = new TTree("MyTree", "MyTree");
  MyTree->ReadFile("Rout_proton65mevOnSolenoid3LIV.txt","track:edep");
  MyTree->SetMarkerStyle(kCircle);
  MyTree->SetMarkerColorAlpha(kRed, 0.3);
  MyTree->Draw("edep:track","","same");
  TTree *MyTreeB = new TTree("MyTreeB", "MyTreeB");
  MyTreeB->ReadFile("Rout_proton65mevOnSolenoid3STD.txt","trackB:edepB");
  MyTreeB->SetMarkerStyle(kFullCircle);
  MyTreeB->SetMarkerColorAlpha(kRed, 0.3);
  MyTreeB->Draw("edepB:trackB","","same");
  Tl.DrawLatex(5, 25,   "p_{s} = 65 MeV/c @ vtx=0.0");
 
  
  TTree *MyTree1 = new TTree("MyTree1", "MyTree1");
  MyTree1->ReadFile("Rout_proton72mevOnSolenoid3LIV.txt","track1:edep1");
  MyTree1->SetMarkerStyle(kOpenSquare);
  MyTree1->SetMarkerColorAlpha(kBlue, 0.3);
  MyTree1->Draw("edep1:track1","","same");
  TTree *MyTree1B = new TTree("MyTree1B", "MyTree1B");
  MyTree1B->ReadFile("Rout_proton72mevOnSolenoid3STD.txt","track1B:edep1B");
  MyTree1B->SetMarkerStyle(kFullSquare);
  MyTree1B->SetMarkerColorAlpha(kBlue, 0.3);
  MyTree1B->Draw("edep1B:track1B","","same");
  Tl.DrawLatex(10, 43,   "p_{s} = 72 MeV/c");

  TTree *MyTree2 = new TTree("MyTree2", "MyTree2");
  MyTree2->ReadFile("Rout_proton90mevOnSolenoid3LIV.txt","track2:edep2");
  MyTree2->SetMarkerStyle(kOpenTriangleUp);
  MyTree2->SetMarkerColorAlpha(kGreen, 0.3);
  MyTree2->Draw("edep2:track2","","same");
  TTree *MyTree2B = new TTree("MyTree2B", "MyTree2B");
  MyTree2B->ReadFile("Rout_proton90mevOnSolenoid3STD.txt","track2B:edep2B");
  MyTree2B->SetMarkerStyle(kFullTriangleUp);
  MyTree2B->SetMarkerColorAlpha(kGreen, 0.3);
  MyTree2B->Draw("edep2B:track2B","","same");
  Tl.DrawLatex(15, 70,   "p_{s} = 90 MeV/c");

  TTree *MyTree3 = new TTree("MyTree3", "MyTree3");
  MyTree3->ReadFile("Rout_proton120mevOnSolenoid3LIV.txt","track3:edep3");
  MyTree3->SetMarkerStyle(kOpenDiamond);
  MyTree3->SetMarkerColorAlpha(kMagenta, 0.3);
  MyTree3->Draw("edep3:track3","","same");
  TTree *MyTree3B = new TTree("MyTree3B", "MyTree3B");
  MyTree3B->ReadFile("Rout_proton120mevOnSolenoid3STD.txt","track3B:edep3B");
  MyTree3B->SetMarkerStyle(kFullDiamond);
  MyTree3B->SetMarkerColorAlpha(kMagenta, 0.3);
  MyTree3B->Draw("edep3B:track3B","","same");
  Tl.DrawLatex(20, 105,   "p_{s} = 120 MeV/c");

  TTree *MyTree4 = new TTree("MyTree4", "MyTree4");
  MyTree4->ReadFile("Rout_proton180mevOnSolenoid3LIV.txt","track4:edep4");
  MyTree4->SetMarkerStyle(kOpenTriangleDown);
  MyTree4->SetMarkerColorAlpha(kBlack, 0.3);
  MyTree4->Draw("edep4:track4","","same");
  TTree *MyTree4B = new TTree("MyTree4B", "MyTree4B");
  MyTree4B->ReadFile("Rout_proton180mevOnSolenoid3STD.txt","track4B:edep4B");
  MyTree4B->SetMarkerStyle(kFullTriangleDown);
  MyTree4B->SetMarkerColorAlpha(kBlack, 0.3);
  MyTree4B->Draw("edep4B:track4B","","same");
  Tl.DrawLatex(25, 170,   "p_{s} = 180 MeV/c");

  Tl.DrawLatex(10, 142,   "open:LIV, full:STD,(0.1% lower)");

}

//.........这里部分代码省略.........
    fcint = fc.GetInterval();  // that was easy.

    RooFitResult* fit = modelWithConstraints->fitTo(*data, Save(true));

    // Third, use a Calculator based on Markov Chain monte carlo
    // Before configuring the calculator, let's make a ProposalFunction
    // that will achieve a high acceptance rate
    ProposalHelper ph;
    ph.SetVariables((RooArgSet&)fit->floatParsFinal());
    ph.SetCovMatrix(fit->covarianceMatrix());
    ph.SetUpdateProposalParameters(true);
    ph.SetCacheSize(100);
    ProposalFunction* pdfProp = ph.GetProposalFunction();  // that was easy

    MCMCCalculator mc(*data, modelConfig);
    mc.SetNumIters(20000); // steps to propose in the chain
    mc.SetTestSize(.05); // 95% CL
    mc.SetNumBurnInSteps(40); // ignore first N steps in chain as "burn in"
    mc.SetProposalFunction(*pdfProp);
    mc.SetLeftSideTailFraction(0.5);  // find a "central" interval
    MCMCInterval* mcInt = (MCMCInterval*)mc.GetInterval();  // that was easy


    // Get Lower and Upper limits from Profile Calculator
    cout << "Profile lower limit on s = " << ((LikelihoodInterval*) lrinterval)->LowerLimit(*s) << endl;
    cout << "Profile upper limit on s = " << ((LikelihoodInterval*) lrinterval)->UpperLimit(*s) << endl;

    // Get Lower and Upper limits from FeldmanCousins with profile construction
    if (fcint != NULL) {
        double fcul = ((PointSetInterval*) fcint)->UpperLimit(*s);
        double fcll = ((PointSetInterval*) fcint)->LowerLimit(*s);
        cout << "FC lower limit on s = " << fcll << endl;
        cout << "FC upper limit on s = " << fcul << endl;
        TLine* fcllLine = new TLine(fcll, 0, fcll, 1);
        TLine* fculLine = new TLine(fcul, 0, fcul, 1);
        fcllLine->SetLineColor(kRed);
        fculLine->SetLineColor(kRed);
        fcllLine->Draw("same");
        fculLine->Draw("same");
        dataCanvas->Update();
    }

    // Plot MCMC interval and print some statistics
    MCMCIntervalPlot mcPlot(*mcInt);
    mcPlot.SetLineColor(kMagenta);
    mcPlot.SetLineWidth(2);
    mcPlot.Draw("same");

    double mcul = mcInt->UpperLimit(*s);
    double mcll = mcInt->LowerLimit(*s);
    cout << "MCMC lower limit on s = " << mcll << endl;
    cout << "MCMC upper limit on s = " << mcul << endl;
    cout << "MCMC Actual confidence level: "
         << mcInt->GetActualConfidenceLevel() << endl;

    // 3-d plot of the parameter points
    dataCanvas->cd(2);
    // also plot the points in the markov chain
    RooDataSet * chainData = mcInt->GetChainAsDataSet();

    assert(chainData);
    std::cout << "plotting the chain data - nentries = " << chainData->numEntries() << std::endl;
    TTree* chain =  RooStats::GetAsTTree("chainTreeData","chainTreeData",*chainData);
    assert(chain);
    chain->SetMarkerStyle(6);
    chain->SetMarkerColor(kRed);

    chain->Draw("s:ratioSigEff:ratioBkgEff","nll_MarkovChain_local_","box"); // 3-d box proportional to posterior

    // the points used in the profile construction
    RooDataSet * parScanData = (RooDataSet*) fc.GetPointsToScan();
    assert(parScanData);
    std::cout << "plotting the scanned points used in the frequentist construction - npoints = " << parScanData->numEntries() << std::endl;
    // getting the tree and drawing it -crashes (very strange....);
    // TTree* parameterScan =  RooStats::GetAsTTree("parScanTreeData","parScanTreeData",*parScanData);
    // assert(parameterScan);
    // parameterScan->Draw("s:ratioSigEff:ratioBkgEff","","goff");
    TGraph2D *gr = new TGraph2D(parScanData->numEntries());
    for (int ievt = 0; ievt < parScanData->numEntries(); ++ievt) {
        const RooArgSet * evt = parScanData->get(ievt);
        double x = evt->getRealValue("ratioBkgEff");
        double y = evt->getRealValue("ratioSigEff");
        double z = evt->getRealValue("s");
        gr->SetPoint(ievt, x,y,z);
        // std::cout << ievt << "  " << x << "  " << y << "  " << z << std::endl;
    }
    gr->SetMarkerStyle(24);
    gr->Draw("P SAME");


    delete wspace;
    delete lrinterval;
    delete mcInt;
    delete fcint;
    delete data;

    // print timing info
    t.Stop();
    t.Print();
}