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

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 AnalyseEvents(ExRootTreeReader *treeReader, MyPlots *plots)
{
  TClonesArray *branchJet = treeReader->UseBranch("Jet");
  TClonesArray *branchElectron = treeReader->UseBranch("Electron");
  TClonesArray *branchMissingET = treeReader->UseBranch("MissingET");

  Long64_t allEntries = treeReader->GetEntries();

  cout << "** Chain contains " << allEntries << " events" << endl;

  Jet *jet[2];
  MissingET *met;
  Electron *electron;

  Long64_t entry;

  Int_t i;

  // Loop over all events
  for(entry = 0; entry < allEntries; ++entry)
  {
    // Load selected branches with data from specified event
    treeReader->ReadEntry(entry);

    // Analyse two leading jets
    if(branchJet->GetEntriesFast() >= 2)
    {
      jet[0] = (Jet*) branchJet->At(0);
      jet[1] = (Jet*) branchJet->At(1);

      plots->fJetPT[0]->Fill(jet[0]->PT);
      plots->fJetPT[1]->Fill(jet[1]->PT);
    }

    // Analyse missing ET
    if(branchMissingET->GetEntriesFast() > 0)
    {
      met = (MissingET*) branchMissingET->At(0);
      plots->fMissingET->Fill(met->MET);
    }

    // Loop over all electrons in event
    for(i = 0; i < branchElectron->GetEntriesFast(); ++i)
    {
      electron = (Electron*) branchElectron->At(i);
      plots->fElectronPT->Fill(electron->PT);
    }
  }
}

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

Int_t printDigits()
{

    AliRunLoader *rl = AliRunLoader::Open("galice.root");

    AliPHOSLoader *prl = (AliPHOSLoader*)rl->GetDetectorLoader("PHOS");

    prl->LoadDigits("READ");
    //prl->LoadDigits();

    Int_t nDigits = 0;
    Int_t nSimEvents = rl->GetNumberOfEvents();
    for (Int_t ev = 0; ev < nSimEvents; ev++)
    {
        rl->GetEvent(ev);

        Int_t nPhosDigits = prl->Digits()->GetEntries();

        //Int_t nDigsFound = 0;
        std::cout << "Number of digits found: " << nPhosDigits << std::endl;
        TClonesArray *phosDigits = prl->Digits();

        for (Int_t iDig = 0; iDig < nPhosDigits; iDig++)
        {
            //const AliPHOSDigit *digit = prl->Digit(iDig);
            AliPHOSDigit *digit = (AliPHOSDigit*)phosDigits->At(iDig);
            nDigits++;
	    //if(digit->GetTime() > 1.4e-08 && digit->GetTime() < 1.6e-08)
            std::cout <<"#: " << iDig << " ID: " << digit->GetId() << " " << "Energy: " << digit->GetEnergy() << " Time: " << digit->GetTime() << " N_prim: " << digit->GetNprimary() <<  " " << digit->GetTimeR() << std::endl;
        }
    }
  return 0;
}

void decayAndFill(int const kf, TLorentzVector* b, double const weight, TClonesArray& daughters)
{
   pydecay->Decay(kf, b);
   pydecay->ImportParticles(&daughters);

   TLorentzVector p1Mom;
   TLorentzVector p2Mom;
   TVector3 v00;

   int nTrk = daughters.GetEntriesFast();
   for (int iTrk = 0; iTrk < nTrk; ++iTrk)
   {
     TParticle* ptl0 = (TParticle*)daughters.At(iTrk);

     switch (ptl0->GetPdgCode())
     {
       //Will only have pions
       case 211:
         ptl0->Momentum(p1Mom);
         v00.SetXYZ(ptl0->Vx() * 1000., ptl0->Vy() * 1000., ptl0->Vz() * 1000.); // converted to μm
         break;
       case -211:
         ptl0->Momentum(p2Mom);
         break;
       default:
         break;
     }
   }
   daughters.Clear();

   fill(kf, b, weight, p1Mom, p2Mom, v00);
}

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

void LHCOWriter::AnalyseEvent()
{
  Event *element;

  element = static_cast<Event*>(fBranchEvent->At(0));

  fprintf(fOutputFile, "%4d %13lld %8d\n", 0, element->Number, 0);

  ++fIntParam[0];
}

//_____________________________________________________________________________
void THaVDC::CorrectTimeOfFlight(TClonesArray& tracks)
{
  const static Double_t v = 3.0e-8;   // for now, assume that everything travels at c

  // get scintillator planes
  THaScintillator* s1 = static_cast<THaScintillator*>
    ( GetApparatus()->GetDetector("s1") );
  THaScintillator* s2 = static_cast<THaScintillator*>
    ( GetApparatus()->GetDetector("s2") );

  if( (s1 == NULL) || (s2 == NULL) )
    return;

  // adjusts caluculated times so that the time of flight to S1
  // is the same as a track going through the middle of the VDC
  // (i.e. x_det = 0) at a 45 deg angle (theta_t and phi_t = 0)
  // assumes that at least the coarse tracking has been performed

  Int_t n_exist = tracks.GetLast()+1;
  //cerr<<"num tracks: "<<n_exist<<endl;
  for( Int_t t = 0; t < n_exist; t++ ) {
    THaTrack* track = static_cast<THaTrack*>( tracks.At(t) );
    
    // calculate the correction, since it's on a per track basis
    Double_t s1_dist, vdc_dist, dist, tdelta;
    if(!s1->CalcPathLen(track, s1_dist))
      s1_dist = 0.0;
    if(!CalcPathLen(track, vdc_dist))
      vdc_dist = 0.0;

    // since the z=0 of the transport coords is inclined with respect
    // to the VDC plane, the VDC correction depends on the location of
    // the track
    if( track->GetX() < 0 )
      dist = s1_dist + vdc_dist;
    else
      dist = s1_dist - vdc_dist;
    
    tdelta = ( fCentralDist - dist) / v;
    //cout<<"time correction: "<<tdelta<<endl;

    // apply the correction
    Int_t n_clust = track->GetNclusters();
    for( Int_t i = 0; i < n_clust; i++ ) {
      THaVDCUVTrack* the_uvtrack = 
	static_cast<THaVDCUVTrack*>( track->GetCluster(i) );
      if( !the_uvtrack )
	continue;
      
      //FIXME: clusters guaranteed to be nonzero?
      the_uvtrack->GetUCluster()->SetTimeCorrection(tdelta);
      the_uvtrack->GetVCluster()->SetTimeCorrection(tdelta);
    }
  }
}

//______________________________________________________________________________
void alice_esd_read()
{
   // Read tracks and associated clusters from current event.

   AliESDRun    *esdrun = (AliESDRun*)    esd->fESDObjects->FindObject("AliESDRun");
   TClonesArray *tracks = (TClonesArray*) esd->fESDObjects->FindObject("Tracks");

   // This needs further investigation. Clusters not shown.
   // AliESDfriend *frnd   = (AliESDfriend*) esd->fESDObjects->FindObject("AliESDfriend");
   // printf("Friend %p, n_tracks:%d\n", frnd, frnd->fTracks.GetEntries());

   if (track_list == 0) {
      track_list = new TEveTrackList("ESD Tracks"); 
      track_list->SetMainColor(6);
      //track_list->SetLineWidth(2);
      track_list->SetMarkerColor(kYellow);
      track_list->SetMarkerStyle(4);
      track_list->SetMarkerSize(0.5);

      gEve->AddElement(track_list);
   }

   TEveTrackPropagator* trkProp = track_list->GetPropagator();
   trkProp->SetMagField( 0.1 * esdrun->fMagneticField ); // kGaus to Tesla

   gProgress->Reset();
   gProgress->SetMax(tracks->GetEntriesFast());
   for (Int_t n=0; n<tracks->GetEntriesFast(); ++n)
   {
      AliESDtrack* at = (AliESDtrack*) tracks->At(n);

      // If ITS refit failed, take track parameters at inner TPC radius.
      AliExternalTrackParam* tp = at;
      if (! trackIsOn(at, kITSrefit)) {
         tp = at->fIp;
      }

      TEveTrack* track = esd_make_track(trkProp, n, at, tp);
      track->SetAttLineAttMarker(track_list);
      track_list->AddElement(track);

      // This needs further investigation. Clusters not shown.
      // if (frnd)
      // {
      //     AliESDfriendTrack* ft = (AliESDfriendTrack*) frnd->fTracks->At(n);
      //     printf("%d friend = %p\n", ft);
      // }
      gProgress->Increment(1);
   }

   track_list->MakeTracks();
}

//_____________________________________________________________________________
Int_t THaVDC::FindVertices( TClonesArray& tracks )
{
  // Calculate the target location and momentum at the target.
  // Assumes that CoarseTrack() and FineTrack() have both been called.

  Int_t n_exist = tracks.GetLast()+1;
  for( Int_t t = 0; t < n_exist; t++ ) {
    THaTrack* theTrack = static_cast<THaTrack*>( tracks.At(t) );
    CalcTargetCoords(theTrack, kRotatingTransport);
  }

  return 0;
}

void testTDime(Int_t nev = 100) {

  gSystem->Load("libEVGEN");
  gSystem->Load("libTDime");
  gSystem->Load("libdime");

  TDime* dime = new TDime();
  dime->SetEnergyCMS(7000.0);
  dime->SetYRange(-2.0, 2.0);   // Set rapidity range of mesons
  dime->SetMinPt(0.1);          // Minimum pT of mesons
  dime->Initialize();

  // (pi+pi-) histograms
  TH1* hM = new TH1D("hM", "DIME #pi^{+}#pi^{-};M_{#pi^{+}#pi^{-}} #[]{GeV/#it{c}^{2}}", 100,  0.0, 5.0);

  TClonesArray* particles = new TClonesArray("TParticle", 6);
  TParticle* part = NULL;
  TLorentzVector v[2];
  TLorentzVector vSum;

  // Event loop
  for (Int_t i = 0; i < nev; ++i) {

    dime->GenerateEvent();
    Int_t np = dime->ImportParticles(particles, "All");
    printf("\n DIME Event %d: Imported %3d particles \n", i, np);

    Int_t nPrimary = 0;

    // Loop over pion (j = 4,5) tracks
    for (Int_t j = 4; j < 6; ++j) {
      part = (TParticle*) particles->At(j); // Choose the particle
      part->Print();
      part->Momentum(v[nPrimary]);          // Copy content to v
      nPrimary++;
    }
    //particles.Clear();

    // 4-vector sum
    vSum = v[0] + v[1];

    // Fill pi+pi- histograms
    hM->Fill(vSum.M());
  }

  // Save plots as pdf
  hM->Draw();    c1->SaveAs("massTDime.pdf");

}

void LHCOWriter::AnalyseMissingET()
{
  MissingET *element;

  element = static_cast<MissingET*>(fBranchMissingET->At(0));

  Reset();

  fIntParam[1] = 6;

  fDblParam[1] = element->Phi;
  fDblParam[2] = element->MET;

  Write();
}

//_____________________________________________________________________________
Int_t THaReacPointFoil::Process( const THaEvData& )
{
  // Calculate the vertex coordinates.

  if( !IsOK() ) return -1;

  Int_t ntracks = fSpectro->GetNTracks();
  if( ntracks == 0 ) return 0;

  TClonesArray* tracks = fSpectro->GetTracks();
  if( !tracks ) return -2;

  TVector3 beam_org, beam_ray( 0.0, 0.0, 1.0 );
  if( fBeam ) {
    beam_org = fBeam->GetPosition();
    beam_ray = fBeam->GetDirection();
  }
  static const TVector3 yax( 0.0, 1.0, 0.0 );
  static const TVector3 xax( 1.0, 0.0, 0.0 );
  TVector3 org, v; 
  Double_t t;

  for( Int_t i = 0; i<ntracks; i++ ) {
    THaTrack* theTrack = static_cast<THaTrack*>( tracks->At(i) );
    // Ignore junk tracks
    if( !theTrack || !theTrack->HasTarget() ) 
      continue;  
    org.SetX( 0. );
    org.SetZ( 0. ); 
    org.SetY( 0. );
    if( !IntersectPlaneWithRay( xax, yax, org, 
				beam_org, beam_ray, t, v ))
      continue; // Oops, track and beam parallel?
    theTrack->SetVertex(v);

    // FIXME: preliminary
    if( theTrack == fSpectro->GetGoldenTrack() ) {
      fVertex = theTrack->GetVertex();
      fVertexOK = kTRUE;
    }
    // FIXME: calculate vertex coordinate errors here (need beam errors)


  }
  return 0;
}

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