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

void UEAnalysisMPI::mpiAnalysisMC(Float_t weight,Float_t etaRegion,Float_t ptThreshold, TClonesArray& ChargedJet)
{
  std::vector<TLorentzVector*> JetMC;
  JetMC.clear();
  
  for(int j=0;j<ChargedJet.GetSize();++j){
    TLorentzVector *v = (TLorentzVector*)ChargedJet.At(j);
    if(fabs(v->Eta())<etaRegion){
      JetMC.push_back(v);
    }
  }
  
  std::vector<AssociatedObject> assoJetMC;
  assoJetMC.clear();

  while(JetMC.size()>1){
    int oldSize = JetMC.size();
    std::vector<TLorentzVector*>::iterator itH = JetMC.begin();
    if((*itH)->Pt()>=ptThreshold){
      for(std::vector<TLorentzVector*>::iterator it=JetMC.begin();it!=JetMC.end();it++){
	float azimuthDistanceJet = fabs( (*itH)->Phi() - (*it)->Phi() );
	if((*it)->Pt()/(*itH)->Pt()>=0.3){
	  if( (piG - rangePhi) <  azimuthDistanceJet && azimuthDistanceJet < (piG + rangePhi)) {
	    AssociatedObject tmpPair((*itH),(*it));
	    assoJetMC.push_back(tmpPair);
	    JetMC.erase(it);
	    int newSize = oldSize -1;
	    oldSize = newSize;
	    JetMC.resize(newSize);
	    break;
	  }
	}
      }
    }
    JetMC.erase(itH);
    int newSize = oldSize -1;
    JetMC.resize(newSize);
  }
  
  if(assoJetMC.size()){
    fNumbMPIMC->Fill(assoJetMC.size());
    std::vector<AssociatedObject>::iterator at= assoJetMC.begin();
    
    const TLorentzVector* leadingJet((*at).first);
    const TLorentzVector* secondJet((*at).second);

    pPtRatio_vs_PtJleadMC->Fill(leadingJet->Pt(),(secondJet->Pt()/leadingJet->Pt()));
    pPtRatio_vs_EtaJleadMC->Fill(fabs(leadingJet->Eta()),(secondJet->Pt()/leadingJet->Pt()));
    pPtRatio_vs_PhiJleadMC->Fill(leadingJet->Phi(),(secondJet->Pt()/leadingJet->Pt()));
    
    fdEtaLeadingPairMC->Fill(leadingJet->Eta()-secondJet->Eta());
    float dPhiJet = fabs(leadingJet->Phi()-secondJet->Phi());
    if(dPhiJet> piG) dPhiJet = 2*piG -dPhiJet;
    dPhiJet = (180*dPhiJet)/piG;
    fdPhiLeadingPairMC->Fill(dPhiJet);
    fptRatioLeadingPairMC->Fill(secondJet->Pt()/leadingJet->Pt());
  }
}

int saModuleDimuonDYDarshana::GetNumberofTracklets(DSTReader *fvtx_trk_map, const DiMuon *dimuon)
{
	if(!fvtx_trk_map){
	cout<<"EXCEPTION: "<<PHWHERE<<endl;
	return NULL;
	}
	
	int ntrklets = 0;
	TClonesArray *array = fvtx_trk_map->get_FvtxCompactTrk();
	for (int i = 0; i < array->GetSize(); i++) {
	TFvtxCompactTrk *tracklet = dynamic_cast<TFvtxCompactTrk*> (array->At(i));
	
	if(!tracklet){
	//cout<<"No tracklet"<<__LINE__<<"size: "<<array->GetSize()<<endl;
	break;
	}
	
	if(_use_cut_tracklet_chi2 && (tracklet->get_chi2_ndf() > _cut_tracklet_chi2)) continue;
	if(!_use_2_hit_tracklet && tracklet->get_nhits() <= 2) continue;

	SingleMuon *muon0 = singlemuoncontainer->get_SingleMuon(0);
	SingleMuon *muon1 = singlemuoncontainer->get_SingleMuon(1);

	float xx0 = tracklet->get_fvtx_vtx().getX()-(tracklet->get_fvtx_vtx().getZ()- dimuons->get_Evt_fvtxZ())*
	tan(tracklet->get_fvtx_theta())*cos(tracklet->get_fvtx_phi());
	float yy0 = tracklet->get_fvtx_vtx().getY()-(tracklet->get_fvtx_vtx().getZ()- dimuons->get_Evt_fvtxZ())*
	tan(tracklet->get_fvtx_theta())*sin(tracklet->get_fvtx_phi());
	float x0y0=sqrt((xx0 - dimuons->get_Evt_fvtxX())*(xx0 - dimuons->get_Evt_fvtxX()) +
	(yy0 - dimuons->get_Evt_fvtxY())*(yy0 - dimuons->get_Evt_fvtxY()));

	if (fabs(TMath::ATan2(sqrt(muon0->get_px_fvtxmutr()*muon0->get_px_fvtxmutr()+
	muon0->get_py_fvtxmutr()*muon0->get_py_fvtxmutr()),muon0->get_pz_fvtxmutr())-
	tracklet->get_fvtx_theta()) >0.001 && fabs(TMath::ATan2(sqrt(muon1->get_px_fvtxmutr()*
	muon1->get_px_fvtxmutr()+muon1->get_py_fvtxmutr()*muon1->get_py_fvtxmutr()),
	muon1->get_pz_fvtxmutr())-tracklet->get_fvtx_theta()) >0.001 && fabs(TMath::ATan2(muon0->get_py_fvtxmutr(),
	muon0->get_px_fvtxmutr())-tracklet->get_fvtx_phi())>0.001 && fabs(TMath::ATan2(muon1->get_py_fvtxmutr(),
	muon1->get_px_fvtxmutr())-tracklet->get_fvtx_phi())>0.001 && tracklet->get_fvtx_theta()+0!=0 && 
	x0y0 < 1.5){
	ntrklets++;
	}

	}

	return ntrklets;
        
}

int main(int argc, char* argv[])
{ 
  //Upload the file with the data
  TFile* file = TFile::Open("/Users/Fer/Documents/traajo/samples/NeroNtuples_9.root"); // TFile::Open() instead of a constructor since it works over xrootd etc.
  //Upload the tree with the event data
  TTree *tree=(TTree*)file->Get("nero/events");

  //Create the vector to store all the particle identifiers
  std::vector<Int_t> * lepPdgId;

  //Create a variable to store all the lepton event data
  TClonesArray *leptondata = new TClonesArray("leptondata");

  //Specify where all the lepton event data will be stores
  tree->SetBranchAddress("lepP4", &leptondata);
  //Specify where all the lepton identifiers will be stored
  tree->SetBranchAddress("lepPdgId", &lepPdgId);

  //Get how many events we have to loop through
  int nentries = tree->GetEntries();

  //Loop through all the events
  for(int ientry = 0; ientry < nentries; ientry++) 
  {
    //Reset the lepton data 
    leptondata->Clear();
    //This line stores the proper data both in "leptondata" and in "lepPdgId"
    tree->GetEntry(ientry);
    
    //Only if "leptondata" is not empty continue, this is to avoid segmentation errors
    if(leptondata->GetSize() == 0) continue;

    //Loop through all the entries in the current event 
    for(int j=0; j<leptondata->GetEntriesFast()-1; j++) 
    {
        //Only if the identifier of the particle is + or - 11 (electron or antielectron) store the data in electrondata
        if(abs(lepPdgId->at(j))==11) continue;
        //Store all the data of the electron in this variable
        TLorentzVector *electrondata = (TLorentzVector *)leptondata->At(j);
        //Get some specific property such as momentum, position or energy
        cout << electrondata->E() << endl;
    }
  }
  return 0; 
}

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

  //Variables for the for loop

  //Get how many events we have to loop through
  int nentries = tree->GetEntries();

  //Create a variable to store the mass values
  Double_t mass;

  //Loop through all the events
  for(int ientry = 0; ientry < nentries; ientry++) 
  {

    //Variable of the whole data
    TLorentzVector addable_lorentz_wholedata;

    //Reset the lepton data 
    leptondata->Clear();

    //Reset the met data 
    metdata->Clear();

    //Reset the jet data 
    jetdata->Clear();

    //This line stores the proper data in the variables qe specified
    tree->GetEntry(ientry);

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

//Implementation of lepton criteria 

    //If "leptondata" is empty it skips and the for loop continues, this is to avoid segmentation errors
    if(leptondata->GetSize() != 1) continue;

      //Only if the identifier of the particle is ±11 or ± 13 (electrons muons or their anti particles) we use the data
      if((abs(lepPdgId->at(0))!=11)&&(abs(lepPdgId->at(0))!=13)) continue; 

      //Store all the data of the electron in this lorentz vector
      TLorentzVector * lorentz_leptondata = (TLorentzVector *)leptondata->At(0);

      //We create another lorentz vector that isn't a pointer for the same reasons that before
      TLorentzVector addable_lorentz_leptondata = *lorentz_leptondata;

      //Get the Eta value of that Lorentz vector
      eta=addable_lorentz_leptondata.Eta();

      //Get the transverse mass of that lorentz vector
      mass=addable_lorentz_leptondata.Mt();

      //Implementing slection criteria
      if (mass<40 || abs(eta)>2.5) continue;

      //Fill the histogram with the current data
      lepmass->Fill(mass);

/////////////////////////////////////////////////////
    
//Implementation of met criteria 

    //Create a lorentz vector with the matdata of the current entry
    TLorentzVector * lorentz_metdata = (TLorentzVector *) metdata->At(0);

    //We cannot use math with pointers for some reason, so we create a lorentz vectors that isn't a pointers
    TLorentzVector addable_lorentz_metdata = *lorentz_metdata;

//.........这里部分代码省略.........
      // show current run/iCat processed
      if(iEvent==0) {
	nCurrentRun = nRun ;
	outlog << "nRun=" << nRun << endl;
      }
      else if(nRun!=nCurrentRun) {
	nCurrentRun=nRun ;
	outlog << "nRun=" << nRun << endl;
      }

//      // run selection (using both json files)
//      //int iJson = detJson(nRun);
//      //int iJson = 4;
//      if(debug) cout << "iJson = " << iJson << endl;
//      outlog << "iJson = " << iJson << endl;
//      if( iJson>-1 && iJson<9) {
//	isGoodRun = AcceptEventByRunAndLumiSection(nRun, nLumi, jsonMap[iJson]);
//	if(!isGoodRun) {
//	  outlog << "failed JSON" << endl;
//	  continue;
//	}
//      }
//      else {
//	outlog << "no proper JSON file" << endl;
//	//continue;
//      }

      // at least 2 electrons
      if(ele_N<2) continue;
      else outlog << "ele_N=" << ele_N << endl;

      // LOOP OVER ELECTRONS //
      if(debug) cout << "<-- ele_N=" << ele_N << endl
		     << "--- electrons.size=" << electrons->GetSize() << endl;
      for( int iEle1=0 ; iEle1<ele_N ; iEle1++ ) {
	if(debug) cout << "--- get ele #" << iEle1 << endl;
	cand[0] = (TLorentzVector*) (electrons->At (iEle1)) ;
	if(debug) cout << "--- got it" << endl;

	// severity selection
	if( ele_severityLevelSeed[iEle1] >= 3 ) continue;
	
	// check whether electrons of the pair pass HLT_Ele27 Id/Iso cuts
	if(debug) cout << "--- checks VBTF cuts" << endl;
	cut_HLT_Ele27[0] = VBTFcuts( "HLT_Ele27", RunPhase,
				     cand[0]->Pt(), cand[0]->Et(), ele_sclEta[iEle1], cand[0]->Eta(), ele_tkSumPt_dr03[iEle1], ele_ecalRecHitSumEt_dr03[iEle1], 
				     ele_hcalDepth1TowerSumEt_dr03[iEle1], ele_hcalDepth2TowerSumEt_dr03[iEle1], ele_expected_inner_hits[iEle1],
				     ele_deltaphiin[iEle1], ele_deltaetain[iEle1], ele_he[iEle1], ele_sigmaietaieta[iEle1],
				     ele_conv_dist[iEle1], ele_conv_dcot[iEle1], ele_fbrem[iEle1], ele_isConversion[iEle1] ) ;

	// check if ele is a good tag candidate : pass VBTF 95 and has pT>5 GeV
	cutEle[0] = 0;
	cutEle[0] = whichCuts( RunPhase, cand[0]->Pt(), cand[0]->Et(), ele_sclEta[iEle1], cand[0]->Eta(), ele_tkSumPt_dr03[iEle1], ele_ecalRecHitSumEt_dr03[iEle1], 
			       ele_hcalDepth1TowerSumEt_dr03[iEle1], ele_hcalDepth2TowerSumEt_dr03[iEle1], ele_expected_inner_hits[iEle1],
			       ele_deltaphiin[iEle1], ele_deltaetain[iEle1], ele_he[iEle1], ele_sigmaietaieta[iEle1],
			       ele_conv_dist[iEle1], ele_conv_dcot[iEle1], ele_fbrem[iEle1], ele_isConversion[iEle1] ) ;
	fidu[0] = 0;
	if ( fabs(ele_sclEta[iEle1]) < 2.5 && ( fabs(ele_sclEta[iEle1]) > 1.566 || fabs(ele_sclEta[iEle1])<1.4442 ) ) 
	  fidu[0] = 1 ;

	if( cutEle[0]>0 && cand[0]->Et()>=5. ) {
	  if(debug) cout << "--- ele #" << iEle1 << " is a good tag candidate" << endl;
	  
	  // loop to find probe candidates
	  for( int iEle2=0 ; iEle2<ele_N ; iEle2++ ) {
	    if(debug) cout << "----- looks Ele #" << iEle2 << endl;

    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();

    }

///
///________________________________________________________________________________
///
Bool_t
UEJetAreaFinder::find( TClonesArray& Input, vector<UEJetWithArea>& _jets )
{
  /// return if no four-vectors are provided
  if ( Input.GetSize() == 0 ) return kFALSE;

  /// prepare input
  std::vector<fastjet::PseudoJet> fjInputs;
  fjInputs.reserve ( Input.GetSize() );

  int iJet( 0 );
  for( int i(0); i < Input.GetSize(); ++i )
    {
      TLorentzVector *v = (TLorentzVector*)Input.At(i);

      if ( TMath::Abs(v->Eta()) > etaRegionInput ) continue;
      if ( v->Pt()              < ptThreshold    ) continue;

      fjInputs.push_back (fastjet::PseudoJet (v->Px(), v->Py(), v->Pz(), v->E()) );
      fjInputs.back().set_user_index(iJet);
      ++iJet;
    }

  /// return if no four-vectors in visible phase space
  if ( fjInputs.size() == 0 ) return kFALSE;
  
  /// print out info on current jet algorithm
  //   cout << endl;
  //   cout << mJetDefinition->description() << endl;
  //   cout << theAreaDefinition->description() << endl;

  /// return if active area is not chosen to be calculated
  if ( ! theAreaDefinition ) return kFALSE;

  //  cout << "fastjet::ClusterSequenceActiveArea* clusterSequence" << endl;
  
  fastjet::ClusterSequenceArea* clusterSequence
    = new fastjet::ClusterSequenceArea (fjInputs, *mJetDefinition, *theAreaDefinition );

  //  cout << "retrieve jets for selected mode" << endl;

  /// retrieve jets for selected mode
  double mJetPtMin( 1. );
  std::vector<fastjet::PseudoJet> jets( clusterSequence->inclusive_jets (mJetPtMin) );
  unsigned int nJets( jets.size() );

  if ( nJets == 0 ) 
    {
      delete clusterSequence;
      return kFALSE;
    }
  //Double_t ptByArea[ nJets ];

  //   int columnwidth( 10 );
  //cout << "found " << jets.size() << " jets" << endl;
  //   cout.width( 5 );
  //   cout << "jet";
  //   cout.width( columnwidth );
  //   cout << "eta";
  //   cout.width( columnwidth );
  //   cout << "phi";
  //   cout.width( columnwidth );
  //   cout << "pT";
  //   cout.width( columnwidth );
  //   cout << "jetArea";
  //   cout.width( 15 );
  //   cout << "pT / jetArea";
  //   cout << endl;
  
  _jets.reserve( nJets );

  vector< fastjet::PseudoJet > sorted_jets ( sorted_by_pt( jets ));
  for ( int i(0); i<nJets; ++i )
    {
      //ptByArea[i] = jets[i].perp()/clusterSequence->area(jets[i]);

      //       cout.width( 5 );
      //       cout << i;
      //       cout.width( columnwidth );
      //       cout << jets[i].eta();
      //       cout.width( columnwidth );
      //       cout << jets[i].phi();
      //       cout.width( columnwidth );
      //       cout << jets[i].perp();
      //       cout.width( columnwidth );
      //       cout << clusterSequence->area(jets[i]);
      //       cout.width( 15 );
      //       cout << ptByArea[i];
      //       cout << endl;

      /// save
      ///
      /// TLorentzVector
      /// area
      /// nconstituents
      
      fastjet::PseudoJet jet( sorted_jets[i] );
//.........这里部分代码省略.........

void 
UEJetAreaHistograms::fill( vector<UEJetWithArea>& theJets, TClonesArray& acceptedTriggers )
{
  //  cout << "UEJetAreaHistograms::fill( vector<UEJetWithArea>& theJets, TClonesArray& acceptedTriggers )" << endl;

  ///
  /// Histo filler for reco-only analysis
  /// HL trigger bits *are* available
  ///
  
  ///
  /// 11 HLT bits :
  /// 4 Min-Bias (Pixel, Hcal, Ecal, general), Zero-Bias, 6 Jet (30, 50, 80, 110, 180, 250)
  ///
  unsigned int iHLTbit(0);
  for ( ; iHLTbit<11; ++iHLTbit )
    {
      ///
      /// ask if trigger was accepted
      ///
      bool hltAccept( false );
      unsigned int nAcceptedTriggers( acceptedTriggers.GetSize() );
      for ( unsigned int itrig(0); itrig<nAcceptedTriggers; ++itrig )
        {
	  std::string filterName( acceptedTriggers.At(itrig)->GetName() );
          if ( filterName == HLTBitNames[iHLTbit] ) hltAccept = true;
        }
      if ( ! hltAccept ) continue;
      
      Double_t ptarray  [ theJets.size() ];
      Double_t areaarray[ theJets.size() ];
      Double_t ptByArea [ theJets.size() ];
      
      for ( unsigned int ijet(0); ijet < theJets.size(); ++ijet )
	{
	  double pTLeadingJet       ( theJets[0].GetMomentum()->Pt() );
	  
	  double pTJet  ( theJets[ijet].GetMomentum()->Pt()  );
	  double etaJet ( theJets[ijet].GetMomentum()->Eta() );
	  double areaJet( theJets[ijet].GetArea()            );
	  
	  ptarray  [ijet] = pTJet;
	  areaarray[ijet] = areaJet;
	  ptByArea [ijet] = pTJet / areaJet;
	  
	  unsigned int nConstituents( theJets[ijet].GetNConstituents() );
	  
	  h_pTAllJets           [iHLTbit]->Fill( pTJet          );
	  h_etaAllJets          [iHLTbit]->Fill( etaJet         );
	  h_areaAllJets         [iHLTbit]->Fill( areaJet        );
	  h_ptByAreaAllJets     [iHLTbit]->Fill( ptByArea[ijet] );
	  h_nConstituentsAllJets[iHLTbit]->Fill( nConstituents  );

	  h2d_pTAllJets_vs_pTjet           [iHLTbit]->Fill( pTLeadingJet, pTJet          );
	  h2d_areaAllJets_vs_pTjet         [iHLTbit]->Fill( pTLeadingJet, areaJet        );
	  h2d_ptByAreaAllJets_vs_pTjet     [iHLTbit]->Fill( pTLeadingJet, ptByArea[ijet] );
	  h2d_nConstituentsAllJets_vs_pTjet[iHLTbit]->Fill( pTLeadingJet, nConstituents  );

	  ///
	  /// histograms for leading jet
	  ///
	  if ( ijet == 0 )
	    {
	      h_pTJet           [iHLTbit]->Fill( pTJet           );
	      h_etaJet          [iHLTbit]->Fill( etaJet          );
	      h_areaJet         [iHLTbit]->Fill( areaJet         );
	      h_ptByAreaJet     [iHLTbit]->Fill( ptByArea[ijet]  );
	      h_nConstituentsJet[iHLTbit]->Fill( nConstituents   );

	      h2d_areaJet_vs_pTjet         [iHLTbit]->Fill( pTJet, areaJet         );
	      h2d_ptByAreaJet_vs_pTjet     [iHLTbit]->Fill( pTJet, ptByArea[ijet]  );
	      h2d_nConstituentsJet_vs_pTjet[iHLTbit]->Fill( pTJet, nConstituents   );
	    }
	}

      double medianPt       ( TMath::Median( theJets.size(), areaarray ) );
      double medianArea     ( TMath::Median( theJets.size(), ptarray   ) );
      double medianPtPerArea( TMath::Median( theJets.size(), ptByArea  ) );
      
      h_medianPt      [iHLTbit]->Fill( medianArea      );
      h_medianArea    [iHLTbit]->Fill( medianPt        );
      h_medianPtByArea[iHLTbit]->Fill( medianPtPerArea );
      
      h2d_medianPt_vs_pTjet      [iHLTbit]->Fill( ptarray[0], medianPt        );
      h2d_medianArea_vs_pTjet    [iHLTbit]->Fill( ptarray[0], medianArea      );
      h2d_medianPtByArea_vs_pTjet[iHLTbit]->Fill( ptarray[0], medianPtPerArea );
    }
}

void UEAnalysisGAM::gammaAnalysisMC(Float_t weight,Float_t etaRegion,Float_t ptThreshold, TClonesArray& MCGamma, TClonesArray& ChargedJet)
{
  
  vector<TLorentzVector*> JetMC;
  vector<TLorentzVector*> GamMC;
  GamMC.clear();
  JetMC.clear();

  for(int j=0;j<MCGamma.GetSize();++j){
        TLorentzVector *g = (TLorentzVector*)MCGamma.At(j);
        if(fabs(g->Eta())<etaRegion){
	  GamMC.push_back(g);
	  if(GamMC.size()==1) JetMC.push_back(g);
	}
   }

  if(JetMC.size() != 0){
    for(int j=0;j<ChargedJet.GetSize();++j){
      TLorentzVector *w = (TLorentzVector*)ChargedJet.At(j);
      if(fabs(w->Eta())<etaRegion){
        JetMC.push_back(w);
      }
    }
    if(JetMC.size()>=2){
      float dPhiJet1 = fabs(JetMC[0]->Phi()-JetMC[1]->Phi());    
      if(dPhiJet1> piG) dPhiJet1 = 2*piG -dPhiJet1;
      dPhiJet1 = (180*dPhiJet1)/piG;
      fdPhiGamma1JetMC->Fill(dPhiJet1);
    }
    if(JetMC.size()>=3){
      float dPhiJet2 = fabs(JetMC[0]->Phi()-JetMC[2]->Phi());    
      if(dPhiJet2> piG) dPhiJet2 = 2*piG -dPhiJet2;
      dPhiJet2 = (180*dPhiJet2)/piG;
      fdPhiGamma2JetMC->Fill(dPhiJet2);
    }
    if(JetMC.size()>=4){
      float dPhiJet3 = fabs(JetMC[0]->Phi()-JetMC[3]->Phi());    
      if(dPhiJet3> piG) dPhiJet3 = 2*piG -dPhiJet3;
      dPhiJet3 = (180*dPhiJet3)/piG;
      fdPhiGamma3JetMC->Fill(dPhiJet3);
    }
 

    vector<AssociatedObject> assoJetMC;
    assoJetMC.clear();

    while(JetMC.size()>1){
      int oldSize = JetMC.size();
      vector<TLorentzVector*>::iterator itH = JetMC.begin();
      if((*itH)->Pt()>=ptThreshold){
	for(vector<TLorentzVector*>::iterator it=JetMC.begin();it!=JetMC.end();it++){
	  float azimuthDistanceJet = fabs( (*itH)->Phi() - (*it)->Phi() );
	  if((*it)->Pt()/(*itH)->Pt()>=0.3){
	    if( (piG - rangePhi) <  azimuthDistanceJet && azimuthDistanceJet < (piG + rangePhi)){
	      AssociatedObject tmpPair((*itH),(*it));
	      assoJetMC.push_back(tmpPair);
	      JetMC.erase(it);
	      int newSize = oldSize -1;
	      oldSize = newSize;
	      JetMC.resize(newSize);
	      break;
	    }
	  }
	}
      }
      JetMC.erase(itH);
      int newSize = oldSize -1;
      JetMC.resize(newSize);
    }
  
    if(assoJetMC.size()){
      fNumbMPIMC->Fill(assoJetMC.size());
      vector<AssociatedObject>::iterator at= assoJetMC.begin();
    
      const TLorentzVector* leadingJet((*at).first);
      const TLorentzVector* secondJet((*at).second);

      pPtRatio_vs_PtJleadMC->Fill(leadingJet->Pt(),(secondJet->Pt()/leadingJet->Pt()));
      pPtRatio_vs_EtaJleadMC->Fill(fabs(leadingJet->Eta()),(secondJet->Pt()/leadingJet->Pt()));
      pPtRatio_vs_PhiJleadMC->Fill(leadingJet->Phi(),(secondJet->Pt()/leadingJet->Pt()));
    
      fdEtaLeadingPairMC->Fill(leadingJet->Eta()-secondJet->Eta());
      float dPhiJet = fabs(leadingJet->Phi()-secondJet->Phi());
      if(dPhiJet> piG) dPhiJet = 2*piG -dPhiJet;
      dPhiJet = (180*dPhiJet)/piG;
      fdPhiLeadingPairMC->Fill(dPhiJet);
      fptRatioLeadingPairMC->Fill(secondJet->Pt()/leadingJet->Pt());
    }
    
    fPhiLeadingGammaMC->Fill(GamMC[0]->Phi());
    fPtLeadingGammaMC->Fill(GamMC[0]->Pt());
    fEtaLeadingGammaMC->Fill(GamMC[0]->Eta());
  }

}