///////////////////////////////////////////////////////////////////////////////
//
// 
// File: D0ChunkAnalyze.cpp
//
// Purpose:  skeleton for performing data analysis at the D0 chunk level
//
// Created:    24-JUL-2002  Marco Verzocchi
//
// History:    05-SEP-2003 Andre Sznajder ( Adapt for Diffractve Analysis )
//             
//
///////////////////////////////////////////////////////////////////////////////

// Dependencies (#includes)
#include "analysis_example/D0ChunkAnalyze.hpp"

#include "rcp/RCP.hpp"

#include "framework/Registry.hpp"

#include "edm/IDSelector.hpp"
#include "edm/LinkPtr.hpp" 

#include "prod_history/HistorySelector.hpp"
#include "prod_history/HistorySelectorLatest.hpp"
#include "prod_history/HistoryChunk.hpp"


#include "em_evt/EMparticleChunkSelector.hpp"
#include "em_evt/EMparticle.hpp"
#include "em_evt/EMCluster.hpp"
#include "em_evt/EMQualityInfo.hpp"

#include "chpart_evt/ChargedParticleChunk.hpp"
#include "chpart_evt/ChargedParticle.hpp"
#include "vertexutil/TrackDCA.hpp"
#include "d0track/D0HitMask.hpp"

#include "muo_evt/MuonParticleChunk.hpp"
#include "muo_evt/MuonParticle.hpp"
#include "muo_evt/MuonQualityInfo.hpp"
#include "muo_evt/MuTmbObj.hpp"

#include "jet_evt/JetChunk.hpp"
#include "jet_evt/JetAlgoInfo.hpp"
#include "jet_evt/Jet.hpp"

#include "AATrack/AA.hpp"
#include "AATrack/Ana.hpp"
#include "AATrack/Trk.hpp"
#include "AATrack/Ptl.hpp"
#include "AATrack/V0Finder.hpp"


#include "missingET/MissingET.hpp"
#include "missingET/MissingETChunk.hpp"

#include "vertex_evt/VertexCollChunk.hpp"
#include "vertexutil/Vertex.hpp"

#include "l1l2_evt/L1L2Chunk.hpp"

#include "l3fchunk/L3Chunk.hpp"
#include "l3fchunk/L3ChunkSelector.hpp"
#include "l3fresults/L3PhysicsResults.hpp"
#include "l3femtools/L3ElePhysicsResults.hpp"
#include "l3fMuonTools/L3MuonPhysicsResults.hpp"
#include "l3fjettools/L3JetsPhysicsResults.hpp"
#include "l3fbtag_ip/L3BTagIPPhysicsResults.hpp"
#include "l3fJetMEt/L3MEtPhysicsResults.hpp"
#include "l3femtools/L3PhotonPhysicsResults.hpp"
#include "l3fTauTools/L3TauPhysicsResults.hpp"
#include "l3fJetMEt/L3HtPhysicsResults.hpp"

#include <iostream>
#include <iomanip>
#include <vector>

using namespace edm;
using namespace std;
using namespace fwk;
using namespace emid;
using namespace muonid;
using namespace jetid;
using namespace trf;
using namespace vertex;
using namespace D0analysis;

// The following lines are here for educational purposes only.

#include "ZMutility/iostream"
#ifndef ZMEXCEPTION_H
#include "Exceptions/ZMexception.h"
#endif // ZMEXCEPTION_H
#include "Exceptions/ZMthrow.h"

ZM_USING_NAMESPACE (zmex); 
ZMexStandardDefinition(ZMexception,ZMxD0analysisElvisIsAlive);
ZMexClassInfo ZMxD0analysisElvisIsAlive::_classInfo("ZMxD0analysis","D0ChunkAnalyze",ZMexERROR );

///////////////////////////////////////////////////////////////////////////////

namespace D0example {

  FWK_REGISTRY_IMPL(D0ChunkAnalyze,"$Name:  $")

  D0ChunkAnalyze::D0ChunkAnalyze(Context* context): 
        Package(context),    Process(context),    Analyze(context), 
        FileOpen(context),   FileClose(context), JobSummary(context) {
  // Default constructor.

    cout << endl
         << "+------------------------------------------+" << endl
         << "| D0ChunkAnalyze - package initialization  |" << endl
         << "+------------------------------------------+" << endl
         << endl;

    // Access the RCP information.
    edm::RCP rcp = packageRCP();


   // Select corrected Muon chunk after MuoCandidate
   edm::RCP muonidAlgoRCP = rcp.getRCP("Muonid_Algo");
   muAlgoSel = MuonParticleChunkSelector(muonidAlgoRCP);

  
   // Select one of the Jet chunks.
    std::string jetType               = rcp.getString("JetAlgo_type");
    std::vector<float> jetValues      = rcp.getVFloat("JetAlgo_values");
    std::vector<std::string> jetNames = rcp.getVString("JetAlgo_names");

    if ( jetNames.size() > 0 )
      while ( *jetNames.rbegin() == "" ) jetNames.pop_back(); 

    if ( jetValues.size() < jetNames.size() ) {
      cout << "D0ChunkAnalyze: inconsistent RCP values for the jet selection, disable it" << endl;
      jetNames.clear();
      jetValues.clear();
    } 

    JetAlgoInfo jetAlgoInfo(jetType,jetValues,jetNames);
    cout << "D0ChunkAnalyze: select jets matching the following criteria: " << endl
         << jetAlgoInfo << endl;

    _jetCSel = JetChunkSelector(jetAlgoInfo);

   // Name of the tag assigned to selected events.
    _event_tag = rcp.getString("Event_tag");
  
   // Initialise event counters.
    _processed_events = 0;
    _triggered_events = 0;
    _selected_events = 0;

   // Trigger selection.
    _selectTriggers = D0TriggerDecoder::Instance();
  

   // Initialize ROOT.


     // Open the ROOT file which will contain the Tree/Histograms.
      string histoFile = rcp.getString("HistogramFile");
      cout << "D0ChunkAnalyze: Tree/Histos will be written to: " << histoFile << endl;

      Int_t nevent = 400;     // by default create 400 events
      Int_t comp   = 1;       // by default file is compressed
      Int_t split  = 1;       // by default, split Event in sub branches
      Int_t bufsize = 256000; // buffer size

      phfile = new TFile("Mujet.root","RECREATE","TTree ROOT file");
      phfile->SetCompressionLevel(comp);


      // Create a simple ROOT tuple.
      ptree = new TTree("Tree","ROOT tuple with muons jets and vertices");
      ptree->SetAutoSave(1000000000);  // autosave when 1 Gbyte written

//      if (split)  bufsize /= 4;
//      TBranch *b = ptree->Branch("event", "Event", &event,  bufsize,split);

     // Book Event Block
      ptree->Branch("nrun",&nrun,"nrun/I");
      ptree->Branch("nevt",&nevt,"nevt/I");

     // Book Muon Block
      ptree->Branch("nmu",&nmu,"nmu/I");      
      ptree->Branch("mu_index",&mu_index,"mu_index[nmu]/I");
      ptree->Branch("mu_pt",&mu_pt,"mu_pt[nmu]/F");
      ptree->Branch("mu_eta",&mu_eta,"mu_eta[nmu]/F");
      ptree->Branch("mu_phi",&mu_phi,"mu_phi[nmu]/F");
      ptree->Branch("mu_e",&mu_e,"mu_e[nmu]/F");
      ptree->Branch("mu_px",&mu_px,"mu_px[nmu]/F");
      ptree->Branch("mu_py",&mu_py,"mu_py[nmu]/F");
      ptree->Branch("mu_pz",&mu_pz,"mu_pz[nmu]/F");
      ptree->Branch("mu_ip",&mu_ip,"mu_ip[nmu]/F"); 
      ptree->Branch("mu_iperr",&mu_iperr,"mu_iperr[nmu]/F");
      ptree->Branch("mu_zpca",&mu_zpca,"mu_zpca[nmu]/F");
      ptree->Branch("mu_tanlam",&mu_tanlam,"mu_tanlam[nmu]/F");
      ptree->Branch("mu_trklnk",&mu_trklnk,"mu_trklnk[nmu]/I");
      ptree->Branch("mu_vtxlnk",&mu_vtxlnk,"mu_vtxlnk[nmu]/I");
      ptree->Branch("mu_whitsA",&mu_whitsA,"mu_whitsA[nmu]/I");
      ptree->Branch("mu_whitsBC",&mu_whitsBC,"mu_whitsBC[nmu]/I");
      ptree->Branch("mu_shitsA",&mu_shitsA,"mu_shitsA[nmu]/I");
      ptree->Branch("mu_shitsBC",&mu_shitsBC,"mu_shitsBC[nmu]/I");
      ptree->Branch("mu_sctimeA",&mu_sctimeA,"mu_sctimeA[nmu]/F");
      ptree->Branch("mu_sctimeB",&mu_sctimeB,"mu_sctimeB[nmu]/F");
      ptree->Branch("mu_sctimeC",&mu_sctimeC,"mu_sctimeC[nmu]/F");
      ptree->Branch("mu_drjet5",&mu_drjet5,"mu_drjet5[nmu]/F");
      ptree->Branch("mu_ntrk5",&mu_ntrk5,"mu_ntrk5[nmu]/I");
      ptree->Branch("mu_etrk5",&mu_etrk5,"mu_etrk5[nmu]/F");
      ptree->Branch("mu_ethalo",&mu_ethalo,"mu_ethalo[nmu]/F");
      ptree->Branch("mu_octant",&mu_octant,"mu_octant[nmu]/I");
      ptree->Branch("mu_region",&mu_region,"mu_region[nmu]/I");
      ptree->Branch("mu_chi2",&mu_chi2,"mu_chi2[nmu]/F");
      ptree->Branch("mu_chi2loc",&mu_chi2loc,"mu_chi2loc[nmu]/F");
      ptree->Branch("mu_caleta",&mu_caleta,"mu_caleta[nmu]/F");
      ptree->Branch("mu_calphi",&mu_calphi,"mu_calphi[nmu]/F");
      ptree->Branch("mu_e33",&mu_e33,"mu_e33[nmu]/F");
      ptree->Branch("mu_e55",&mu_e55,"mu_e55[nmu]/F");   
      ptree->Branch("mu_econe1",&mu_econe1,"mu_econe1[nmu]/F");
      ptree->Branch("mu_econe15",&mu_econe15,"mu_econe15[nmu]/F");
      ptree->Branch("mu_econe2",&mu_econe2,"mu_econe2[nmu]/F");
      ptree->Branch("mu_econe4",&mu_econe4,"mu_econe4[nmu]/F");
      ptree->Branch("mu_econe6",&mu_econe6,"mu_econe6[nmu]/F");


     // Book Jet Block
      ptree->Branch("njet",&njet,"njet/I");
      ptree->Branch("jet_et",&jet_et,"jet_et[njet]/F");
      ptree->Branch("jet_eta",&jet_eta,"jet_eta[njet]/F");
      ptree->Branch("jet_phi",&jet_phi,"jet_phi[njet]/F");
      ptree->Branch("jet_px",&jet_px,"jet_px[njet]/F");
      ptree->Branch("jet_py",&jet_py,"jet_py[njet]/F");
      ptree->Branch("jet_pz",&jet_pz,"jet_pz[njet]/F");
      ptree->Branch("jet_n90",&jet_n90,"jet_n90[njet]/I");
      ptree->Branch("jet_emf",&jet_emf,"jet_emf[njet]/F");
      ptree->Branch("jet_chf",&jet_chf,"jet_chf[njet]/F");
      ptree->Branch("jet_hcr",&jet_hcr,"jet_hcr[njet]/F");

     // Book Vertex Block
      ptree->Branch("nvtx",&nvtx,"nvtx/I");
      ptree->Branch("vtx_x",&vtx_x,"vtx_x[nvtx]/F");
      ptree->Branch("vtx_y",&vtx_y,"vtx_y[nvtx]/F");
      ptree->Branch("vtx_z",&vtx_z,"vtx_z[nvtx]/F");
      ptree->Branch("vtx_chi2",&vtx_chi2,"vtx_chi2[nvtx]/F");
      ptree->Branch("vtx_ndof",&vtx_ndof,"vtx_ndof[nvtx]/I");


    // Finish initialization.
    cout << endl
         << "+------------------------------------------+" << endl
         << "| D0ChunkAnalyze - initialization finished |" << endl
         << "+------------------------------------------+" << endl
         << endl;
  
  }

  Result D0ChunkAnalyze::processEvent(Event &event) {
  // Process one event.
  
    // Increase the counter for the number of events processed.
    ++_processed_events;
  
    // Check whether the event satisfies the user trigger requirements
    // (assume we want to count the events in which the EM_HI trigger fires).
    if ( _selectTriggers->TriggerFired("EM_HI") ) ++_triggered_events;

    return Result::success;
  
  }

  Result D0ChunkAnalyze::analyzeEvent(const Event &event) {
  
    // Run and event number.
    nevt = static_cast<int> (event.collisionID().eventNumber());
    nrun = static_cast<int> (event.collisionID().runNumber());

    // Get the D0reco processing version, extracting it from the history chunk.
    HistorySelectorLatest hs("D0reco");
    edm::TKey<HistoryChunk> key(hs);
    THandle<HistoryChunk> history=key.find(event);
    string vrec = history->prod_version();


    // ---------------------------------------------------------------------------------
    // Access the muon particle chunk
    // ---------------------------------------------------------------------------------


    // Select the certified MuonParticleChunk ( has to call muocand before analyze)
    const edm::TKey<muonid::MuonParticleChunk> mupKey(muAlgoSel);
    THandle<MuonParticleChunk> muChunk=mupKey.find(event);


    // Check whether a non empty muon particle chunk was found
    int numMuo = 0;
    if ( !muChunk.isValid() ) {
      cout << "SHIT !!! no muchunk pointer !!!! " << endl;
    }
    else{

      vector<MuonParticle>::const_iterator mupptr;
      const vector<MuonParticle>* pvector=muChunk->getParticles();

      // Number of muons
      numMuo = pvector->size();

      nmu=0; // initialize Rootree muon index

      // Loop on the muon candidates.
      for ( mupptr=pvector->begin(); mupptr!=pvector->end(); ++mupptr ) {

        const muonid::MuonQualityInfo* muquality = mupptr->qualInfo(); // pointer to the quality information

        // Muon quality words
        int nSegMu = muquality->nseg();  // number of track segments (+ for central track matches)
        int nHitMu = muquality->nhit();  // packed word containing the hits in the wires and scintillators
        float chiSqMu = muquality->chisqloc();  // chi square of the local muon fit

      
        // Decode # hits informations
        int nWScHit[5];
        for ( int ihit=0; ihit!=5; ++ihit ) 
          nWScHit[ihit] = 0;
        int ii = 0;
        while ( nHitMu>0 ) {      // decode the information on the number of hits in
          nWScHit[ii] = nHitMu%10;                        // the A layer scintillator and wires hits, and in the 
          nHitMu /= 10;                                   // B+C layers scintillator and wires hits:
          ++ii;                                           
        } // Decoding muon quality information
        int nAWHit = nWScHit[0];                          // A layer wires
        int nBCWHit = nWScHit[1] + nWScHit[2]*10;         // B+C layers wires
        int nASHit = nWScHit[3];                          // A layer scintillators
        int nBCSHit = nWScHit[4];                         // B+C layers scintillators


        // Apply the tight muon selection cuts
//        if ( nWScHit[0]>=2    && nBCWHit>=2 && 
//             nWScHit[3]>=1 && nWScHit[4]>=1 && chiSqMu>0. ) {

          float globalpT = mupptr->pT();                  // best measurement of the muon pT
          globalpT = globalpT*mupptr->charge();           // sign the muon pT with the charge of the muon
          
          // For muons matched to central tracks compare the two possible muon momentum measurements
          float localpT= 0.0;
          if ( nSegMu>0 ) {

            localpT = muquality->qptloc();       // local measurement of the muon pT
          } // Muons matched to central tracks

             float ip = 0.0;
             float iperr = 0.0;
             mupptr->impPar( ip , iperr ); // Impact parameter


          // Fill muon tree buffers
  


//          mu_index[nmu] = mupptr->index();
          mu_index[nmu] = 0;

          mu_pt[nmu] = globalpT;
          mu_eta[nmu]= mupptr->eta();
          mu_phi[nmu]= mupptr->phi();
          mu_e[nmu]  = mupptr->E();
          mu_px[nmu] = mupptr->px();
          mu_py[nmu] = mupptr->py();
          mu_pz[nmu] = mupptr->pz();
          mu_ip[nmu]    = ip;
          mu_iperr[nmu] = iperr;
          mu_zpca[nmu]  = mupptr->zAtPca();
          mu_tanlam[nmu] = mupptr->tanlam();

//          mu_trklnk[nmu] = mupptr->vtxindex();
//          mu_vtxlnk[nmu] = mupptr->chpLink();
//          mu_local[nmu] = muquality->hasLocal();
//          mu_central[nmu] = muquality->hasCentral();
//          mu_cal[nmu] = muquality->hasCal();
//          mu_axial[nmu] = muquality->isAxialMatched();

          mu_trklnk[nmu] = 0;
          mu_vtxlnk[nmu] = 0;
          mu_whitsA[nmu]  = muquality->whitsA();  // # of wire hits
          mu_whitsBC[nmu] = muquality->whitsBC();
          mu_shitsA[nmu]  = muquality->shitsA();  // # of scint hits
          mu_shitsBC[nmu] = muquality->shitsBC();
          mu_sctimeA[nmu] = muquality->sctimeA(); //
          mu_sctimeB[nmu] = muquality->sctimeB();
          mu_sctimeC[nmu] = muquality->sctimeC();
          mu_octant[nmu]  = muquality->octant();
          mu_region[nmu]  = muquality->region();
          mu_chi2[nmu]    = muquality->chisq();
          mu_chi2loc[nmu] = muquality->chisqloc();
          mu_drjet5[nmu]  = muquality->drJet5();
          mu_ntrk5[nmu]   = muquality->nTrk5();
          mu_etrk5[nmu]   = muquality->etTrkCone5();
          mu_ethalo[nmu]  = muquality->etHalo();
          mu_caleta[nmu]  = muquality->CalEta();
          mu_calphi[nmu]  = muquality->CalPhi();
          mu_e33[nmu]     = muquality->e33();
          mu_e55[nmu]     = muquality->e55();
          mu_econe1[nmu]  = muquality->EInCone1();
          mu_econe15[nmu] = muquality->EInCone15();
          mu_econe2[nmu]  = muquality->EInCone2();
          mu_econe4[nmu]  = muquality->EInCone4();
          mu_econe6[nmu]  = muquality->EInCone6();

          nmu++; // increment muon index

          // Give up the event if there are too many muons
          if ( nmu == 20 ) return Result::success;

//        } // Tight muon selection

      } // Loop on muon candidates

    } // Valid muon particle chunk

    // Give up the event if there are no muons
    if ( nmu == 0 ) return Result::success;
    


    // ---------------------------------------------------------------------------------
    // Access the jet particle chunk
    // ---------------------------------------------------------------------------------
    const edm::TKey<JetChunk> jetKey(_jetCSel);
    edm::THandle<JetChunk> jhandle = jetKey.find(event);

    njet=0;


    // Check whether a non empty jet particle chunk was found
    int numJet = 0;
    if ( jhandle.isValid() ) {

      JetChunk* jetChunk=jhandle.ptr();
      const vector<Jet>* pjetvec = jetChunk->getParticles();

      // Number of jets
      numJet = pjetvec->size();

      // Loop over jets
      vector<Jet>::const_iterator itj;
      list<const Jet*> listjet;        // list of pointers to jets
      for ( itj = pjetvec->begin() ; itj != pjetvec->end() ; ++itj ){
        listjet.push_back(&(*itj));
      }


      // Loop on all the jets and copy the pointers into a list
//      list<const Jet*> listjet;
//      for ( int i=0; i!=numJet; ++i ) {
//        const Jet* jetptr = &(*pjetvec)[i];
//        listjet.push_back(jetptr);
//      }

      if ( numJet) {
        // Sort the jets in order of decreasing Pt
        listjet.sort(JetPtOrder());

        // Apply some simple quality cuts on jets
        list<const Jet*>::iterator ijet; 
        for( ijet = listjet.begin() ; ijet != listjet.end() ; ++ijet ){

          int jetN90 = (*ijet)->n90();   // # of towers containing 90% of the energy         
          float jetEMfrac = (*ijet)->emETfraction();  // fraction of the energy in the EM calorimeter
          if ( jetN90>0 && jetEMfrac<0.95 ) {
            
            jet_et[njet] = (*ijet)->pT();     // jet uncorrected Et
            jet_eta[njet] = (*ijet)->eta();
            jet_phi[njet] = (*ijet)->phi();
            jet_px[njet]  = (*ijet)->px();
            jet_py[njet]  = (*ijet)->py();
            jet_pz[njet]  = (*ijet)->pz();
            jet_n90[njet] = jetN90;
            jet_emf[njet] = jetEMfrac;
            jet_chf[njet] = (*ijet)->chETfraction();
            jet_hcr[njet] = (*ijet)->hotcellratio();

            njet++; // increment jet index

          } // Jet quality requirements

        } // Jets list

      } // Jet(s) found

    } // Valid jet chunk

    // ---------------------------------------------------------------------------------
    // Access the calorimeter chunk chunk for Rap Gap info
    // ---------------------------------------------------------------------------------



// Eta range
float eta1=2.6;
float eta2=5.3;

// Cell energy thresholds (GeV)
float EM_cellcut=0.1;
float FH_cellcut=0.2;

  for (int i=0; i<17; i++) {
     if (i<=6) cellcut[i] = EM_cellcut;
     else if (i>=10 && i<=13) cellcut[i] = FH_cellcut;
     else cellcut[i] = 0.0;
   }


// Log (cell energy sum) threshold
float E_sumcut=1.0;




  //__ ACCESS CHUNKS ____________________________
  
  // Access L1L2 Chunk
  edm::TKey<L1L2Chunk> const l1l2Key;
  edm::THandle<L1L2Chunk> l1l2Chunk = l1l2Key.find(event);
  
  // Access CalDataChunk
  edm::TKey<CalDataChunk> const calKey;
  edm::THandle<CalDataChunk> calChunk = calKey.find(event);
  
  // Exit if bad chunks
  if (!l1l2Chunk.isValid() || !calChunk.isValid()) 
    return Result::success;

  //__ LUMINOSITY MONITOR _______________________
 
  // To require that the North LM is off:
  if (!and_or(event,221) 
      && !and_or(event,217) 
      && !and_or(event,222) && !and_or(event,223)) {
    // LM North gap tag
  }
  
  // To require that the South LM is off:
  if (!and_or(event,220) 
      && !and_or(event,217) && !and_or(event,222) 
      && !and_or(event,223)) {
    // LM South gap tag
  }
  
  //__ CALORIMETER ___________________________

  // Get vertex of event
  float vtx[3]={primary_vertex(event,0),
primary_vertex(event,1),
primary_vertex(event,2)};
  
  // Calorimeter cells
  list<const CalCell *> cell_list;
  CalData* calData = calChunk->getCalData(vtx);
  calData->getCells(cell_list);
  
  // Reject if no cells in readout
  if (cell_list.empty()) return Result::success;
  
  // Loop through cells in candidate events
  std::list<const CalCell*>::const_iterator icell;
  for (icell = cell_list.begin();icell != cell_list.end();icell++) {
    
    int celllayer = (*icell)->address().layer();
    float cellenergy = (*icell)->E();
    float celleta = (*icell)->eta();
    
    int forward = 0;
    int cell_depth = 0;
    int cell_end = 0;
    
    // Categorise cell
    if (abs(celleta) > eta1 && abs(celleta) < eta2) 
      forward = 1; // FORWARD
    if (celllayer >= 1 && celllayer <= 7) cell_depth = 1; // EM
    if (celllayer >= 11 && celllayer <= 14) cell_depth = 2; // FH
    if (celleta < 0) cell_end = 1; // NORTH
    if (celleta > 0) cell_end = 2; // SOUTH
    
    //- CELL ENERGY SUM 
    
    // Adding cells to energy sum
    if (forward && cell_depth && cell_end) {
      // Apply energy threshold  
      if (cellenergy > cellcut[celllayer-1]) {
cell_energy_sum[cell_depth-1][cell_end-1] += cellenergy;
      }
    } 
  } // cell loop
  
  //__ ENERGY SUM CUT _____________________
  
  // North
  if (log10(cell_energy_sum[0][0] + cell_energy_sum[1][0]) 
      < Esumcut) {
    // Esum North gap tag
  }
  
  // South
  if (log10(cell_energy_sum[0][1] + cell_energy_sum[1][1])
      < E_sumcut) {
    // Esum South gap tag
  }
  
  // Gap North event if LM North tag + Esum North tag

  // Gap South event if LM South tag + Esum South tag


// Function to access L1 and/or terms
bool and_or(const Event& event, const int& number)
{
  bool and_or=0;
  int number2=number;
  edm::TKey<L1L2Chunk> const l1l2Key;
  edm::THandle<L1L2Chunk> l1l2Chunk = l1l2Key.find(event);
  if ( l1l2Chunk.isValid() ) {
    l1l2info trigmap = (l1l2Chunk.ptr())->Ret_l1l2info();
    and_or = trigmap.l1andor(number2);
  }    
  return and_or;
}



    // ---------------------------------------------------------------------------------
    // Access the track chunk
    // ---------------------------------------------------------------------------------
//        // Loop over muons to find a matching track
//        int imu;
//        float PI = 3.14159265; 
//        float drmin = 999.0;
//        for (imu = 0 ; imu <= nmu ; ++imu ){
//          float dphi = abs( phi - mu_phi[imu] );
//          if ( dphi > (PI/2) ) dphi = ( PI - dphi );
//          float deta = ( eta - mu_eta[imu] );
//          float dr = sqrt( dphi*dphi + deta*deta );
//          if ( drmin > dr ) {
//            drmin = dr;
//            mu_index= imu
//          }
//        }
//        
//   cout << "DRMIN=" << drmin;

       
    //------------------------------------------------------------------------
    //  Access AATrack vertex information
    //------------------------------------------------------------------------

    //    params();
    //    AA::displayParameters();

  //define beam-spot and field value
    //    AA::setField();
    //    AA::spot.set();

 if(AA::trkBox.tracks()->size() > 400) continue;

 //refit of tracks
    //  for(AA::TrkLstIt p = AA::trkBox.tracks()->begin(); 
    //	 p != AA::trkBox.tracks()->end(); ++p){
    //	    if(!(*p)->valid()) continue;
    //	    (*p)->fit(true);
    //	    (*p)->debug();
    //      }

 //Make vertex of these 3 particles
   //    TrkLst lstT;
   //    lstT.push_back(ptl1);
   //    lstT.push_back(ptl2);
   //    lstT.push_back(ptl3);
   //    Vrt vrt;
   //    if(!vrt.fill(ptl->decayVertex()->x(),&lstT)) continue;
   //    if(!vrt.filter()) continue;
   //    if(vrt.size() != 3) continue;

 

  // Get the list of primary VTX
  const VrtLst* pvl = AA::vrtPBox.vertices();

  // Loop over list of primary VTX
  for(VrtLstCIt p = pvl->begin(); p != pvl->end(); ++p){  
  
    // Get the list of secundary VTX associated to a given primary VTX 
    const SVLst* plsv = AA::svFinder.getSV(*p); 
    if(plsv == 0) continue;
    
    // Loop over list of secundary VTX
    for(SVLstCIt psv = plsv->begin(); psv != plsv->end(); ++psv){  

      // Get the list of particles associated to a given secundary VTX
      const ChainPLst* pcl = psv->ptl->children();
      
      // Loop over particles attached to a given secundary VTX
      for(ChainPLstCIt pc1 = pcl->begin(); pc1 != pcl->end(); ++pc1){
	Ptl* ptl1 = AA::particle((*pc1)->track());
	for(ChainPLstCIt pc2 = pc1; pc2 != pcl->end(); ++pc2){
	  if(pc2 == pc1) continue;
	  Ptl* ptl2 = AA::particle((*pc2)->track());
	  
	  // Exclude particles with large impact parameter error ???
	  if(ptl1->imp(1)*ptl1->imp(1) < 9.*ptl1->vimp(1,1) &&
             ptl2->imp(1)*ptl2->imp(1) < 9.*ptl2->vimp(1,1)) continue;

          // Exlude neutral particles
	  if(ptl1->q() * ptl2->q() >= 0) continue;

          // Calculate the hypothesys 2body invariant mass 
	  HepVector mom(3);
	  double mmm = AA::xMass(psv->pve,*pc1,*pc2,AA::MU_PLUS,AA::MU_MINUS,mom);

	  //Try to find some B+
	  if(mmm > 2.90 && mmm < 3.25) {
	    
	    // Loop over a 3rd particle attached to the same VTX to get the B+ mass
	    for(ChainPLstCIt pc3 = pcl->begin(); pc3 != pcl->end(); ++pc3){
	      if(pc3 == pc1) continue;
	      if(pc3 == pc2) continue;
	      Ptl* ptl3 = AA::particle((*pc3)->track());

	      ChainPLst lstc;
	      lstc.push_back(*pc1);
	      lstc.push_back(*pc2);
	      lstc.push_back(*pc3);
	      
	      // Calculate hypotheses B mass
	      double mb = AA::xMass(psv->pve,&lstc,&types) - mmm + 3.097;

              // Impose minmum Pt treshold
	      if(ptl1->pt()>1.5 && ptl2->pt()>1.5 && ptl3->pt()>0.5){

                 // Fiil histos
	         HFILL(1630,mb,0.,1.);

              }

	    }  // end loop over 3rd particle


	  } // end if over B+ search
	  	  
	} // end loop over 2nd particle

      } // end loop over 1st particle
        
    } // end loop over secundary VTX

  } // end loop over primary VTX





// JPsi studies

 vector<PType> types;
  types.push_back(AA::MU_PLUS);
  types.push_back(AA::MU_MINUS);
  types.push_back(AA::K_PLUS);

  //Get list of all J/psi candidates
  const PtlLst* plst = AA::jPsiFinder.particles();

  //Make loop and analyze them
  for(AA::PtlLstCIt p = plst->begin(); p != plst->end(); ++p){
    AA::Ptl* ptl = *p;

    //Selection of J/psi signal
    //    if(ptl->ptot() < 5.) continue;

    //decay length for J/psi and its mass
    double dxy, vxy;
    ptl->decayVertex()->distanceXY(ptl->primaryVertex(),dxy,vxy);
    double cxy = ptl->cosXY();
    double mjpsi = ptl->mass(AA::MU_PLUS, AA::MU_MINUS);

    //Pointers to 2 muons from J/psi
    const list<AA::Chain*>* plsc = ptl->children();
    AA::Ptl* ptl1 = AA::particle(plsc->front()->track());
    AA::Ptl* ptl2 = AA::particle(plsc->back()->track());

    if(ptl1->pt() > 2.0 && ptl2->pt() > 2.0) {
      if(mjpsi>2.95 && mjpsi < 3.25 && ptl1->nSMT()>0 && ptl2->nSMT()>0) {
	if(cxy >= 0.) HFILL(1405,dxy,0.,1.);
	else          HFILL(1405,-dxy,0.,1.);
	if(cxy >= 0.) HFILL(1406,dxy/sqrt(vxy),0.,1.);
	else          HFILL(1406,-dxy/sqrt(vxy),0.,1.);
	if(cxy >= 0.9)       HFILL(1407,dxy/sqrt(vxy),0.,1.);
	else if(cxy <= -0.9) HFILL(1407,-dxy/sqrt(vxy),0.,1.);
      }
    }

    //Selected signal of J/psi
    HFILL(2001,mjpsi,0.,1.);

    //Special selections for high mass mu-mu pairs
    if(mjpsi > 20.) {
      MuoGlb* pmu1 = ptl1->muon();
      MuoGlb* pmu2 = ptl2->muon();
      bool ok = true;
      if(pmu1 != 0 && pmu2 != 0) {
	if(abs(pmu1->nseg())<1 || abs(pmu1->nseg()) < 1) ok = false;
	if(ptl1->pt()<2.5 || ptl2->pt()<2.5) ok = false;
	if(ok) HFILL(2000,mjpsi,0.,1.);
      }
    }

    //J/psi with CFT-only tracks
    if(ptl1->nSMT() < 2 || ptl2->nSMT() < 2) HFILL(2004,mjpsi,0.,1.);

    //Mass constraint on J/psi
    //dp contains corrections of momenta from the mass constraint
    HepVector dp(5);
    if(mjpsi<2.80 || mjpsi>3.35) continue;
    if(!ptl->massConstraint(AA::JPSI,AA::MU_PLUS,AA::MU_MINUS,dp)) dp *= 0.;
    dp(3) = 0.;
    dp(4) = 0.;
    dp(5) = 0.;

    //List of all particles in event
    PtlLst lst = *AA::ptlBox.particles();

    //Remove particles identified as products of K0, Lambda, gamma decays
    v0Finder.filter(lst);

    findBToJK(ptl,&lst,dp);
    findBToJKStar(ptl,&lst,dp);
    findBToJK0(ptl,dp);
  }



//Search for decay B+ -> J/psi K+

  //Pointers to 2 muons from J/psi -> mu+ mu-
  const list<AA::Chain*>* plsc = ptl->children();
  Ptl* ptl1 = AA::particle(plsc->front()->track());
  Ptl* ptl2 = AA::particle(plsc->back()->track());

  //types of analysed particles (will be used for mass computation)
  vector<PType> types;
  types.push_back(AA::MU_PLUS);
  types.push_back(AA::MU_MINUS);
  types.push_back(AA::K_PLUS);

  //Loop over all particles
  for(PtlLstCIt p = plst->begin(); p != plst->end(); ++p){
    Ptl* ptl3 = *p;
    if(ptl3 == ptl1) continue;
    if(ptl3 == ptl2) continue;

    //Cut on pt of track
    if(ptl3->pt() < 0.5) continue;

    //At least 2 tracks with SMT measurements 
    //(required for precise vertex measurement)
    int count = 0;
    if(ptl1->nSMT() < 2) ++count;
    if(ptl2->nSMT() < 2) ++count;
    if(ptl3->nSMT() < 2) ++count;
    if(count > 1) continue;

    //Make vertex of these 3 particles
    TrkLst lstT;
    lstT.push_back(ptl1);
    lstT.push_back(ptl2);
    lstT.push_back(ptl3);
    Vrt vrt;
    if(!vrt.fill(ptl->decayVertex()->x(),&lstT)) continue;
    if(!vrt.filter()) continue;
    if(vrt.size() != 3) continue;

    //Find primary vertex with minimal distance to JK vertex
    Vrt* pv = 0;
    if(!vrt.associateToVrt(&AA::vrtPBox,pv))    continue;
    double lv,vlv;
    vrt.distance(pv,lv,vlv);
    double lxy, vlxy;
    vrt.distanceXY(pv,lxy,vlxy);

    //Construct new particle (B+) from these 3 tracks
    int q = ptl1->q() + ptl2->q() + ptl3->q();
    Ptl jk;
    if(!jk.combine(&vrt,q)) continue;
    if(!jk.associateToVrt(&AA::vrtPBox)) continue;
    if(jk.chi2Vrt() > 64.) continue;

    //compute mass of this new particle and angle between momentum and direction
    jk.fillMomentumCorrection(dp);
    double mjk = jk.massCorrected(types);
    double cxyjk = jk.cosXY();

    //Compute mass (jpsi - pi) (for Bc search)
    types.back() = AA::PI_PLUS;
    double mjpi = jk.massCorrected(types);
    types.back() = AA::K_PLUS;

    //Final selections (should be similar to bc.f)
    if(ptl3->pt() < 0.5) continue;
    if(ptl3->ptot() < 0.7) continue;
    if(ptl->ptot() < 5.0) continue;
    if(cxyjk < 0.95) continue;
    if(lxy * lxy < vlxy*20.25) continue; //20.25 = 4.5*4.5
    if(ptl3->missForw()->size() > 2) continue;
    if(ptl3->nSMT() < 2) continue;
    if(ptl3->chi2Vrt() < 9.) continue;
    if(vrt.chi2() > 16.) continue;
    if(jk.chi2Vrt() > 49.) continue;
    if(ptl3->pt() < 1. ) {
      if(lxy*lxy < vlxy*30.25) continue; //30.25 = 5.5*5.5
      if(vrt.chi2() > 9.) continue;
      if(ptl3->jet() != ptl1->jet() && ptl3->jet() != ptl2->jet()) continue;
    }

    //Resulting Mass spectrum
    HFILL(2002,mjk,0.,1.);
    HFILL(2003,mjpi,0.,1.);
    HFILL(2009,mjpi,0.,1.);

  }





    // Fill the ROOT tuple




//_______________________________________________________________________
//
//   Fill the Root Tree
//_______________________________________________________________________    



    // Fill all ROOT Tree
    ptree->Fill();

    return Result::success;
  }
  
  Result D0ChunkAnalyze::jobSummary() {
  // Print a summary page and save the histograms.
    cout << endl
         << "+------------------------------------------+" << endl
         << "| D0ChunkAnalyze - package termination     |" << endl
         << "+------------------------------------------+" << endl
         << "| Number of events processed " << setw(13) << _processed_events << " |" << endl
         << "| Number of events triggered " << setw(13) << _triggered_events << " |" << endl
         << "| Number of events selected  " << setw(13) << _selected_events << " |" << endl
         << "+------------------------------------------+" << endl;
  
    SaveHistograms();
  
    cout << "| Histograms saved                         |" << endl
         << "+------------------------------------------+" << endl;
  
    return Result::success;
  }
  
  D0ChunkAnalyze::~D0ChunkAnalyze() {
  // Default destructor, ensure that histograms have been saved.
//    if ( !_histo_saved ) SaveHistograms();
  }

  void D0ChunkAnalyze::SaveHistograms() {
  // Save the ROOT Tree and Histograms.

  cout << "D0ChunkAnalyze::SaveHistograms  closing file Mujet.root" << endl;
  ptree->Print();
  phfile->Write();

  }

  Result D0ChunkAnalyze::fileOpen(const FileInfo &fileinfo) {
  // Executed whenever a new input file is opened
    cout << "D0ChunkAnalyze::fileOpen opening file " << fileinfo.filename() << endl;
    return Result::success;
  }

  Result D0ChunkAnalyze::fileClose(const FileInfo &fileinfo) {
  // Executed whenever an input file is closed
    cout << "D0ChunkAnalyze::fileOpen closing file " << fileinfo.filename() << endl;
    return Result::success;
  }

} // namespace D0example