C++ TClonesArray::GetEntries方法代码示例

void
PrintAlignment()
{
  AliCDBManager* cdb   = AliCDBManager::Instance();
  cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
  AliCDBEntry*   align = cdb->Get("FMD/Align/Data");
  if (!align) {
    Error("PrintAlignment","didn't alignment data from CDB");
    return;
  }
  
  TClonesArray* array = dynamic_cast<TClonesArray*>(align->GetObject());
  if (!array) {
    Warning("PrintAlignement", "Invalid align data from CDB");
    return;
  }
  Int_t nAlign = array->GetEntries();
  for (Int_t i = 0; i < nAlign; i++) {
    AliAlignObjParams* a = static_cast<AliAlignObjParams*>(array->At(i));
    Double_t ang[3];
    Double_t trans[3];
    a->GetAngles(ang);
    a->GetTranslation(trans);
    std::cout << a->GetVolPath() << "\n" 
	      << "  translation: "
	      << "(" << std::setw(12) << trans[0] 
	      << "," << std::setw(12) << trans[1] 
	      << "," << std::setw(12) << trans[2] << ")\n"
	      << "  rotation:    "
	      << "(" << std::setw(12) << ang[0] 
	      << "," << std::setw(12) << ang[1] 
	      << "," << std::setw(12) << ang[2]  << ")" << std::endl;
    // a->Print();
  }
}

void analyzer()
{

  TString processName = "ZJets";

  TFile* f = TFile::Open(Form("hist_%s.root", processName.Data()), "recreate");

  // Create chain of root trees
  TChain chain("DelphesMA5tune");
  //kisti
  //chain.Add("/cms/home/tjkim/fcnc/sample/ZToLL50-0Jet_sm-no_masses/events_*.root");
  //hep
  chain.Add("/Users/tjkim/Work/fcnc/samples/ZToLL50-0Jet_sm-no_masses/events_*.root");  

  // Create object of class ExRootTreeReader
  ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
  Long64_t numberOfEntries = treeReader->GetEntries();
  
  // Get pointers to branches used in this analysis
  TClonesArray *branchMuon = treeReader->UseBranch("DelphesMA5tuneMuon");
  
  // Book histograms
  TH1 *histDiMuonMass = new TH1F("dimuon_mass","Di-Muon Invariant Mass (GeV)",100, 50, 150);
 
  // Loop over all events

  for(Int_t entry = 0; entry < numberOfEntries; ++entry)
  {
    // Load selected branches with data from specified event
    treeReader->ReadEntry(entry);
 
    int nmuon = 0; 
    for( int i = 0; i < branchMuon->GetEntries(); i++)
    {
      Muon *muon = (Muon*) branchMuon->At(i);
      if( muon->PT <= 20 || abs(muon->Eta) >= 2.4 ) continue;
      nmuon = nmuon + 1 ;
    }

    if( nmuon >= 2){

      Muon *mu1 = (Muon*) branchMuon->At(0);
      Muon *mu2 = (Muon*) branchMuon->At(1);

      // Plot di-muon invariant mass 
      histDiMuonMass->Fill(((mu1->P4()) + (mu2->P4())).M());
    }
 
  }

  // Show resulting histograms
  histDiMuonMass->Write();

  f->Close();

}

void digitsTOF(Int_t nevents, Int_t nfiles){


  TH1F *hadc     = new TH1F("hadc","ADC [bin]",200, -100., 10000.);
  TH1F *hadclog     = new TH1F("hadclog","ADC [bin]",200, -1., 7.);
  
    TTree *treeD=0x0;
  
    TClonesArray *digits =0x0;
  
    for (Int_t event=0; event<nevents; event++) {
      cout << "Event " << event << endl;

      treeD = GetTreeD(event, "TOF", nfiles);
      if ( ! treeD ) {
        cerr << "Event directory not found in " << nfiles <<  " files" << endl;
        exit(1);
      }      

      digits = NULL;
      treeD->SetBranchAddress("TOF", &digits);

      for(Int_t iev=0; iev<treeD->GetEntries(); iev++){
	treeD->GetEntry(iev);

	for (Int_t j = 0; j < digits->GetEntries(); j++) {
	  IlcTOFdigit* dig = dynamic_cast<IlcTOFdigit*> (digits->At(j));
	  hadc->Fill(dig->GetAdc());
	  if(dig->GetAdc()>0)hadclog->Fill(TMath::Log10(dig->GetAdc()));
	}

      }
    }

   TFile fc("digits.TOF.root","RECREATE");
   fc.cd();
   
   hadc->Write();
   hadclog->Write();

   fc.Close();

}

void digitsPHOS(Int_t nevents, Int_t nfiles)
{

 TH1F * hadc = new TH1F ("hadc", "hadc", 100, -10, 200);
 TH1F * hadcLog = new TH1F ("hadclog", "hadclog", 100, -2, 4);
 IlcRunLoader* runLoader = IlcRunLoader::Open("gilc.root","Event","READ");
 IlcPHOSLoader * phosLoader = dynamic_cast<IlcPHOSLoader*>(runLoader->GetLoader("PHOSLoader"));
 
 for (Int_t ievent=0; ievent <nevents; ievent++) {
  // for (Int_t ievent = 0; ievent < runLoader->GetNumberOfEvents(); ievent++) {
   runLoader->GetEvent(ievent) ;
   phosLoader->CleanDigits() ; 
   phosLoader->LoadDigits("READ") ;
   TClonesArray * digits    = phosLoader->Digits() ;
   printf("Event %d contains %d digits\n",ievent,digits->GetEntriesFast());
   
   
   for (Int_t j = 0; j < digits->GetEntries(); j++) {
     
     IlcPHOSDigit* dig = dynamic_cast<IlcPHOSDigit*> (digits->At(j));
     //cout << dig->GetEnergy() << endl;
     hadc->Fill(dig->GetEnergy());
     if(dig->GetEnergy()>0)
       hadcLog->Fill(TMath::Log10(dig->GetEnergy()));
     
   }
   
 }
 
 TFile fc("digits.PHOS.root","RECREATE");
 fc.cd();
 hadc->Write();
 hadcLog->Write();
 fc.Close();


}

    void fill(const JMETree& jme) {
        eventinfo.runNum = jme.run;
        eventinfo.evtNum = jme.evt;
        eventinfo.lumiSec = jme.lumi;
        eventinfo.nPU =  jme.npus->size() ? (*jme.npus)[0] : 0;
        eventinfo.nPUmean = jme.tnpus->size() ? (*jme.tnpus)[0] : 0;
        eventinfo.rhoIso = 0;
        eventinfo.rhoJet = jme.rho;
        eventinfo.triggerBits = 0;
        for(auto i = jme.paths->cbegin() ;  i != jme.paths->cend() ; ++i) {
            //cout << "testing " << *i << endl;
            for(auto i2 = triggernames->cbegin(); i2 != triggernames->cend() ; ++i2) {
                if(*i2 > *i) break;
                //cout << "trying " << *i2 << endl;
                if(i->compare(0,min(i->size(),i2->size()),*i2) == 0) {
                    //std::cout << "setting" << *i << ", " << *i2 << " " << i2-triggernames->cbegin() << endl;
                    eventinfo.triggerBits[i2-triggernames->cbegin()] = true;
                }
            }
            //std::cout << *i << endl;
        }
        //eventinfo.triggerBits = jme.prescales->size();
        eventinfo.pfMET =  (*jme.met_p4)[0].Pt();// not stored atm, to be done later
        eventinfo.pfMETphi = (*jme.met_p4)[0].Phi();
        geneventinfo.id_1 = jme.pdf_id->first;
        geneventinfo.id_2 = jme.pdf_id->second;
        geneventinfo.x_1 = jme.pdf_x->first;
        geneventinfo.x_2 = jme.pdf_x->second;
        geneventinfo.weight = jme.weight;
        //geneventinfo.pthat = jme.pthat; // is missing in the header - but anyway not needed (at the moment)


        vertices->Clear();
        jets->Clear();
        addjets->Clear();


        for(unsigned int j = 0 ; j < jme.position->size() ; ++j) {
            assert(vertices->GetEntries() < vertices->GetSize());
            const int index = vertices->GetEntries();
            new((*vertices)[index]) baconhep::TVertex();
            baconhep::TVertex    *pVertex = (baconhep::TVertex*)(*vertices)[index];
            //pVertex->nTracksFit =
            pVertex->ndof = (*jme.ndof)[j];
            pVertex->chi2 = (*jme.normalizedChi2)[j];
            pVertex->x = (*jme.position)[j].X();
            pVertex->y = (*jme.position)[j].Y();
            pVertex->z = (*jme.position)[j].Z();
        }



        for(unsigned int j = 0 ; j < jme.p4->size() ; ++j) {
            assert(jets->GetEntries() < jets->GetSize());
            const int index = jets->GetEntries();
            new((*jets)[index]) baconhep::TJet();
            baconhep::TJet    *pJet = (baconhep::TJet*)(*jets)[index];
            pJet->pt    = (*jme.p4)[j].pt();
            pJet->eta   = (*jme.p4)[j].eta();
            pJet->phi   = (*jme.p4)[j].phi();
            pJet->mass  = (*jme.p4)[j].mass();
            pJet->ptRaw = (*jme.p4)[j].pt() * (*jme.jec_toraw)[j];
            //pJet->csv = (*jme.pfCombinedSecondaryVertexV2BJetTags)[j];  // only stored in AK4PFCHS, to be done later
            pJet->area  = (*jme.jtarea)[j];
            pJet->genpt    = (*jme.gen_p4)[j].pt();
            pJet->geneta   = (*jme.gen_p4)[j].eta();
            pJet->genphi   = (*jme.gen_p4)[j].phi();
            pJet->genm  = (*jme.gen_p4)[j].mass();
            pJet->betaStar = (*jme.betaStar)[j];
            //std::cout << pJet->pt << std::endl;
        }


        tree->Fill();

    }

void run_radius_correction ()
{
    TStopwatch timer;
    timer.Start();

    gStyle->SetPalette(1,0);
    gStyle->SetHistLineWidth(2);

    // ----  Load libraries   -------------------------------------------------
    gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
    basiclibs();
    gROOT->LoadMacro("$VMCWORKDIR/macro/rich/cbmlibs.C");
    cbmlibs();

   // gROOT->LoadMacro("$VMCWORKDIR/macro/rich/setstyle.C");
   // setphdStyle();
    SetStyles();

    char fileMC[200], fileRec[200];

    sprintf(fileMC,"/d/cbm02/slebedev/rich/JUL09/correction/mc.00.root");
    cout<<fileMC<<endl;
    TFile *f1 = new TFile(fileMC,"R");
    TTree* t1 = f1->Get("cbmsim");
    TFolder *fd1 = f1->Get("cbmroot");
    TClonesArray* fMCTracks = (TClonesArray*) fd1->FindObjectAny("MCTrack");
    t1->SetBranchAddress(fMCTracks->GetName(),&fMCTracks);

    sprintf(fileRec, "/d/cbm02/slebedev/rich/JUL09/correction/reco.00.root");
    TFile *f = new TFile(fileRec,"R");
    TTree* t = f->Get("cbmsim");
    TFolder *fd = f->Get("cbmout");
    TClonesArray *fRichRings = (TClonesArray*) fd->FindObjectAny("RichRing");
    t->SetBranchAddress(fRichRings->GetName(),&fRichRings);
    TClonesArray *fRichMatches = (TClonesArray*) fd->FindObjectAny("RichRingMatch");
    t->SetBranchAddress(fRichMatches->GetName(),&fRichMatches);

    //Int_t fNofBinsX = 40;
    //Int_t fNofBinsY = 50;
    Int_t fNofBinsX = 25;
    Int_t fNofBinsY = 25;
    ///A axis
    TH2D* fh_axisAXYCount;
    TH2D* fh_axisAXYW;
    TH2D* fh_axisAXY;
    TH2D* fh_axisASigma;
    TH2D* mapaxisAXY;


    ///B axis
    TH2D* fh_axisBXYCount;
    TH2D* fh_axisBXYW;
    TH2D* fh_axisBXY;
    TH2D* fh_axisBSigma;
    TH2D* mapaxisBXY;

    mapaxisAXY = new TH2D("fh_mapaxisAXY","dA distribution (x,y);X, [cm];Y, [cm]",fNofBinsX,-200,200,fNofBinsY,-250,250);
    mapaxisBXY = new TH2D("fh_mapaxisBXY","dB distribution (x,y);X, [cm];Y, [cm]",fNofBinsX,-200,200,fNofBinsY,-250,250);
    fh_axisAXYCount = new TH2D("fh_axisAXYCount","A Count",fNofBinsX,-200,200,fNofBinsY,-250,250);
    fh_axisAXYW = new TH2D("fh_axisAXYW","",fNofBinsX,-200,200,fNofBinsY,-250,250);
    fh_axisBXYCount = new TH2D("fh_axisBXYCount","B Count",fNofBinsX,-200,200,fNofBinsY,-250,250);
    fh_axisBXYW = new TH2D("fh_axisBXYW","",fNofBinsX,-200,200,fNofBinsY,-250,250);
    fh_axisAXY = new TH2D("fh_axisAXY","A distribution (x,y);X, [cm];Y, [cm]",fNofBinsX,-200,200,fNofBinsY,-250,250);
    fh_axisBXY = new TH2D("fh_axisBXY","B distribution (x,y);X, [cm];Y, [cm]",fNofBinsX,-200,200,fNofBinsY,-250,250);

    Double_t fMinAaxis = 4.5;
    Double_t fMaxAaxis = 7.5;

    ///Set Mean value of A and B axeses, Compact RICH
    //Double_t fMeanAaxis = 5.06;
    //Double_t fMeanBaxis = 4.65;

    ///Set Mean value of A and B axeses, Large RICH
    Double_t fMeanAaxis = 6.17;
    Double_t fMeanBaxis = 5.6;

    Int_t nEvents=t->GetEntries();
    cout<<" nEvents ="<<nEvents<<endl;
    for(Int_t ievent=0;ievent<nEvents; ievent++ ) {
        cout<<"ievent = "<<ievent;
        CbmRichRing *ring=NULL;
        CbmRichRingMatch *match=NULL;
        t->GetEntry(ievent);
        t1->GetEntry(ievent);
        Int_t nofRings = fRichRings->GetEntries();
        cout<<"  nofRings = "<<nofRings;
        cout<<"  nofMatches = "<< fRichMatches->GetEntries() ;
        cout<<"  nofMCTracks = "<<fMCTracks->GetEntries() << endl;

        for(Int_t iRing=0; iRing < nofRings; iRing++){
            ring = (CbmRichRing*)fRichRings->At(iRing);
            if (!ring) continue;
            match = (CbmRichRingMatch*)fRichMatches->At(iRing);
            if (!match) continue;

            Int_t trackId = match->GetMCTrackID();
            if (trackId == -1) continue;
            if (trackId > fMCTracks->GetEntries()) continue;

            CbmMCTrack* mcTrack = (CbmMCTrack*)fMCTracks->At(trackId);
//.........这里部分代码省略.........

//!PG main function
int 
  selector (TChain * tree, histos & plots, int if_signal)
{
 
 plots.v_hardTAGPt = -99;
 plots.v_softTAGPt = -99;
 plots.v_TAGDProdEta = -99;
 plots.v_TAGDeta = -99;
 plots.v_TAGMinv = -99;
 plots.v_LepLep = -99;
 plots.v_hardLEPPt = -99;
 plots.v_softLEPPt = -99;
 plots.v_LEPDPhi = -99;
 plots.v_LEPDEta = -99;
 plots.v_LEPDR = -99;
 plots.v_LEPMinv = -99;
 plots.v_LEPProdCharge = -99;
 plots.v_hardLEPCharge = -99;
 plots.v_softLEPCharge = -99;
 plots.v_MET = -99;
 
 plots.v_ojets = -99 ;
 plots.v_ojetsCJV = -99 ;
 plots.v_ojetsRegionalCJV = -99 ;
 
 plots.v_ojetsZepp_01 = -99 ;
 plots.v_ojetsZepp_02 = -99 ;
 plots.v_ojetsZepp_03 = -99 ;
 plots.v_ojetsZepp_04 = -99 ;
 plots.v_ojetsZepp_05 = -99 ;
 plots.v_ojetsZepp_06 = -99 ;
 plots.v_ojetsZepp_07 = -99 ;
 plots.v_ojetsZepp_08 = -99 ;
 plots.v_ojetsZepp_09 = -99 ;
 plots.v_ojetsZepp_10 = -99 ;
 plots.v_ojetsZepp_11 = -99 ;
 plots.v_ojetsZepp_12 = -99 ;
 plots.v_ojetsZepp_13 = -99 ;
 plots.v_ojetsZepp_14 = -99 ;
 
 plots.v_decay_Channel_e = -99 ;
 plots.v_decay_Channel_mu = -99 ;
 plots.v_decay_Channel_tau = -99 ;
 
 
 TClonesArray * genParticles = new TClonesArray ("TParticle") ;
 tree->SetBranchAddress ("genParticles", &genParticles) ;
 
 
//  TClonesArray * tagJets = new TClonesArray ("TLorentzVector") ; 
//  tree->SetBranchAddress ("tagJets", &tagJets) ;
 TClonesArray * otherJets_temp = new TClonesArray ("TLorentzVector") ;
 tree->SetBranchAddress (g_KindOfJet.c_str(), &otherJets_temp) ;
//  tree->SetBranchAddress ("otherJets", &otherJets_temp) ;
 
 
 TClonesArray * electrons = new TClonesArray ("TLorentzVector") ;
 tree->SetBranchAddress ("electrons", &electrons) ;
 TClonesArray * muons = new TClonesArray ("TLorentzVector") ;
 tree->SetBranchAddress ("muons", &muons) ;
 TClonesArray * MET = new TClonesArray ("TLorentzVector") ;
 tree->SetBranchAddress ("MET", &MET) ;
 TClonesArray * tracks = new TClonesArray ("TLorentzVector") ;
 tree->SetBranchAddress ("tracks", &tracks) ;
 
 TClonesArray * tagJets = new TClonesArray ("TLorentzVector") ; 
 TClonesArray * otherJets = new TClonesArray ("TLorentzVector") ;
 
  
 int EleId[100];
 float IsolEleSumPt_VBF[100];
 int nEle;
 int EleCharge[30];
 tree->SetBranchAddress ("nEle", &nEle) ;
 tree->SetBranchAddress ("EleId",EleId ) ;
 tree->SetBranchAddress ("IsolEleSumPt_VBF",IsolEleSumPt_VBF ) ;
 tree->SetBranchAddress ("EleCharge",EleCharge ) ;
 
 float IsolMuTr[100];
 int nMu ;
 int MuCharge[30];
 tree->SetBranchAddress ("nMu", &nMu) ;
 tree->SetBranchAddress ("IsolMuTr",IsolMuTr ) ;
 tree->SetBranchAddress ("MuCharge", MuCharge) ;
 

 int IdEvent;
 tree->SetBranchAddress ("IdEvent", &IdEvent) ;
 
 
 int nentries = (int) tree->GetEntries () ;

 
 plots.passedJetAndLepNumberSelections = 0;
 plots.analyzed = 0;
 
 plots.analyzed_ee = 0;
 plots.analyzed_mumu = 0;
 plots.analyzed_tautau = 0;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
			strcpy(head_folder,"_Tops_Events_WI_Matching/");
		else
			strcpy(head_folder,"_Tops_Events_WI/");
		head_folder[0] = channel;
		head_folder[13] = ISR_or_NOT[0];
		head_folder[14] = ISR_or_NOT[1];

		Char_t current_folder[] = "_Tops_MG_1K_AG_WI_003/";
		current_folder[0] = channel;
		current_folder[15] = ISR_or_NOT[0];
		current_folder[16] = ISR_or_NOT[1];

		Char_t unidad = 0x30 + iRun%10;
		Char_t decena = 0x30 + int(iRun/10)%10;
		Char_t centena = 0x30 + int(iRun/100)%10;

		current_folder[18] = centena;
		current_folder[19] = decena;
		current_folder[20] = unidad;

		Char_t *file_delphes;
		file_delphes = (Char_t*) malloc(512*sizeof(Char_t));
		strcpy(file_delphes,local_path);
		strcat(file_delphes,head_folder);
		strcat(file_delphes,current_folder);
		strcat(file_delphes,"Events/run_01/output_delphes.root");

		cout << "\nReading the file: \nDelphes: " << file_delphes << endl;

		chain_Delphes.Add(file_delphes);
		// Objects of class ExRootTreeReader for reading the information
		ExRootTreeReader *treeReader_Delphes = new ExRootTreeReader(&chain_Delphes);

		Long64_t numberOfEntries = treeReader_Delphes->GetEntries();

		// Get pointers to branches used in this analysis
		TClonesArray *branchJet = treeReader_Delphes->UseBranch("Jet");
		TClonesArray *branchMissingET = treeReader_Delphes->UseBranch("MissingET");

		cout << endl;
		cout << " Number of Entries Delphes = " << numberOfEntries << endl;
		cout << endl;

		// particles, jets and vectors
		MissingET *METpointer;
		TLorentzVector *vect_currentJet = new TLorentzVector;
		TLorentzVector *vect_auxJet = new TLorentzVector;
		TLorentzVector *vect_leading = new TLorentzVector;
		Jet *currentJet = new Jet;
		Jet *auxJet = new Jet;
		TRefArray array_temp;

		// Temporary variables
		Double_t MET = 0.0; // Missing transverse energy
		Double_t delta_phi = 0.0; // difference between the phi angle of MET and the jet
		Double_t transverse_mass = 0.0; // Transverse mass
		Double_t HT = 0.0; // Sum of jets' PT
		Double_t HT_R1 = 0.0; // Sum of jets' PT which are in the same hemisphere of the ISR jet hemisphere
		Double_t HT_R2 = 0.0; // Sum of jets' PT which are in the opposite hemisphere of the ISR jet hemisphere
		Double_t ISR_Eta = 0.0; // Pseudorapidity of the ISR jet
		Int_t number_Btags = 0; // Number of B jets per event
		Int_t ISR_Btags = 0; // Number of BTags which are also ISR jets
		Double_t delta_PT_jet = 0.0; // |PT-<PT>|
		Double_t PT_sum = 0.0; // sum(PT)
		Double_t PT_aver = 0.0; // <PT>
		Double_t Delta_eta_aver = 0.0; // sum_i|eta-eta_i|/(Nj-1)

void GeneralExample(const char *inputFile, const char *outputFile)
{
  //  gSystem->Load("libDelphes");

  // Create chain of root trees
  TChain chain("Delphes");
  chain.Add(inputFile);
  
  // Create object of class ExRootTreeReader
  ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
  Long64_t numberOfEntries = treeReader->GetEntries();
  
  // Get pointers to branches used in this analysis
  TClonesArray *branchGenJet = treeReader->UseBranch("GenJet");
  TClonesArray *branchJet = treeReader->UseBranch("Jet");
  TClonesArray *branchRho = treeReader->UseBranch("Rho");
  TClonesArray *branchGlobalRho = treeReader->UseBranch("GlobalRho");
  TClonesArray *branchNPU = treeReader->UseBranch("NPU");
  
  TClonesArray *branchElectron = treeReader->UseBranch("Electron");
  TClonesArray *branchMuon = treeReader->UseBranch("Muon");
  TClonesArray *branchPhoton = treeReader->UseBranch("Photon");

  // Constituents will be 0 otherwise
  TClonesArray *branchEFlowTrack = treeReader->UseBranch("EFlowTrack");
  TClonesArray *branchEFlowTower = treeReader->UseBranch("EFlowTower");
  TClonesArray *branchEFlowMuon = treeReader->UseBranch("EFlowMuon");
  TClonesArray *branchGenParticle = treeReader->UseBranch("Particle");
  TClonesArray *branchBeamSpotParticle = treeReader->UseBranch("BeamSpotParticle");

  TClonesArray *branchMissingET = treeReader->UseBranch("MissingET");
  TClonesArray *branchGenMissingET = treeReader->UseBranch("GenMissingET");
  TClonesArray *branchPileUpJetIDMissingET = treeReader->UseBranch("PileUpJetIDMissingET");
  TClonesArray *branchPuppiMissingET = treeReader->UseBranch("PuppiMissingET");


  bool verbose = true;
  bool listJetTowers = false;
  bool listMET = true;
  bool listRho = false;

  bool treeMET = true;

  TFile *f;
  TTree *t;
    //  TH1F *hdmet, hdpujidmet, hdpuppimet;
    //  TH1F *hfdmet, hfdpujidmet, hfdpuppimet;
    //  TH2F *h2met, h2pujidmet, h2puppimet;
  float genmet, met, puppimet, pujidmet;
  float genmetx, metx, puppimetx, pujidmetx;
  float genmety, mety, puppimety, pujidmety;


  if (treeMET) {
    f = new TFile(outputFile,"RECREATE");
    t = new TTree("t","t");
    t->Branch("genmet",&genmet,"genmet/F");
    t->Branch("met",&met,"met/F");
    t->Branch("pujidmet",&pujidmet,"pujidmet/F");
    t->Branch("puppimet",&puppimet,"puppimet/F");
    t->Branch("genmetx",&genmetx,"genmetx/F");
    t->Branch("metx",&metx,"metx/F");
    t->Branch("pujidmetx",&pujidmetx,"pujidmetx/F");
    t->Branch("puppimetx",&puppimetx,"puppimetx/F");
    t->Branch("genmety",&genmety,"genmety/F");
    t->Branch("mety",&mety,"mety/F");
    t->Branch("pujidmety",&pujidmety,"pujidmety/F");
    t->Branch("puppimety",&puppimety,"puppimety/F");

  }

  // Loop over all events
  for(Int_t entry = 0; entry < numberOfEntries; ++entry)
  {
    // Load selected branches with data from specified event
    treeReader->ReadEntry(entry);
  
    if (listMET||verbose||entry%5000==0) cout << "Event " << entry << " / " << numberOfEntries << endl;
    
    for (int i = 0 ;  i < branchGenMissingET->GetEntries() ; i++) {
      MissingET *m = (MissingET*) branchGenMissingET->At(i);
      if (verbose || listMET) cout << "Gen MissingET: " << m->MET << endl;
      genmet = m->MET;
      genmetx = m->MET*cos(m->Phi);
      genmety = m->MET*sin(m->Phi);
    }

    for (int i = 0 ;  i < branchMissingET->GetEntries() ; i++) {
      MissingET *m = (MissingET*) branchMissingET->At(i);
      if (verbose || listMET) cout << "MissingET: " << m->MET << endl;
      met = m->MET;
      metx = m->MET*cos(m->Phi);
      mety = m->MET*sin(m->Phi);

    }

    for (int i = 0 ;  i < branchPileUpJetIDMissingET->GetEntries() ; i++) {
      MissingET *m = (MissingET*) branchPileUpJetIDMissingET->At(i);
      if (verbose || listMET) cout << "MissingET using PileUpJetID: " << m->MET << endl;
      pujidmet = m->MET;
//.........这里部分代码省略.........

void RecoTrack_wSks(char *ridffile="./ridf/sm_com/sdaq02/run0208.ridf", Int_t nevent=100000000000){
//void RecoTrack_wSks(char *ridffile="./ridf/sm_com/sdaq02/run0208.ridf", Int_t nevent=100){

  TArtStoreManager * sman = TArtStoreManager::Instance();
  TArtEventStore *estore = new TArtEventStore;
  estore->Open(ridffile);

  TArtBigRIPSParameters * bripsparameters = TArtBigRIPSParameters::Instance();
  bripsparameters->LoadParameter((char*)"db/BigRIPSPlastic.xml");

  TArtSAMURAIParameters * samuraiparameters = TArtSAMURAIParameters::Instance();
  samuraiparameters->LoadParameter((char*)"db/SAMURAIFDC1.xml");
  samuraiparameters->LoadParameter((char*)"db/SAMURAIFDC2.xml");

  TArtCalibPlastic * calpla     = new TArtCalibPlastic;
  TArtCalibSAMURAI * calsamurai = new TArtCalibSAMURAI;

  TArtRecoFragment *reco = new TArtRecoFragment();
  TClonesArray *fdc1trks = (TClonesArray*)sman->FindDataContainer("SAMURAIFDC1Track");
  TClonesArray *fdc2trks = (TClonesArray*)sman->FindDataContainer("SAMURAIFDC2Track");
  TClonesArray *frags = (TClonesArray*)sman->FindDataContainer("SAMURAIFragment");

  calsamurai->LoadDCTDCDistribution((char*)"db/dc/run0203.root"); // for 3T Brho scan analysis

  Double_t fdc1chi2, fdc2chi2;
  Double_t delta, brho, trchi2, trpos[3], trl2hod;
  Double_t f3plat, f3plaq;
  Double_t f13plat, f13plaq;
  Int_t trstatus, trnhit;

  TFile *fout = new TFile("fout.root","RECREATE");
  TTree *otree = new TTree("otree","otree");
  otree->Branch("f3plaq",&f3plaq,"f3plaq/D");
  otree->Branch("f3plat",&f3plat,"f3plat/D");
  otree->Branch("f13plaq",&f13plaq,"f13plaq/D");
  otree->Branch("f13plat",&f13plat,"f13plat/D");
  otree->Branch("tr1chi2",&fdc1chi2,"tr1chi1/D");
  otree->Branch("tr2chi2",&fdc2chi2,"tr2chi1/D");
  otree->Branch("delta",&delta,"delta/D");

  otree->Branch("trstatus",&trstatus,"trstatus/I");
  otree->Branch("brho",&brho,"brho/D");
  otree->Branch("trchi2",&trchi2,"trchi2/D");
  otree->Branch("trnhit",&trnhit,"trnhit/I");
  otree->Branch("trpos",trpos,"trpos[3]/D");
  otree->Branch("trl2hod",&trl2hod,"trl2hod/D");

  int neve = 0;
  while(estore->GetNextEvent()&&neve < nevent){

    if(neve%1000==0) cout << "event: " << neve << endl;

    //    if(neve<33410) {neve ++; continue;} if(neve%100==0) cout << "event: " << neve << endl;

    calpla->ClearData();
    calsamurai->ClearData();
    reco->ClearData();

    calpla->ReconstructData();
    calsamurai->ReconstructData();
    reco->ReconstructData();

    f3plat = -9999; f3plaq = -9999;
    f13plat = -9999; f13plaq = -9999;
    fdc1chi2 = -9999; fdc2chi2 = -9999;
    delta = -9999; 
    brho = -9999;
    trchi2 = -9999;
    for(int i=0;i<3;i++) trpos[i] = -9999;
    trl2hod = -9999;
    trstatus = 0;
    trnhit = -1;

    TArtPlastic *pla = calpla->FindPlastic((char*)"F3pl");
    if(pla){
      if(pla->GetTLRaw() > 0 && pla->GetTRRaw() > 0)
	f3plat = ((Double_t)pla->GetTLRaw() + (Double_t)pla->GetTRRaw())/2;
      if(pla->GetQLRaw() > 0 && pla->GetQRRaw() > 0)
	f3plaq = pla->GetQAveRaw();
    }
    else{
      cout << "cannot find f3pla object!!" << endl;
    }

    pla = calpla->FindPlastic((char*)"F13pl-1");
    if(pla){
      if(pla->GetTLRaw() > 0 && pla->GetTRRaw() > 0)
	f13plat = ((Double_t)pla->GetTLRaw() + (Double_t)pla->GetTRRaw())/2;
      if(pla->GetQLRaw() > 0 && pla->GetQRRaw() > 0)
	f13plaq = pla->GetQAveRaw();
    }
    else{
      cout << "cannot find f13pla-1 object!!" << endl;
    }

    int nfdc1trks = fdc1trks->GetEntries();
    int nfdc2trks = fdc2trks->GetEntries();
    int nfrags = frags->GetEntries();

    cout << "num of fdc1trks in event: " << nfdc1trks << endl;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  // Variables to store acceptances and uncertainties (per input file)
  vector<Double_t> nEvtsv, nSelv;
  vector<Double_t> nSelCorrv, nSelCorrVarv;
  vector<Double_t> accv, accCorrv;
  vector<Double_t> accErrv, accErrCorrv;

  const baconhep::TTrigger triggerMenu("../../BaconAna/DataFormats/data/HLT_50nsGRun");

  //
  // loop through files
  //
  for(UInt_t ifile=0; ifile<fnamev.size(); ifile++) {  

    // Read input file and get the TTrees
    cout << "Processing " << fnamev[ifile] << " ..." << endl;
    infile = TFile::Open(fnamev[ifile]); 
    assert(infile);

    eventTree = (TTree*)infile->Get("Events"); assert(eventTree);  
    eventTree->SetBranchAddress("Info",              &info); TBranch *infoBr     = eventTree->GetBranch("Info");
    eventTree->SetBranchAddress("GenEvtInfo",         &gen); TBranch *genBr      = eventTree->GetBranch("GenEvtInfo");
    eventTree->SetBranchAddress("GenParticle", &genPartArr); TBranch *genPartBr  = eventTree->GetBranch("GenParticle");
    eventTree->SetBranchAddress("Electron",   &electronArr); TBranch *electronBr = eventTree->GetBranch("Electron");
    eventTree->SetBranchAddress("PV",   &vertexArr); TBranch *vertexBr = eventTree->GetBranch("PV");

    nEvtsv.push_back(0);
    nSelv.push_back(0);
    nSelCorrv.push_back(0);
    nSelCorrVarv.push_back(0);

    //
    // loop over events
    //
    for(UInt_t ientry=0; ientry<eventTree->GetEntries(); ientry++) {
      genBr->GetEntry(ientry);
      genPartArr->Clear(); genPartBr->GetEntry(ientry);
      infoBr->GetEntry(ientry);

      Int_t glepq1=-99;
      Int_t glepq2=-99;

      if (fabs(toolbox::flavor(genPartArr, BOSON_ID))!=LEPTON_ID) continue;
      TLorentzVector *vec=new TLorentzVector(0,0,0,0);
      TLorentzVector *lep1=new TLorentzVector(0,0,0,0);
      TLorentzVector *lep2=new TLorentzVector(0,0,0,0);
      toolbox::fillGen(genPartArr, BOSON_ID, vec, lep1, lep2,&glepq1,&glepq2,1);
      if(vec->M()<MASS_LOW || vec->M()>MASS_HIGH) continue;      
      delete vec; delete lep1; delete lep2;

      vertexArr->Clear();
      vertexBr->GetEntry(ientry);
      double npv  = vertexArr->GetEntries();
      Double_t weight=gen->weight;
      if(doPU>0) weight*=h_rw->GetBinContent(h_rw->FindBin(info->nPUmean));

      nEvtsv[ifile]+=weight;        

      // trigger requirement               
      if (!isEleTrigger(triggerMenu, info->triggerBits, kFALSE)) continue;
      
      // good vertex requirement
      if(!(info->hasGoodPV)) continue;
      
      electronArr->Clear();
      electronBr->GetEntry(ientry);

//.........这里部分代码省略.........
      ring_histo[i] = new TH1D(name2,"",60,0.,3.0);
      ring_histo[i]->SetXTitle("E/p");
      char namerpve[100];
      sprintf(namerpve,"ring_pve_%i",i);
      ring_pve[i] = new TH2F(namerpve,"",20,0,20.0,20,0,20.0);
      ring_pve[i]->SetXTitle("E (GeV)");
      ring_pve[i]->SetYTitle("p (GeV)");
    }

  for(int i = 0; i < 6; i++){
    char jname[100];
    sprintf(jname,"jan_pve_%i",i);
    jan_pve[i] = new TH2F(jname,"",120,0,3.0,20,0,20.0);
    jan_pve[i]->SetXTitle("E/p");
    jan_pve[i]->SetYTitle("p (GeV)");
  }
  //global graphics functions
  gStyle->SetOptStat("oue");
  gStyle->SetOptFit(111);
  gStyle->SetCanvasColor(10);
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetOptTitle(0);
  gStyle->SetPalette(1);
  gStyle->SetStatColor(0);

  StEmcOfflineCalibrationTrack* track = new StEmcOfflineCalibrationTrack();
  StEmcOfflineCalibrationCluster* cluster = new StEmcOfflineCalibrationCluster();
  tree->SetBranchAddress("clusters",&cluster);

  //**********************************************//
  //Loop Over Tracks, Fill Histograms             //
  //**********************************************//

  int nentries = tree->GetEntries();
  cout<<nentries<<endl;
  int ngoodhit = 0;
  int nplt10 = 0;
  int nnosis = 0;
  int nfinal = 0;
  int nbsmdgood = 0;
  int nnottrig = 0;
  int nfidu = 0;
  int nenterexit = 0;
  for(int j=0; j<nentries; j++){
    tree->GetEntry(j);
    track = &(cluster->centralTrack);
    TClonesArray *tracks = cluster->tracks;
    if(j % 500000 == 0) cout<<"reading "<<j<<" of "<<nentries<<endl;

    httrig->Fill((float)track->htTrig);

    if(track->charge > 0)pplus->Fill(track->p);
    if(track->charge < 0)pminus->Fill(track->p);

    int bsmdeadctot = 0;
    int bsmdpadctot = 0;
    for(int i = 0; i < 11; i++){
      if(track->smde_adc[i] > track->smde_pedestal[i])bsmdeadctot += track->smde_adc[i] - track->smde_pedestal[i];
      if(track->smdp_adc[i] > track->smdp_pedestal[i])bsmdpadctot += track->smdp_adc[i] - track->smdp_pedestal[i];
    }

    double dR = TMath::Sqrt(track->deta*track->deta + track->dphi*track->dphi);


    double scaled_adc = (track->tower_adc[0] - track->tower_pedestal[0]) / track->p;

bdump(){
  /* 
  * Loading the evio-root shared library. Change the path if running the script
  * from a different directory.
  */
  gSystem->Load("../lib/libEvioRoot.so");  
  TEvioFileReader *reader = new TEvioFileReader();
  TEvioDataEvent  *event  = new TEvioDataEvent();
  reader->open("/home/celentan/WORK/BeamDump/GenoaTest/data/run_JLabFadc_000626.evio.0");
  
  TH1D *hADC0 = new TH1D("hADC0","",500,0.0,500.); 
  TH1D *hADC00 = new TH1D("hADC00","",500,0.0,500.);
  
  TH1D *hADC1 = new TH1D("hADC1","",500,0.0,500.);
  TH1D *hADC2 = new TH1D("hADC2","",500,0.0,500.);
  TH1D *hADC3 = new TH1D("hADC3","",500,0.0,500.);
  TH1D *hADC4 = new TH1D("hADC4","",500,0.0,500.);
  
  
  
  TH1D *hCHARGE0 = new TH1D("hCHARGE0","",2500,-1000.0,20000.);
  TH1D *hCHARGE1 = new TH1D("hCHARGE1","",500,-1000.0,20000.);
  TH1D *hCHARGE2 = new TH1D("hCHARGE2","",500,-1000.0,20000.);
  TH1D *hCHARGE3 = new TH1D("hCHARGE3","",500,-1000.0,20000.);
  TH1D *hCHARGE4 = new TH1D("hCHARGE4","",500,-1000.0,20000.);
  
  
  TH2D *hSTRIP0,*hSTRIP1,*hSTRIP2,*hSTRIP3,*hSTRIP4;
  
  
  
  
  int icounter = 0;
  int Nevents = reader->getEntries();
  reader->next();reader->next();
  
  
  hSTRIP0=new TH2D("hSTRIP0","",Nevents/600,0,Nevents,500,-1000,15000.);
  hSTRIP1=new TH2D("hSTRIP1","",Nevents/600,0,Nevents,500,-1000,20000.);
  hSTRIP2=new TH2D("hSTRIP2","",Nevents/600,0,Nevents,500,-1000,20000.);
  hSTRIP3=new TH2D("hSTRIP3","",Nevents/600,0,Nevents,500,-1000,20000.);
  hSTRIP4=new TH2D("hSTRIP4","",Nevents/600,0,Nevents,500,-1000,20000.);
  
  
  hSTRIP00=new TH2D("hSTRIP00","",Nevents/600,0,Nevents,500,-1000,15000.);
  hSTRIP10=new TH2D("hSTRIP10","",Nevents/600,0,Nevents,500,-1000,20000.);
  hSTRIP20=new TH2D("hSTRIP20","",Nevents/600,0,Nevents,500,-1000,20000.);
  hSTRIP30=new TH2D("hSTRIP30","",Nevents/600,0,Nevents,500,-1000,20000.);
  hSTRIP40=new TH2D("hSTRIP40","",Nevents/600,0,Nevents,500,-1000,20000.);
  
  TClonesArray *flash;
  TClonesArray *tdcL;
  TCanvas *c=new TCanvas("c","c");
  while( reader->next() == true ){
    
    flash=0;
    tdcL=0;
    hADC00->Reset();
    if (icounter==(Nevents-10)) break;
    
  
    // if (icounter==10000) break;
    /* reader->getEvent(*event);*/
    tdcL = reader->getTDC(1,57607);//tag mother, tag TDC 
    if (tdcL!=0){
      int nentries = tdcL->GetEntries();   
      for(int loop = 0; loop < nentries; loop++){	
	TTDCClass *tdc = (TTDCClass *) tdcL->At(loop);
	int slot = tdc->GetSlot();
	int channel = tdc->GetChannel();
	int value = tdc->GetValue();
	
	//cout<<slot<<" "<<channel<<" "<<value<<endl;
      }
    }
    
    
    flash= reader->getFlashADC(1,57601); //tag mother, tag FADC
    
    
    if (flash!=0){
      int nentries1 = flash->GetEntries();
      
      vector< int > v0,v1,v2,v3,v4;
      double q0,q1,q2,q3,q4;
      double ped0,ped1,ped2,ped3,ped4; 
      
      
      for(int loop = 0; loop < nentries1; loop++){
	TADCClass *adc = (TADCClass *) flash->At(loop);

	int slot    = adc->GetSlot();
	int channel = adc->GetChannel();
	int nsamples = adc->GetRows();
	int value;
		
	//	cout << " SLOT = " << adc->GetSlot() << "  CHANNEL = " << adc->GetChannel() << " "<<nsamples<<endl;
	for(int i = 0; i < nsamples; i++){
	  value = adc->GetValue(i);
	  if ((slot == 4) && (channel == 1)){
//.........这里部分代码省略.........

int
main (int argc, char *argv[])
{
  if (argc != 3)
    {
      cout << "Usage: " << argv[0] << " INPUT_FILE OUTPUT_FILE" << endl;
      cout << "  or:  " << argv[0] << " INPUT_LIST OUTPUT_FILE" << endl;
      cout << "Run vertexing performance analysis on INPUT_FILE and store output to OUTPUT_FILE. If the" << endl;
      cout << "first argument is a text file, analysis is run on each file listed inside." << endl;
      cout << endl;

      return 0;
    }
  string inputFile = argv[1], outputFile = argv[2], upperInputFile;
  upperInputFile.resize (inputFile.length ());
  transform (inputFile.begin (), inputFile.end (), upperInputFile.begin (), ::toupper);

  // so canvases don't appear on the screen when being created
  // very useful when running on the OSU T3 from CERN
  gROOT->SetBatch();

  //Create chain of root trees
  TChain chain("Delphes");

  if (upperInputFile.length () < 5 || upperInputFile.substr (upperInputFile.length () - 5, 5) != ".ROOT")
  {
    ifstream fin (inputFile);
    string line;
    while(getline(fin, line))
      {
	chain.Add(line.c_str());
      }
    fin.close();
  }
  else
    chain.Add(inputFile.c_str());

  
  // Create object of class ExRootTreeReader
  ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
  Long64_t numberOfEntries = treeReader->GetEntries();
  cout << "processing " << numberOfEntries << " events..." << endl << endl;

  // Get pointers to branches used in this analysis
  TClonesArray *branchTrack = treeReader->UseBranch("Track");
  TClonesArray *branchNPU = treeReader->UseBranch("NPU");
  TClonesArray *branchCluster = treeReader->UseBranch("Cluster");
  TClonesArray *branchVBFquarks = treeReader->UseBranch("VBFquarks");


  ///////////////////////////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////////////////////////
  // Book histograms
  ///////////////////////////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////////////////////////

  TH1::SetDefaultSumw2();

  // VBF Quark Histograms
  TH1D *hVBFQuarkEta   = new TH1D("VBFQuarkEta", ";VBF quark #eta", 1000, -5.0, 5.0);
  TH1D *hVBFQuarkPt   = new TH1D("VBFQuarkPt", ";VBF quark p_{T} [GeV]", 1000, 0.0, 1000.0);
  TH2D *hVBFQuarkEtaVsPhi   = new TH2D("VBFQuarkPt", ";VBF quark #phi;VBF quark #eta", 500, -3.5, 3.5, 500, -5.0, 5.0);

  // Track Histograms
  TH1D *hTrackPt   = new TH1D("TrackPt", "; Track p_{T} [GeV]", 400, 0.0, 200.0);
  TH1D *hVBFTrackPt   = new TH1D("VBFTrackPt", ";VBF track p_{T} [GeV]", 400, 0.0, 200.0);
  TH1D *hPUTrackPt   = new TH1D("PUTrackPt", ";PU track p_{T} [GeV]", 400, 0.0, 200.0);
  TH1D *hTrackX = new TH1D("TrackX", "; Track x [mm]", 300, -3, 3);
  TH1D *hTrackY = new TH1D("TrackY", "; Track y [mm]", 300, -3, 3);
  TH1D *hTrackZ = new TH1D("TrackZ", "; Track z [mm]", 300, -300, 300);

  TH1D *hTrackEX = new TH1D("TrackEX", "; Track #sigma_{x} [mm]", 100, 0, 1);
  TH1D *hTrackEY = new TH1D("TrackEY", "; Track #sigma_{y} [mm]", 100, 0, 1);
  TH1D *hTrackEZ = new TH1D("TrackEZ", "; Track #sigma_{z} [mm]", 100, 0, 20);

  TH1D *hTrackXSmearingShift = new TH1D("TrackXSmearingShift", "; Track x_{smeared} - x_{true} [mm]", 100, -2, 2);
  TH1D *hTrackYSmearingShift = new TH1D("TrackYSmearingShift", "; Track y_{smeared} - y_{true} [mm]", 100, -2, 2);
  TH1D *hTrackZSmearingShift = new TH1D("TrackZSmearingShift", "; Track z_{smeared} - z_{true} [mm]", 100, -50, 50);

  TH1D *hTrackDisplacementClusterX = new TH1D("TrackDisplacementClusterX", "; |Track x - Cluster x| [mm]", 100, 0, 10);
  TH1D *hTrackDisplacementClusterY = new TH1D("TrackDisplacementClusterY", "; |Track y - Cluster y| [mm]", 100, 0, 10);
  TH1D *hTrackDisplacementClusterZ = new TH1D("TrackDisplacementClusterZ", "; |Track z - Cluster z| [mm]", 100, 0, 300);

  TH1D *hTrackDisplacementInteractionX = new TH1D("TrackDisplacementInteractionX", "; |Track x - Interaction x| [mm]", 100, 0, 10);
  TH1D *hTrackDisplacementInteractionY = new TH1D("TrackDisplacementInteractionY", "; |Track y - Interaction y| [mm]", 100, 0, 10);
  TH1D *hTrackDisplacementInteractionZ = new TH1D("TrackDisplacementInteractionZ", "; |Track z - Interaction z| [mm]", 100, 0, 50);

  // Cluster Histograms
  TH1D *hClusterX = new TH1D("ClusterX", "; Cluster x [mm]", 300, -3, 3);
  TH1D *hClusterEX = new TH1D("ClusterEX", "; Cluster #sigmax [mm]", 100, 0, 0.1);
  TH1D *hClusterY = new TH1D("ClusterY", "; Cluster y [mm]", 300, -3, 3);
  TH1D *hClusterEY = new TH1D("ClusterEY", "; Cluster #sigmay [mm]", 100, 0, 0.1);
  TH1D *hClusterZ = new TH1D("ClusterZ", "; Cluster z [mm]", 300, -300, 300);
  TH1D *hClusterEZ = new TH1D("ClusterEZ", "; Cluster #sigmaz [mm]", 100, 0, 0.1);
  TH2D *hClusterTransversePosition = new TH2D("ClusterTransversePosition", "; Cluster x [mm]; Cluster y [mm]", 300, -3,3, 300, -3, 3);
  TH2D *hClusterRZPosition = new TH2D("ClusterRZPosition", "; Cluster z [mm]; Cluster r [mm]", 100, -300, 300, 100, 0, 1.5);

  // Reco PV Histograms
  TH1D *hRecoPVX = new TH1D("RecoPVX", "; Reco PV x [mm]", 300, -3, 3);
  TH1D *hRecoPVY = new TH1D("RecoPVY", "; Reco PV y [mm]", 300, -3, 3);
//.........这里部分代码省略.........

void TestEMCALSDigit()
{
  
  // Getting EMCAL Detector and Geometry.
  
  AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(),"read");
  
  if (rl == 0x0)
    cout<<"Can not instatiate the Run Loader"<<endl;
  
  rl->LoadgAlice();//Needed to get geometry
  
  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>
    (rl->GetDetectorLoader("EMCAL"));
  
  TGeoManager::Import("geometry.root");
  
  AliRun * alirun   = rl->GetAliRun(); // Needed to get Geometry
  AliEMCALGeometry * geom ;
  if(alirun){
    AliEMCAL * emcal  = (AliEMCAL*)alirun->GetDetector("EMCAL");
    geom = emcal->GetGeometry();
  }
  
  if (geom == 0) cout<<"Did not get geometry from EMCALLoader"<<endl;
  else   geom->PrintGeometry();
  
  //Load Digits
  rl->LoadSDigits("EMCAL");
  
  //Get maximum number of events
  Int_t maxevent =  rl->GetNumberOfEvents();
  cout<<"Number of events "<<maxevent<<endl;
  //maxevent = 10 ;
  
  Int_t iEvent  = -1 ;
  Float_t amp   = -1 ;
  Float_t time  = -1 ;
  Int_t id      = -1 ;
  Int_t iSupMod =  0 ;
  Int_t iTower  =  0 ;
  Int_t iIphi   =  0 ;
  Int_t iIeta   =  0 ;
  Int_t iphi    =  0 ;
  Int_t ieta    =  0 ;
  
  AliEMCALDigit * dig;
  
  for ( iEvent=0; iEvent<maxevent; iEvent++)
    {
      cout <<  " ======> Event " << iEvent << endl ;
      //Load Event
      rl->GetEvent(iEvent);
      
      //Fill array of digits
      TClonesArray *digits = emcalLoader->SDigits();
      
      //Get digits from the list      
      for(Int_t idig = 0; idig< digits->GetEntries();idig++){
	//cout<<">> idig "<<idig<<endl;
	dig = static_cast<AliEMCALDigit *>(digits->At(idig)) ;
	
	if(dig != 0){
	  id   = dig->GetId() ; //cell (digit) label
	  amp  = dig->GetAmplitude(); //amplitude in cell (digit)
	  time = dig->GetTime();//time of creation of digit after collision
	  
	  cout<<"Cell ID "<<id<<" Amp "<<amp<<endl;//" time "<<time<<endl;
	  
	  //Geometry methods  
	  if(geom){
	    geom->GetCellIndex(id,iSupMod,iTower,iIphi,iIeta); 
	    //Gives SuperModule and Tower numbers
	    geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
					      iIphi, iIeta,iphi,ieta);
	    //Gives label of cell in eta-phi position per each supermodule
	    cout<< "SModule "<<iSupMod<<"; Tower "<<iTower <<"; Eta "<<iIeta
		<<"; Phi "<<iIphi<<"; Cell Eta "<<ieta<<"; Cell Phi "<<iphi<<endl;
	  }
	}
	else
	  cout<<"Digit pointer 0x0"<<endl;
      }
      
    }


}