C++ TObjArray::Write方法代码示例

int main(int argc, char* argv[]){
  if (argc != 3 && argc != 1) return 0;

  TH1::AddDirectory(false);
  gROOT->ProcessLine(".x /user2/sfarry/WJet/root/lhcbStyle.C");

  if (argc == 1){
    TFile* file0 = TFile::Open("root://hepgrid11.ph.liv.ac.uk///dpm/ph.liv.ac.uk/home/lhcb/electroweak/zmumu.pythia.root");
    TTree* t = (TTree*)file0->Get("Zjet");
    TObjArray* hists = get_mass_v_vars(t);
    TFile* f = new TFile("output.root","RECREATE");
    hists->Write();
    f->Close();

  }
  else if (argc == 3){
    char* input = argv[1];
    char* output = argv[2];
    TFile* f = TFile::Open(input, "READ");
    TTree* t = (TTree*)f->Get("Zjet");
  
    TObjArray *hists = 0;
    if ( t ) {
      hists = get_mass_v_vars(t);
    }
    if (hists){
      TFile* g = TFile::Open(output, "RECREATE");
      hists->Write();
      g->Close();
    }
  }
  return 0;
}

/* *************************************************************************** *
 *                                                                             *
 * Format:                                                                     *
 *                                                                             *
 * Save into a format to be used in main code                                  *
 *                                                                             *
 * *************************************************************************** */
void GeomAcc::Format(TString saveName, TString pol){

  TString plotName;
  if(pol == "NoPol/") plotName = "GeoEff";
  else if (pol == "Long/") plotName = "GeoLongEff";
  else if (pol == "Trans/") plotName = "GeoTransEff"; 

  TFile *f1 = TFile::Open("Output/" + saveName);
  TH2D *GeoEff2 = (TH2D*) f1->Get(plotName);
  GeoEff2->SetName("Geometric Acceptance");

  //
  //hack to use existing description plot
  //

  TFile *f2 = TFile::Open("Description/" + saveName);
  TH1 * Descrip = (TH1*) f2->Get("Description");

  TObjArray * histList = new TObjArray(0);
  histList->Add(GeoEff2);
  histList->Add(Descrip);

  TFile* saveFile = new TFile("Efficiency/"+ pol+ saveName, "recreate");
  histList->Write();
  saveFile->Close();
 

  return;
}

void DrawGlobal(TObjArray List,TFile *in, TFile *out){
  TObjArray List;

  TCanvas *MeanTk = new TCanvas("MeanTk","MeanTk",10,10,900,500);

  TIFTree->Project("h0","Tracks.MeanTrack:number","Clusters.entries_all>2000");
  
  h0->SetTitle("Number of Events with at least on Track/Total number of Events");
  h0->GetXaxis()->SetTitle("Run Number");
  h0->GetXaxis()->CenterTitle(1);
  h0->GetYaxis()->SetTitle("Fraction of Events with nTracks non zero");
  h0->GetYaxis()->CenterTitle(1); 
  h0->SetMarkerStyle(20);
  h0->SetStats(0);
  h0->Draw();
  List.Add(MeanTk);

  TCanvas *Events = new TCanvas("Events","Events",10,10,900,500);

  TIFTree->Project("h3","Clusters.entries_all:number","Clusters.entries_all>2000");
  
  h3->SetTitle("Total Number of Events");
  h3->GetXaxis()->SetTitle("Run Number");
  h3->GetXaxis()->CenterTitle(1);
  h3->GetYaxis()->SetTitle("Events");
  h3->GetYaxis()->CenterTitle(1); 
  h3->SetMarkerStyle(20);
  h3->SetStats(0);
  h3->Draw();
  List.Add(Events);

  TCanvas *nTracks = new TCanvas("nTracks","nTracks",10,10,900,500);

  TIFTree->Project("h","Tracks.mean:number","Clusters.entries_all>2000");
  
  h->SetTitle("Mean Track number per Event");
  h->GetXaxis()->SetTitle("Run Number");
  h->GetXaxis()->CenterTitle(1);
  h->GetYaxis()->SetTitle("Mean Track number per Event");
  h->GetYaxis()->CenterTitle(1); 
  h->SetMarkerStyle(20);
  h->SetStats(0);
  h->Draw();
  List.Add(nTracks);
  
  TCanvas *nRecHits= new TCanvas("nRecHits","nRecHits",10,10,900,500);  
  TIFTree->Project("h1","RecHits.mean:number","Clusters.entries_all>2000");
  
  h1->SetTitle("Mean RecHits Number per Event");
  h1->GetXaxis()->SetTitle("Run Number");
  h1->GetXaxis()->CenterTitle(1);
  h1->GetYaxis()->SetTitle("Mean RecHits Number per Event");
  h1->GetYaxis()->CenterTitle(1); 
  h1->SetMarkerStyle(20);
  h1->SetStats(0);
  h1->Draw();
  List.Add(nRecHits);

  TCanvas *nCluster= new TCanvas("nClusters","nClusters",10,10,900,500);
  TIFTree->Project("h2","Clusters.mean_corr:number","Clusters.entries_all>2000");
  
  h2->SetTitle("Mean Number of Cluster per Event");
  h2->GetXaxis()->SetTitle("Run Number");
  h2->GetXaxis()->CenterTitle(1);
  h2->GetYaxis()->SetTitle("Mean Number of Cluster per Event");
  h2->GetYaxis()->CenterTitle(1); 
  h2->SetMarkerStyle(20);
  h2->SetStats(0);
  h2->Draw();
  List.Add(nClusters);

 out->cd();
 List.Write();
 // out->Close();

 in->cd();
}

//.........这里部分代码省略.........
  TCut cut("muminus_Refit_PT > 20000 && muplus_Refit_PT > 20000 && muplus_Refit_ETA > 2.0 && muplus_Refit_ETA < 4.5 && muminus_Refit_ETA > 2.0 && muminus_Refit_ETA < 4.5 && Z0_Refit_M > 60000 && Z0_Refit_M < 120000");

  RooRealVar m( "Z0_Refit_M", "Z0_Refit_M", 60000., 120000. );

  RooRealVar pt1("muplus_Refit_PT"   , "muplus_Refit_PT"   , 20000, 1000000);
  RooRealVar pt2("muminus_Refit_PT"  , "muminus_Refit_PT"  , 20000, 1000000);
  RooRealVar eta1("muplus_Refit_ETA" , "muplus_Refit_ETA"  , 2, 4.5);
  RooRealVar eta2("muminus_Refit_ETA", "muminus_Refit_ETA" , 2, 4.5);
  RooRealVar phi1("muplus_Refit_PHI" , "muplus_Refit_PHI"  , -TMath::Pi(), TMath::Pi());
  RooRealVar phi2("muminus_Refit_PHI", "muminus_Refit_PHI" , -TMath::Pi(), TMath::Pi());
  RooRealVar polarity("Polarity", "Polarity", -2, 2);

  RooDataSet* dsUp = new RooDataSet("dsUp", "dsUp", RooArgSet(m, pt1, pt2, eta1, eta2, phi1, phi2, polarity), RooFit::Import(*t), RooFit::Cut(cut.GetTitle()));
  RooDataSet* dsDwn = new RooDataSet("dsDwn", "dsDwn", RooArgSet(m, pt1, pt2, eta1, eta2, phi1, phi2, polarity), RooFit::Import(*u), RooFit::Cut(cut.GetTitle()));

  return get_refit_mass_v_vars(dsUp, dsDwn);
}

int main(int argc, char* argv[]){
  if (argc != 3 && argc != 1 && argc != 4) return 0;
  



  TH1::AddDirectory(false);
  gROOT->ProcessLine(".x /user2/sfarry/WJet/root/lhcbStyle.C");

  if (argc == 1){
    TFile* file0 = TFile::Open("root://hepgrid11.ph.liv.ac.uk///dpm/ph.liv.ac.uk/home/lhcb/refit/ZMuMu.Refit.2015.root");
    TTree* t = (TTree*)file0->Get("ZMuMu/DecayTree");
    TObjArray* hists = get_mass_v_vars(t);
    TObjArray* refit_hists = get_refit_mass_v_vars(t);
    TFile* f = new TFile("test.root","RECREATE");
    hists->Write();
    refit_hists->Write();
    f->Close();

  }
  else if (argc == 3){
    char* input = argv[1];
    char* output = argv[2];
    TFile* f = TFile::Open(input, "READ");
    RooDataSet* dsUp = (RooDataSet*)f->Get("dataup");
    RooDataSet* dsDown = (RooDataSet*)f->Get("datadown");
    TTree* t = (TTree*)f->Get("ZMuMu/DecayTree");
    TObjArray* branches = t->GetListOfBranches();
    bool has_refit = false;
    for (int i = 0 ; i < branches->GetEntries() ; ++i){
      TBranch* branch = (TBranch*)branches->At(i);
      string name = branch->GetName();
      //if (strncmp((name, "boson", 5) == 0){
      if (name.find("boson") == 0){
	std::string replace = name.replace(0, 5, "Z0");
	branch->SetNameTitle(replace.c_str(), replace.c_str());
      }
      if (name.find("Refit") != std::string::npos) has_refit = true;

    }
  
    TObjArray *hists = 0;
    TObjArray *refit_hists = 0;
    if ( dsUp && dsDown ) {
      hists = get_mass_v_vars(dsUp, dsDown);
      refit_hists = get_refit_mass_v_vars(dsUp, dsDown);
    }
    else if ( t ) {

//.........这里部分代码省略.........
  AliMpDEIterator deit;

  // Loop over Chambers
  for (Int_t ich = 0 ; ich < 10 ; ich++)
  {
    // Compute the ratios for DE vs run
    deit.First(ich);
  
//    while ( !deit.IsDone() )
//    {
//      TString currentDEName = Form("EffDE%dVSrun",deit.CurrentDEId());
//      gData = static_cast<TGraphAsymmErrors*>(listDEEffVSrunData->FindObject(currentDEName.Data()));
//      gSim = static_cast<TGraphAsymmErrors*>(listDEEffVSrunSim->FindObject(currentDEName.Data()));
//      
//      TString name =  Form("RatioEffDE%dVsRun",deit.CurrentDEId()); TString title = Form("DE %d data/sim tracking efficiency versus run",deit.CurrentDEId());
//      if (!gData || !gSim )
//      {
//        printf("Error readig from DE efficiency vs run list \n");
//        return;
//      }
//      //----Eff of DEs vs run
//      TGraphAsymmErrors *ratioEffDEvsRun = CreateRatioGraph(name.Data(),title.Data(),*gData,*gSim);
//      deVSrunRatios.Add(ratioEffDEvsRun);
//      
//      TGraphAsymmErrors* ratioEffDEvsRunCopy = static_cast<TGraphAsymmErrors*>(ratioEffDEvsRun->Clone());
//      
//      TString nameRatio =  Form("RatioEffDE%dVsRunAndEff",deit.CurrentDEId()); TString titleRatio = Form("Comparison Data&MC DE%d tracking efficiency versus run",deit.CurrentDEId());
//      deVSrunRatiosAndEff.Add(DrawRatio(nameRatio.Data(),titleRatio.Data(),gData,gSim,ratioEffDEvsRunCopy));
//      //-------
//
////      deVSrunRatios.Add(CreateRatioGraph(name.Data(),title.Data(),*gData,*gSim));
//      
//      deit.Next();
//    }
  
    // Compute the ratios for Ch vs DE
    TString hname3 = integrated ? "integratedEff" : "eff";
    gData = static_cast<TGraphAsymmErrors*>(listChEffVSDEData->FindObject(Form("%sCh%dVsDE",hname3.Data(),ich+1)));
    gSim = static_cast<TGraphAsymmErrors*>(listChEffVSDESim->FindObject(Form("%sCh%dVsDE",hname3.Data(),ich+1)));
    
    if (!gData || !gSim )
    {
      printf("Error reading from Chamber efficiency per DE list \n");
      return;
    }
    TString name =  Form("RatioEffCh%dVsDE",ich+1); TString title = Form("Chamber %d data/sim tracking efficiency versus DE",ich+1);
    //----Eff of CHs vs DE
    TGraphAsymmErrors *ratioEffChvsDE = CreateRatioGraph(name.Data(),title.Data(),*gData,*gSim);
    chamberVSdeRatios.Add(ratioEffChvsDE);
    
    TGraphAsymmErrors* ratioEffChvsDECopy = static_cast<TGraphAsymmErrors*>(ratioEffChvsDE->Clone());
    
    TString nameRatio =  Form("RatioEffCh%dVsDEAndEff",ich+1); TString titleRatio = Form("Comparison Data&MC Ch%d tracking efficiency versus DE",ich+1);
    chamberVSdeRatiosAndEff.Add(DrawRatio(nameRatio.Data(),titleRatio.Data(),gData,gSim,ratioEffChvsDECopy));
    //-------

    
//    chamberVSdeRatios.Add(CreateRatioGraph(name.Data(),title.Data(),*gData,*gSim));
    
  }

  //Beautify graphs
  BeautifyGraphs(globalRatios,"EffData/EffSim");
//  BeautifyGraphs(deVSrunRatios,"EffData/EffSim");
//  BeautifyGraphs(chamberVSrunRatios,"EffData/EffSim");
  BeautifyGraphs(chamberVSdeRatios,"EffData/EffSim");

//  BeautifyGraphs(globalRatiosAndEff,"EffData/EffSim");
//  BeautifyGraphs(deVSrunRatiosAndEff,"EffData/EffSim");
//  BeautifyGraphs(chamberVSrunRatiosAndEff,"EffData/EffSim");
//  BeautifyGraphs(chamberVSdeRatiosAndEff,"EffData/EffSim");

  // set bin labels
//  SetRunLabel(deVSrunRatios,irun,runs);
//  SetRunLabel(chamberVSrunRatios,irun,runs);
//  SetRunLabel(globalRatios,irun,runs,1); //Write it in such a way the number is the position on the list of the graph you want to label
//  
  // save output
  TFile* file = new TFile("EffComparison.root","update");
  
  globalRatios.Write("GlobalEffRatios", TObject::kOverwrite | TObject::kSingleKey);
//  chamberVSrunRatios.Write("ChambersEffVSrunRatios", TObject::kOverwrite | TObject::kSingleKey);
//  deVSrunRatios.Write("DEEffVSrunRatios", TObject::kOverwrite | TObject::kSingleKey);
  chamberVSdeRatios.Write("ChamberEffperDERatios", TObject::kOverwrite | TObject::kSingleKey);
  
  globalRatiosAndEff.Write("GlobalEffRatiosAndEffs", TObject::kOverwrite | TObject::kSingleKey);
//  chamberVSrunRatiosAndEff.Write("ChambersEffVSrunRatiosAndEff", TObject::kOverwrite | TObject::kSingleKey);
//  deVSrunRatiosAndEff.Write("DEEffVSrunRatiosAndEff", TObject::kOverwrite | TObject::kSingleKey);
  chamberVSdeRatiosAndEff.Write("ChamberEffperDERatiosAndEff", TObject::kOverwrite | TObject::kSingleKey);

  
   
  file->Close();
 
  fileData->Close();
  fileSim->Close();
   


}

int rdMu2bbc(
 TString fullName="ccX",
 int nEve=10000,
 Int_t nFiles  = 4,
 char* file="inp/R5112017.lis", // min-b
 //char* file="inp/R5118053.lis", // prodPP
 char* inDir   = "./",
 char* outDir   = "outBbc/"
){ 
  
  gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
  loadSharedLibraries();
  cout << " loading done " << endl;
  gSystem->Load("myTest1");  

  // create chain    
  chain = new StChain("StChain"); 

  printf("adding muDst from '%s' ....\n",file);
  // Now we add Makers to the chain...   
  muMk = new StMuDstMaker(0,0,inDir,file,"MuDst.root",nFiles);
  TChain* tree=muMk->chain(); assert(tree); 
  int nEntries=tree->GetEntries();
  printf("total eve in chain =%d\n",nEntries);

  myMk3=new StBbcVertexMaker("bbcVertex","MuDst");
  TObjArray  HList;
  myMk3->SetHList(&HList);
  myMk3->readCalib("outBbc/bbcEcalib2.dat");
  myMk3->readCalib("outBbc/bbcWcalib2.dat");
  myMk3->setTdcCalib(2.0); // cm/tdc ch

  gMessMgr->SwitchOff("D");
  gMessMgr->SwitchOn("I");
 
  chain->Init();
  chain->ls(3);
  //  return;
  int eventCounter=0;
  int stat=0;
  int t1=time(0);
  //---------------------------------------------------
  while ( stat==0 ) {// loop over events
    if(eventCounter>=nEve) break;
    eventCounter++;
    chain->Clear();
    stat = chain->Make();

    if(eventCounter%1000!=0)continue;

    printf("\n\n ====================%d  processing  ==============\n", eventCounter);

  }
  printf("sorting done, nEve=%d of %d\n",nEve, nEntries);
  int t2=time(0);
  float rate=1.*eventCounter/(t2-t1);
  float nMnts=(t2-t1)/60.;
  printf("sorting done %d of   nEve=%d, CPU rate=%.1f Hz, tot %.1f minutes\n",eventCounter,nEntries,rate,nMnts);

   chain->Finish();
   //   HList.ls();
   fullName+=".hist.root";
   TFile f( outDir+fullName,"recreate");
   assert(f.IsOpen());
   printf("%d histos are written  to '%s' ...\n",HList.GetEntries(),fullName.Data());
   HList.Write();
   f.Close();
   assert(!f.IsOpen());

}

//.........这里部分代码省略.........
        
        if (isSoftMuon(varValP,nTest))
        {
          if (issig) hists_sig_distr_CutTest[varname[ihist]]->Fill(varValP[ihist],weight); // Variable distribution with the rest of soft MuID cuts applied (all the soft muon cuts if the variable does't belong to the soft cut)
          if (isbkg) hists_bkg_distr_CutTest[varname[ihist]]->Fill(varValP[ihist],weight);
        }
        if ( !isQQVar && isSoftMuon(varValM,nTest) )
        {
          if (issig) hists_sig_distr_CutTest[varname[ihist]]->Fill(varValM[ihist],weight);
          if (isbkg) hists_bkg_distr_CutTest[varname[ihist]]->Fill(varValM[ihist],weight);
        }
        ////_______
        
        
        dummyA = static_cast<TObjArray*>(aMinvCuts->FindObject(Form("mInv_%s",varname[ihist])));
        
        TH1* dummyH(0x0);
        Double_t nBins = hists_sig[varname[ihist]]->GetNbinsX();
        for ( Int_t j=1 ; j<=nBins ; j++)
        {
          //// Cuts on MuID variables
          Double_t cutVal = hists_sig[varname[ihist]]->GetBinLowEdge(j);
          if (Cut(varname[ihist],i,cutVal))
          {
            if (issig) hists_sig[varname[ihist]]->Fill(hists_sig[varname[ihist]]->GetBinCenter(j),weight); // Variable cut distribution
            if (isbkg) hists_bkg[varname[ihist]]->Fill(hists_bkg[varname[ihist]]->GetBinCenter(j),weight);
            
            if ( !(dummyH = static_cast<TH1*>(dummyA->FindObject(Form("mInv_%s_Cut%2.4f",varname[ihist],cutVal)))) ) // Creates the minv histogram and store in array (only once)
            {
              dummyH = new TH1D(Form("mInv_%s_Cut%2.4f",varname[ihist],cutVal),Form("Dimuon m_{inv} distribution (cut on %s: %2.4f)",varname[ihist],cutVal),600,0.,15.0);
              dummyA->Add(dummyH);
            }
            
            dummyH->Fill(mass,weight); // Invariant mass histo for each variable cut
            
            if ( isSoftMuon(varValP,nTest) && isSoftMuon(varValM,nTest) ) // If the variable is not in the soft MuID cut, the full soft cut is applied and the variable cut distribution is filled
            {
              if (issig) hists_sig_SofMu[varname[ihist]]->Fill(hists_sig[varname[ihist]]->GetBinCenter(j),weight);
              if (isbkg) hists_bkg_SofMu[varname[ihist]]->Fill(hists_bkg[varname[ihist]]->GetBinCenter(j),weight);
            }
          }
        }
      }
      //=======
      i++;
    }
  }
  
  //====== Save histos ======
  TFile *f = new TFile(Form("histos_%s_%s_%s_Pt%2.1f_%2.1f_Y%2.1f_%2.1f.root",fileType,particle,fIncludeHighPurity ? "HPincl" : "NoHPincl",Ptmin,Ptmax,Ymin,Ymax),"RECREATE");
  
  for (int i=0; i<nvar; i++)
  {
    aSigC->Add(hists_sig[varname[i]]);
    aBkgC->Add(hists_bkg[varname[i]]);
    
    aSig->Add(hists_sig_distr[varname[i]]);
    aBkg->Add(hists_bkg_distr[varname[i]]);
  
    aDistr->Add(hists_distr[varname[i]]);
    
    aSigSoftMuCutsButOne->Add(hists_sig_distr_CutTest[varname[i]]);
    aBkgSoftMuCutsButOne->Add(hists_bkg_distr_CutTest[varname[i]]);
    
    aSigSoftMuCuts->Add(hists_sig_SofMu[varname[i]]);
    aBkgSoftMuCuts->Add(hists_bkg_SofMu[varname[i]]);
  }
  
  aSigC->Write("SigSingleCut", TObject::kOverwrite | TObject::kSingleKey);
  aBkgC->Write("BkgSingleCut", TObject::kOverwrite | TObject::kSingleKey);
  
  aSig->Write("SigDistr", TObject::kOverwrite | TObject::kSingleKey);
  aBkg->Write("BkgDistr", TObject::kOverwrite | TObject::kSingleKey);
  
  aDistr->Write("Distr", TObject::kOverwrite | TObject::kSingleKey);
  
  aSigSoftMuCutsButOne->Write("SigDistrSoftMuCutButOne", TObject::kOverwrite | TObject::kSingleKey);
  aBkgSoftMuCutsButOne->Write("BkgDistrSoftMuCutButOne", TObject::kOverwrite | TObject::kSingleKey);
  
  aSigSoftMuCuts->Write("SigSoftMuCut", TObject::kOverwrite | TObject::kSingleKey);
  aBkgSoftMuCuts->Write("BkgSoftMuCut", TObject::kOverwrite | TObject::kSingleKey);
  
  
  aGlob->Add(hMinvw);
  aGlob->Add(hRapw);
  aGlob->Add(hPtw);
  aGlob->Add(hpt);
  aGlob->Add(hptw);
  aGlob->Add(hCent);
  aGlob->Add(hCentw);

  aGlob->Write("Global", TObject::kOverwrite | TObject::kSingleKey);
  
  
  aMinvCuts->Write("MinvWCuts", TObject::kOverwrite | TObject::kSingleKey);
  
  f->Close();
  //=========
  
}

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;
}

void saveHistRoot(const char *cfile){

  cout<<"File to be saved = "<<cfile<<endl;
  TFile *outputFile = new TFile(cfile,"recreate");   

  hdN_dPhiTrkNJet->Write(), hdN_dPhiTrkAJet->Write();
  hdN_dEtaTrkNJet->Write(), hdN_dEtaTrkAJet->Write();
  hdN_dRTrkNJet->Write(), hdN_dRTrkAJet->Write();

  hdN_dPhiTrkTrk_dPhiTrkJet->Write();
  hdN_dEtaTrkTrk_dEtaTrkJet->Write();
  hdN_dRTrkTrk_dRTrkJet->Write();
  hdN_dPhiTrkTrk_dRTrkJet->Write();

  hdN_dEtadPhiTrkNJet->Write(), hdN_dEtadPhiTrkAJet->Write();

  hdNdPt_NJet->Write(), hdNdPt_AJet->Write();
  
  // dN/dpT nomalization
  /*
  for(int i=0; i<ptBins.size(); i++){
    float dbin = hdNdPt_NJet_vbin->GetBinWidth(i+1);

    float dn_N = hdNdPt_NJet_vbin->GetBinContent(i+1);
    float dn_A = hdNdPt_AJet_vbin->GetBinContent(i+1);
    
    float edn_N = hdNdPt_NJet_vbin->GetBinError(i+1);
    float edn_A = hdNdPt_AJet_vbin->GetBinError(i+1);

    hdNdPt_NJet_vbin->SetBinContent(i+1,dn_N/dbin);
    hdNdPt_AJet_vbin->SetBinContent(i+1,dn_A/dbin);

    hdNdPt_NJet_vbin->SetBinError(i+1,edn_N/dbin);
    hdNdPt_AJet_vbin->SetBinError(i+1,edn_A/dbin);
  }
  */
  hdNdPt_NJet_vbin->Write(), hdNdPt_AJet_vbin->Write();

  hdNdZ_NJet->Write(), hdNdZ_AJet->Write();
  hdNdNR_NJet->Write(), hdNdNR_AJet->Write();
  hdNdNRW_NJet->Write(), hdNdNRW_AJet->Write();

  hdNdJetEt_NJet->Write(),  hdNdJetEt_AJet->Write();
  hdNdTrkSPt_NJet->Write(), hdNdTrkSPt_AJet->Write();

  hdNPtEtRatio_NJet->Write(), hdNPtEtRatio_AJet->Write();
  hdNPtEtAsymm_NJet->Write(), hdNPtEtAsymm_AJet->Write();

  dJetDist = outputFile->mkdir("jetDist"), dJetDist->cd();
  hdEtadPhiJetArray.Write();
  
  dTrkDist = outputFile->mkdir("trkDist"), dTrkDist->cd();
  hdEtadPhiTrkArray.Write();
  
  dTrkWDist = outputFile->mkdir("trkWDist"), dTrkWDist->cd();
  hdEtadPhiTrkWArray.Write();


  dJetTagDist = outputFile->mkdir("jetTagDist"), dJetTagDist->cd();
  hdEtadPhiJetTagArray.Write();

  dTrkTagDist = outputFile->mkdir("trkTagDist"), dTrkTagDist->cd();
  hdEtadPhiTrkTagArray.Write();

  dTrkWTagDist = outputFile->mkdir("trkWTagDist"), dTrkWTagDist->cd();
  hdEtadPhiTrkWTagArray.Write();


  
}