void fitH(TH1F *h, TF1 *g)
{
  Double_t rmax, rmin;

  Int_t i = h->GetMaximumBin();
  while(h->GetBinContent(i) > 0.5*h->GetMaximum()) {
    rmax = h->GetBinCenter(i+1);
    i++;
  }
  
  i = h->GetMaximumBin();
  while(h->GetBinContent(i) > 0.5*h->GetMaximum()) {
    rmin = h->GetBinLowEdge(i-1);
    i--;
  }

  //cout << "########" << endl << rmin << " "  << rmax << endl;
  g->SetParLimits(1,rmin,rmax);
  g->SetRange(rmin,rmax);
  h->Fit("g","R");
  //ofile << g->GetParameter(1) << endl; 
}

// main macro function
// ===================
void drawData()
{
  Int_t par1,par2,par3,par4;

  par3=1;

  TString inFilename("/home/muzalevsky/work/exp1803/data/exp1804/h5_14/out40cal.root");
  TString inTreeName("tree");
  TString outPdfDir("pdf");

  // Open the input file
  // ========================================================================
  TFile *f1 = new TFile(inFilename, "READ");
  if (f1->IsZombie()) {
    cerr << "Could not open file '" << inFilename << "'" << endl;
    return;
  } else {
    cerr << "Input file '" << inFilename << "' opened." << endl;
  }
  // ========================================================================

  // Get the input tree
  // ========================================================================
  TTree *t1 = (TTree*)f1->Get(inTreeName);
  if (t1 == NULL) {
    cerr << "Could not find tree '" << inTreeName << "' in file '" << inFilename << "'" << endl;
    return;
  } else {
    cerr << "Input tree '" << inTreeName << "' found." << endl;
  }
  ////t1->SetMakeClass(1);
  // ========================================================================

/*
  if (par1) {

    cerr << "Processing 'par1' part." << endl;

    TFile *f = new TFile("/home/muzalevsky/work/exp1803/data/exp1804/h5_14/out.root");
    TTree *t = (TTree*)f->Get("tree");
    TCanvas *c1 = new TCanvas();
    c1->Divide(2,2);
    c1->cd(1);
    c1_1->SetLogz();
//    t->Draw("F5[0]+F3[0] : tF5[0]-tF3[0] >> (500,100,700,500,100,4000)","","col");
    t->Draw("F5[]+F3[] : tF5[]-tF3[] >> (500,100,700,500,100,4000)","","col");

    c1->cd(2);
    c1_2->SetLogz();
    t->Draw("SQX_R : CsI_R","SQX_R>1.15 && CsI_R>0","col"); 

    c1->cd(3);
    //c1_3->SetLogz();
//    t->Draw("F5[0]+F3[0] : tF5[0]-tF3[0]","(F5[0]+F3[0])>300 && (F5[0]+F3[0])<700 && (tF5[0]-tF3[0])>540 && (tF5[0]-tF3[0])<640","col"); 
    t->Draw("F5[]+F3[] : tF5[]-tF3[]","(F5[]+F3[])>300 && (F5[]+F3[])<700 && (tF5[]-tF3[])>540 && (tF5[]-tF3[])<640",""); 

    //c1->cd(4);
    //t->Draw("SQX_R:CsI_R","SQX_R>1.15 && CsI_R>0 && (F5[]+F3[])>300 && (F5[]+F3[])<700 && (tF5[]-tF3[])>540 && (tF5[]-tF3[])<640",""); 

  } // end of 'par1' part

  if (par2) {

    cerr << "Processing 'par2' part." << endl;

    TCanvas *c2 = new TCanvas();
    c2->Divide(2,2);
    c1_1->SetLogz();
    c2->cd(1);
    TH2F *h1 = new TH2F("h1","TOF",500,100.,1000.,500,0.,7000.);
    t->Draw("NeEvent.F5[0]+NeEvent.F3[0] : NeEvent.tF5[0]-NeEvent.tF3[0] >> h1","","col");
    c2->Update();

    c2->cd(2);
    TH2F *h2 = new TH2F("h2","TOF cut",100,600.,750.,100,300.,800.);
    c1_2->SetLogz();
    t->Draw("NeEvent.F5[0]+NeEvent.F3[0]:NeEvent.tF5[0]-NeEvent.tF3[0] >> h2","","col");
    c2->Update();

    c2->cd(3);
    t->Draw("NeEvent.SQX_R[0]","NeEvent.SQX_R[0]>0 && NeEvent.SQX_R[0]<2000","");
    c2->Update();

    c2->cd(4);
    t->Draw("NeEvent.SQX_R[]","NeEvent.tF5[]-NeEvent.tF3[] > 600 && NeEvent.tF5[]-NeEvent.tF3[] < 750 && NeEvent.F5[]+NeEvent.F3[]>300 && NeEvent.F5[]+NeEvent.F3[]<800 && NeEvent.SQX_R[]>0","");
    
    TCanvas *c1 = new TCanvas();
    c1->Divide(2,2);
    c1->cd(1);
    t->Draw("NeEvent.SQX_R[]:NeEvent.CsI_R[]","","");
    c1->cd(2);
    t->Draw("NeEvent.SQX_R[0]:NeEvent.CsI_R[0]","","");
    c1->cd(3);
    t->Draw("NeEvent.SQX_R[0]","NeEvent.SQX_R[0]>0 && NeEvent.SQX_R[0]<2000","");
    c1->cd(4);
    t->Draw("NeEvent.SQX_R[0]:NeEvent.CsI_R[0]","NeEvent.SQX_R[0]>200 && NeEvent.CsI_R[0]>100","");

  } // end of 'par2' part
*/
  if (par3) {

    cerr << "Processing 'par3' part." << endl;

    TH1I *hmulX = new TH1I("hmulX","x multiplicity",4,1,5);
    TH1I *hmulY = new TH1I("hmulY","y multiplicity",4,1,5);

    TH2F *h1;
    TH2F *h2;     
   

    TFile *f1 = new TFile("/home/muzalevsky/work/exp1803/data/exp1804/h5_14/out40cal.root");
    TTree *t1 = (TTree*)f1->Get("tree");
    TCanvas *c2 = new TCanvas("c2","Si amp^times",1800,1000);
    c2->Divide(2,2);
    c2->cd(1);
    gPad->SetLogz();
    //t1->Draw("SQY_L[4]:tSQY_L[4]>>htemp1(400,1000,1400,100,0,90)","tSQY_L[4]>0 && SQY_L[4]>0","box",5000000,0);
    t1->Draw("SQX_L[8]:tSQX_L[8]>>htemp1(400,200,600,100,0,90)","tSQX_L[8]>0 && SQX_L[8]>0","box",20000000,0);
    h1 = (TH2F*)gPad->GetPrimitive("htemp1");
    //cout << h1 << endl;
    //c2->cd(3);
    h1->Draw();    


    c2->cd(2);
    gPad->SetLogz();
    //c2_2->SetLogz();
    //t1->Draw("SQY_L[4]:tSQY_L[4]>>htemp2(400,1000,1400,100,0,90)","tSQY_L[4]>0 && SQY_L[4]>0 && multY_L>1","box",5000000,0);
    t1->Draw("SQX_L[8]:tSQX_L[8]>>htemp2(400,200,600,100,0,90)","tSQX_L[8]>0 && SQX_L[8]>0 && multY_L>1","box",20000000,0);
    h2 = (TH2F*)gPad->GetPrimitive("htemp2");
    h2->Draw();   

    //h2->Draw();
    c2->cd(3);
    gPad->SetLogz();
    //t1->Draw("SQY_L[]:tSQY_L[]>>htemp3(400,1000,1400,100,0,90)","tSQY_L[]>0 && SQY_L[]>0","col",5000000,0);
    t1->Draw("SQX_L[]:tSQX_L[]>>htemp3(400,200,600,100,0,90)","tSQX_L[]>0 && SQX_L[]>0","col",500000,0);
    //t1->Draw("SQY_L[12]:tSQY_L[12]","multY_L==2","",4000000,0);

    c2->cd(4);
    gPad->SetLogz();
    //t1->Draw("SQY_L[]:tSQY_L[]>>htemp4(400,1000,1400,100,0,90)","multY_L==1 && tSQY_L[]>0 && SQY_L[]>0","col",5000000,0);
    t1->Draw("SQX_L[]:tSQX_L[]>>htemp4(400,200,600,100,0,90)","multY_L==1 && tSQX_L[]>0 && SQX_L[]>0","col",500000,0);
    
  } // end of 'par3' part



  /*if (par4) {

    cerr << "Processing 'par4' part." << endl;

    // Open output text file
    // ========================================================================
   // TString outTxtFilename("/home/muzalevsky/work/exp1803/data/exp1804/h5_14/positionstest.txt");
    TString outTxtFilename("tParY_L.clb");
    ofstream outcalfile;
    outcalfile.open(outTxtFilename.Data());
    if (!outcalfile.is_open()) {
      cerr << "Could not open output text file'" << outTxtFilename << "'" << endl;
      return;
    } else {
      cerr << "Opened output text file '" << outTxtFilename << "'" << endl;
    }

	    
    // ========================================================================

    //TH1::AddDirectory(kFALSE);

    // Create a Gaussian function
    TF1* g = new TF1("g", "gaus", 1000., 1200.);
    //g->SetParLimits(0, 1., 1e+7);
    g->SetParLimits(1, 1000., 1200.);
    g->SetParLimits(2, 1., 50.);
  
    Int_t nhists = 16;
    Int_t ncanvases = nhists/4;
    Double_t shift,zeroLevel;

    TH1F *hX[32];
    TH1F *hY[16];
    TH2F *h2D[32];  

    TH1F* h = new TH1F("h","temp hist", 140,1020.,1160.);
    TH2F *h2 = new TH2F("h2","temp 2d hist",140,1020.,1160.,400,0.,90.);

    TString hname,h2name;
    //for(Int_t i=0;i<nhists;i++) {
    //  TString hname;
   //   hname.Form("h%d",i+1);
    //  hX[i] = new TH1F(hname.Data(),"title",1200,0,1200);
   // }

    outcalfile << 2 << endl << nhists<< endl;

    // Create canvases
    // ========================================================================
    TCanvas* c[4]; // works only if ncanvases==8
    for(Int_t i=0; i<4; i++) {
      TString cName;
      cName.Form("c_%d",i);
      c[i] = new TCanvas(cName.Data(),"calibrated spectra",1000,1000);
      c[i]->Divide(2,2);
    }

    TCanvas* c2[4]; // works only if ncanvases==8
    for(Int_t i=0; i<4; i++) {
      TString cName;
      cName.Form("c2_%d",i);
      c2[i] = new TCanvas(cName.Data(),"2-dim calibrated spectra",1000,1000);
      c2[i]->Divide(2,2);
    }
    cerr << "Canvases created." << endl;
    // ========================================================================

    for (Int_t i=0;i<nhists;i++)
    {

      cerr << "i=" << i << endl;


      hname.Form("h%d",i);

      TString cut;
      cut.Form("SQY_L[%d]>14",i);
      //cut.Form("tSQX_L[%d]>1000 && tSQX_L[%d]<1600",i,i,i);

      TString vary1;
      //vary1.Form("tSQX_L[%d]>>(200,1000,1200)",i);
      vary1.Form("tSQY_L[%d]>>h%d(200,1000,1200)",i,i);

      Int_t count = i/4;
      Int_t nPad = (i%4)+1;

      c[count]->cd(nPad);
      t1->Draw(vary1.Data(),cut.Data(),"",10000000,0);

      hX[i] = (TH1F*)gPad->GetPrimitive(hname.Data());
      hX[i]->SetTitle(hname.Data());
      fitH(hX[i],g);
      if(i==0) {
        zeroLevel = g->GetParameter(1);
        shift = 0;
      } 
      shift = zeroLevel - g->GetParameter(1);
	    outcalfile
	    << setprecision(4) << shift*0.3 << "\t"
	    << setprecision(4) << "0.3" << "\t\t" //E=a+bN
	    << "0" << endl;

    
      c[count]->Update();

      c2[count]->cd(nPad);      

      h2name.Form("h2_%d",i);
      TString cut1;      
      cut1.Form("tSQX_L[%d]>900 && tSQX_L[%d]<1600",i,i);

      TString vary;
      vary.Form("SQY_L[%d]:tSQY_L[%d]>>h2_%d(400,1000.,1400.,100,0.,90)",i,i,i);
      gPad->SetLogz();
      t1->Draw(vary.Data(),"","col",100000,0);

     
      h2D[i] = (TH2F*)gPad->GetPrimitive(h2name.Data());
      TString title;
      title.Form("2-dim %d",i);
      h2D[i]->SetTitle(title.Data());
      //h2D[i]->Draw();

      c2[count]->Update();

    }

    // Export canvases into output PDF files
    // ========================================================================
    for (Int_t iCanv=0; iCanv<ncanvases; iCanv++) {
      TString outFilename;
      outFilename.Form("%s/canv_%d.pdf", outPdfDir.Data(), iCanv);
      c[iCanv]->SaveAs(outFilename);
      outFilename.Form("%s/canv2_%d.pdf", outPdfDir.Data(), iCanv);
      c2[iCanv]->SaveAs(outFilename);
    }
    cerr << "Exported canvases into output PDF files." << endl;
    // ========================================================================


  } // end of 'par4' part*/

}