C++ TIterator::Reset方法代码示例

void LHCOWriter::AnalyseJets()
{
  Jet *element;
  Track *track;
  Int_t counter;

  fItJet->Reset();
  while((element = static_cast<Jet*>(fItJet->Next())))
  {
    if(element->TauTag != 0) continue;

    Reset();

    counter = 0;
    fItTrack->Reset();
    while((track = static_cast<Track*>(fItTrack->Next())))
    {
      if(element->P4().DeltaR(track->P4()) < 0.5) ++counter;
    }

    fIntParam[1] = 4;

    fDblParam[0] = element->Eta;
    fDblParam[1] = element->Phi;
    fDblParam[2] = element->PT;
    fDblParam[3] = element->Mass;
    fDblParam[4] = counter;
    fDblParam[5] = element->BTag;
    fDblParam[6] = element->EhadOverEem;

    Write();
  }
}

void LHCOWriter::AnalysePhotons()
{
  Photon *element;

  fItPhoton->Reset();
  while((element = static_cast<Photon*>(fItPhoton->Next())))
  {
    Reset();

    fIntParam[1] = 0;

    fDblParam[0] = element->Eta;
    fDblParam[1] = element->Phi;
    fDblParam[2] = element->PT;

    fDblParam[6] = element->EhadOverEem;

    Write();
  }
}

void LHCOWriter::AnalyseElectrons()
{
  Electron *element;

  fItElectron->Reset();
  while((element = static_cast<Electron*>(fItElectron->Next())))
  {
    Reset();

    fIntParam[1] = 1;

    fDblParam[0] = element->Eta;
    fDblParam[1] = element->Phi;
    fDblParam[2] = element->PT;

    fDblParam[4] = element->Charge;

    fDblParam[6] = element->EhadOverEem;

    Write();
  }
}

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

  vector<string> foundChannels;
  vector<string> skippedChannels;  

  cout << "Getting constraints" << endl;
  map<string, RooDataSet*> data_map;
  map<string, RooAbsPdf*> pdf_map;
  RooCategory* decCat = new RooCategory("dec_channel","dec_channel");
  // int i = 0;
  TIterator* catItr = cat->typeIterator();
  RooCatType* type;
  RooArgSet allConstraints;
  while ((type = (RooCatType*)catItr->Next())) {
    RooAbsPdf* pdf =  simPdf->getPdf(type->GetName());

    string typeName(type->GetName());
    if (channelNames.size() && channelNames.find(typeName) == channelNames.end())  {
      skippedChannels.push_back(typeName);
      continue;
    }
    cout << "On channel " << type->GetName() << endl;
    foundChannels.push_back(typeName);

    decCat->defineType(type->GetName());
    // pdf->getParameters(*data)->Print("v");

    RooArgSet nuis_tmp1 = nuis;
    RooArgSet nuis_tmp2 = nuis;
    RooArgSet* constraints = pdf->getAllConstraints(*mc_obs, nuis_tmp1, true);
    constraints->Print();
    allConstraints.add(*constraints);
  }

  catItr->Reset();

  while ((type = (RooCatType*)catItr->Next())) {
    RooAbsPdf* pdf =  simPdf->getPdf(type->GetName());

    string typeName(type->GetName());
    cout << "Considering type " << typeName << endl;
    if (channelNames.size() && channelNames.find(typeName) == channelNames.end()) continue;
    cout << "On channel " << type->GetName() << endl;

    RooArgSet nuis_tmp1 = nuis;
    RooArgSet nuis_tmp2 = nuis;
    RooArgSet* constraints = pdf->getAllConstraints(*mc_obs, nuis_tmp1, true);

    cout << "Adding pdf to map: " << typeName << " = " << pdf->GetName() << endl;
    pdf_map[typeName] = pdf;

    RooProdPdf prod("prod","prod",*constraints);

    RooArgSet* params = pdf->getParameters(*data);
    antiNuis.remove(*params);
    antiGlobs.remove(*params);

    allParams.add(*params);
    // cout << type->GetName() << endl;
  }
  // return;

  RooArgSet decNuis;
  TIterator* nuiItr = mc_nuis->createIterator();
  TIterator* parItr = allParams.createIterator();
  RooAbsArg* nui, *par;
  while ((par = (RooAbsArg*)parItr->Next())) {

int main(int argc, char *argv[])
{
  // Declaration of variables
  ExRootConfReader *confReader = 0;
  Delphes *modularDelphes = 0;
  DelphesFactory *factory = 0;
  TObjArray *allParticleOutputArray = 0;
  TObjArray *stableParticleOutputArray = 0;
  TObjArray *partonOutputArray = 0;

  Int_t event;

  TObjArray *inputArray = 0;
  TIterator *inputIterator = 0;
  Candidate *candidate = 0;
  TLorentzVector momentum;

  JetDefinition *definition = 0;
  vector<PseudoJet> inputList, outputList;
  PseudoJet jet;

  gROOT->SetBatch();

  int appargc = 1;
  char appName[] = "ExternalFastJetBasic";
  char *appargv[] = {appName};
  TApplication app(appName, &appargc, appargv);

  if(argc != 2)
  {
    cout << " Usage: " << appName << " config_file" << endl;
    cout << " config_file - configuration file in Tcl format." << endl;
    return 1;
  }

  try
  {
    // Initialization
    confReader = new ExRootConfReader;
    confReader->ReadFile(argv[1]);

    modularDelphes = new Delphes("Delphes");
    modularDelphes->SetConfReader(confReader);

    factory = modularDelphes->GetFactory();

    allParticleOutputArray = modularDelphes->ExportArray("allParticles");
    stableParticleOutputArray = modularDelphes->ExportArray("stableParticles");
    partonOutputArray = modularDelphes->ExportArray("partons");

    modularDelphes->InitTask();

    
    // fastjet definition
    ClusterSequence::print_banner();
    definition = new JetDefinition(antikt_algorithm, 0.5);
    
    // Define your input candidates to fastjet (by default particle-flow objects).
    // If you want pure calorimeter towers change "EFlowMerger/eflow" into "Calorimeter/towers":
     
    inputArray = modularDelphes->ImportArray("EFlowMerger/eflow");
      
    inputIterator = inputArray->MakeIterator();

    // Event loop
    for(event = 0; event < NEVENTS; ++event)
    {
      modularDelphes->Clear();
      
      // convert EVENT input array into Delphes internal format
      ConvertInput(event, factory, allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
      
      // run Delphes reconstruction
      modularDelphes->ProcessTask();

      inputList.clear();
      inputIterator->Reset();
      
      
      // pass delphes candidates to fastjet clustering  
      while((candidate = static_cast<Candidate*>(inputIterator->Next())))
      {
        momentum = candidate->Momentum;
        jet = PseudoJet(momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E());
        inputList.push_back(jet);
      }
     
      // run clustering 
      ClusterSequence sequence(inputList, *definition);
      outputList.clear();
      outputList = sorted_by_pt(sequence.inclusive_jets(0.0));

      // tell the user what was done
      //  - the description of the algorithm used
      //  - show the inclusive jets as
      //      {index, rapidity, phi, pt}
      //----------------------------------------------------------
      cout << "Ran " << definition->description() << endl;

      // label the columns
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    TIterator* itercompPars;//forward declare so it persists outside the if statement
    RooRealVar* varcompPars;
    if(!(vartPCtitle.Contains(":")||vartPCtitle.Contains("@"))){//only for non-sub-shapes
      while(i==0||i==10||i==4||i==1||i==5||(i>=10&&i<=27))i++;//avoid white and blue and black and yellow and horribleness
      RooArgSet* compPars = vartPC->getParameters(data);//set of the parameters of the component the loop is on
      itercompPars = compPars->createIterator();
      varcompPars = (RooRealVar*) itercompPars->Next();
    
      while(varcompPars){//write and print mean, sig, etc. of sub-shapes
        TString vartitle = varcompPars->GetTitle();
        double varval = varcompPars->getVal();
        TString varvalstring = Form("%f",varval);
        double hi = varcompPars->getErrorHi();
        
        TString varerrorstring = "[exact]";
        if(hi!=-1){
          double lo = varcompPars->getErrorLo();
          double varerror = TMath::Max(fabs(lo),hi);
          varerrorstring = Form("%E",varerror);
        }
        
        outputtext = vartitle+" = "+varvalstring+" +/- "+varerrorstring;
        textfile<<outputtext<<endl;
        cout<<outputtext<<endl;
        
        varcompPars = (RooRealVar*) itercompPars->Next(); 
      }
      totalPdf.plotOn(framex,Name(vartPC->GetName()),LineStyle(kDashed),LineColor(i),Components(vartPC->GetName()));
      leg->AddEntry(framex->findObject(vartPC->GetName()),vartPCtitle,"l");
    
      i++;
    }
    vartPC = (RooAddPdf*) itertPC->Next();
    itercompPars->Reset();//make sure it's ready for the next vartPC
  }
  
  // Calculate chi2/ndf
  RooArgSet *floatpar = totalPdf.getParameters(data);
  int floatpars = (floatpar->selectByAttrib("Constant",kFALSE))->getSize();
  Double_t chi2 = framex->chiSquare("curvetot","Hist",floatpars);
  TString chi2string = Form("%f",chi2);
  //create text box to list important parameters on the plot
  // TPaveText* txt = new TPaveText(0.1,0.5,0.7,0.9,"NBNDC");
  // txt->SetTextSize(0.06);
  // txt->SetTextColor(kBlack);
  // txt->SetBorderSize(0);
  // txt->SetFillColor(0);
  // txt->SetFillStyle(0);
  outputtext = "#chi^{2}/N_{DoF} = "+chi2string;
  cout<<outputtext<<endl;
  textfile<<outputtext<<endl;
  // txt->AddText(outputtext);
  
  // Print stuff
  TIterator* iteryields =  yields.createIterator();
  RooRealVar* varyields = (RooRealVar*) iteryields->Next();//only inherits things from TObject unless class specified
  vector<double> Y, E;//holds yields and associated errors
  vector<TString> YS, ES;//holds strings of the corresponding yields
  int j=0;//count vector position
  int jS=0, jL=0;//these hold the position of the S and L results;initialized in case there is no nsigS or nsigL
  while(varyields){//loop over yields
    TString varname = varyields->GetName();
    TString vartitle = varyields->GetTitle();
    double varval = varyields->getVal();
    Y.push_back(varval);
    double lo = varyields->getErrorLo();

int main(int argc, char *argv[])
{
  char appName[] = "hepmc2pileup";
  stringstream message;
  FILE *inputFile = 0;
  DelphesFactory *factory = 0;
  TObjArray *stableParticleOutputArray = 0, *allParticleOutputArray = 0, *partonOutputArray = 0;
  TIterator *itParticle = 0;
  Candidate *candidate = 0;
  DelphesPileUpWriter *writer = 0;
  DelphesHepMCReader *reader = 0;
  Int_t i;
  Long64_t length, eventCounter;

  if(argc < 2)
  {
    cout << " Usage: " << appName << " output_file" << " [input_file(s)]" << endl;
    cout << " output_file - output binary pile-up file," << endl;
    cout << " input_file(s) - input file(s) in HepMC format," << endl;
    cout << " with no input_file, or when input_file is -, read standard input." << endl;
    return 1;
  }

  signal(SIGINT, SignalHandler);

  gROOT->SetBatch();

  int appargc = 1;
  char *appargv[] = {appName};
  TApplication app(appName, &appargc, appargv);

  try
  {
    writer = new DelphesPileUpWriter(argv[1]);

    factory = new DelphesFactory("ObjectFactory");
    allParticleOutputArray = factory->NewPermanentArray();
    stableParticleOutputArray = factory->NewPermanentArray();
    partonOutputArray = factory->NewPermanentArray();

    itParticle = stableParticleOutputArray->MakeIterator();

    reader = new DelphesHepMCReader;

    i = 2;
    do
    {
      if(interrupted) break;

      if(i == argc || strncmp(argv[i], "-", 2) == 0)
      {
        cout << "** Reading standard input" << endl;
        inputFile = stdin;
        length = -1;
      }
      else
      {
        cout << "** Reading " << argv[i] << endl;
        inputFile = fopen(argv[i], "r");

        if(inputFile == NULL)
        {
          message << "can't open " << argv[i];
          throw runtime_error(message.str());
        }

        fseek(inputFile, 0L, SEEK_END);
        length = ftello(inputFile);
        fseek(inputFile, 0L, SEEK_SET);

        if(length <= 0)
        {
          fclose(inputFile);
          ++i;
          continue;
        }
      }

      reader->SetInputFile(inputFile);

      ExRootProgressBar progressBar(length);

      // Loop over all objects
      eventCounter = 0;
      factory->Clear();
      reader->Clear();
      while(reader->ReadBlock(factory, allParticleOutputArray,
        stableParticleOutputArray, partonOutputArray) && !interrupted)
      {
        if(reader->EventReady())
        {
          ++eventCounter;

          itParticle->Reset();
          while((candidate = static_cast<Candidate*>(itParticle->Next())))
          {
            const TLorentzVector &position = candidate->Position;
            const TLorentzVector &momentum = candidate->Momentum;
            writer->WriteParticle(candidate->PID,
              position.X(), position.Y(), position.Z(), position.T(),
//.........这里部分代码省略.........

void LHCOWriter::AnalyseMuons()
{
  Muon *element;
  Track *track;
  Tower *tower;
  Jet *jet;
  Int_t muonCounter, tauCounter, jetCounter, minIndex;
  Float_t sumPT, sumET, ratET, jetDR, minDR;

  muonCounter = 0;
  fItMuon->Reset();
  while((element = static_cast<Muon*>(fItMuon->Next())))
  {
    Reset();

    sumPT = 0.0;
    fItTrack->Reset();
    while((track = static_cast<Track*>(fItTrack->Next())))
    {
      if(element->P4().DeltaR(track->P4()) < 0.5) sumPT += track->PT;
    }

    sumET = 0.0;
    fItTower->Reset();
    while((tower = static_cast<Tower*>(fItTower->Next())))
    {
      if(element->P4().DeltaR(tower->P4()) < 0.5) sumET += tower->ET;
    }

    tauCounter = 0;
    jetCounter = 0;
    minIndex = -1;
    minDR = 1.0E9;
    fItJet->Reset();
    while((jet = static_cast<Jet*>(fItJet->Next())))
    {
      if(jet->TauTag != 0)
      {
        ++tauCounter;
        continue;
      }

      jetDR = element->P4().DeltaR(jet->P4());
      if(jetDR < minDR)
      {
        minIndex = jetCounter;
        minDR = jetDR;
      }
      ++jetCounter;
    }

    fIntParam[1] = 2;

    fDblParam[0] = element->Eta;
    fDblParam[1] = element->Phi;
    fDblParam[2] = element->PT;

    fDblParam[3] = 0.11;

    fDblParam[4] = element->Charge;

    if(minIndex >= 0)
    {
      fDblParam[5] = fIntParam[0] + fBranchMuon->GetEntriesFast() - muonCounter + tauCounter + minIndex;
    }

    ratET = sumET/element->PT;
    fDblParam[6] = Float_t(TMath::Nint(sumPT)) + (ratET < 1.0 ? ratET : 0.99);

    Write();
    ++muonCounter;
  }
}

int main(int argc, char *argv[])
{
  char appName[] = "root2pileup";
  stringstream message;
  TChain *inputChain = 0;
  ExRootTreeReader *treeReader = 0;
  TClonesArray *branchParticle = 0;
  TIterator *itParticle = 0;
  GenParticle *particle = 0;
  DelphesPileUpWriter *writer = 0;
  Long64_t entry, allEntries;
  Int_t i;

  if(argc < 3)
  {
    cout << " Usage: " << appName << " output_file" << " input_file(s)" << endl;
    cout << " output_file - output binary pile-up file," << endl;
    cout << " input_file(s) - input file(s) in ROOT format." << endl;
    return 1;
  }

  signal(SIGINT, SignalHandler);

  gROOT->SetBatch();

  int appargc = 1;
  char *appargv[] = {appName};
  TApplication app(appName, &appargc, appargv);

  try
  {
    inputChain = new TChain("Delphes");
    for(i = 2; i < argc && !interrupted; ++i)
    {
      inputChain->Add(argv[i]);
    }

    treeReader = new ExRootTreeReader(inputChain);
    branchParticle = treeReader->UseBranch("Particle");
    itParticle = branchParticle->MakeIterator();

    writer = new DelphesPileUpWriter(argv[1]);

    allEntries = treeReader->GetEntries();
    cout << "** Input file(s) contain(s) " << allEntries << " events" << endl;

    if(allEntries > 0)
    {
      ExRootProgressBar progressBar(allEntries - 1);
      // Loop over all events in the input file
      for(entry = 0; entry < allEntries && !interrupted; ++entry)
      {
        if(!treeReader->ReadEntry(entry))
        {
          cerr << "** ERROR: cannot read event " << entry << endl;
          break;
        }

        itParticle->Reset();
        while((particle = static_cast<GenParticle*>(itParticle->Next())))
        {
          writer->WriteParticle(particle->PID,
            particle->X, particle->Y, particle->Z, particle->T,
            particle->Px, particle->Py, particle->Pz, particle->E);
        }
        
        writer->WriteEntry();

        progressBar.Update(entry);
      }
      progressBar.Finish();

      writer->WriteIndex();
    }

    cout << "** Exiting..." << endl;

    delete writer;
    delete itParticle;
    delete treeReader;
    delete inputChain;
    return 0;
  }
  catch(runtime_error &e)
  {
    if(writer) delete writer;
    if(itParticle) delete itParticle;
    if(treeReader) delete treeReader;
    if(inputChain) delete inputChain;
    cerr << "** ERROR: " << e.what() << endl;
    return 1;
  }
}

//.........这里部分代码省略.........
      vertexX = fPrimVertReco.getX();
      vertexY = fPrimVertReco.getY();
      vertexZ = fPrimVertReco.getZ();
      
      hvertexXZ->Fill(vertexZ,vertexX);
      hvertexXY->Fill(vertexY,vertexX);
      hvertexX ->Fill(vertexX);
      hvertexY ->Fill(vertexY);
      hvertexZ ->Fill(vertexZ);
      
      //------Alexandr's-vertex-cuts---------
      
      tgtChi2 = (Float_t) fPrimVertReco.getChi2();

      hTargChi2->Fill(tgtChi2);


      if( tgtChi2<2.0 || tgtChi2>40.) continue; //-skip-events-with-badly-reconstructed-target-------
      if( vertexZ<-59.|| vertexZ>0. ) continue; //-skip-events-with-primary-vertex-outside-of-target-

      hTargX_EVT->Fill(vertexX);
      hTargY_EVT->Fill(vertexY);
      hTargZ_EVT->Fill(vertexZ);

      if( (vertexX-tgtXc)*(vertexX-tgtXc)+(vertexY-tgtYc)*(vertexY-tgtYc)>R2 ) continue; //-skip-events-with-primary-vertex-far-from-target---
      hCutSteps->Fill(5.5);

      hTargXYsel->Fill(vertexX,vertexY);
      hTargXZsel->Fill(vertexZ,vertexX);
      //---------end-of-Alexandr's-vertex-cut--------------
      
      //------HWall-iteration------------------------------
      if (iterWallHit){
    iterWallHit->Reset();
    HWallHit* pWallHit;
        
    while ( (pWallHit = (HWallHit*)iterWallHit->Next()) != 0 ){
         
         cellNum    = pWallHit->getCell();
             cellTime   = pWallHit->getTime();
             cellCharge = pWallHit->getCharge();
             cellPhi    = pWallHit->getPhi();
         
         hWallHitTIME->Fill(cellNum-0.1, cellTime  );
         hWallHitCHRG->Fill(cellNum-0.1, cellCharge);
         
         if (cellCharge > d.dEdxCut(cellNum)/10. && cellTime > 16. && cellTime < 44.){
           hWallHitTIMEcut->Fill(cellNum-0.1, cellTime  );
           hWallHitCHRGcut->Fill(cellNum-0.1, cellCharge);
         }//end-cut
         
    }//------end-HWall-iteration-----------------------
      }//endif

    if( evtinfo->getSumParticleCandMult()<2 ) continue;  //getSumParticleCandMult()>1
      if(!evtinfo->isGoodEvent(
                               Particle::kGoodVertexClust|
                               Particle::kGoodVertexCand|
                               Particle::kNoPileUpSTART|
                               Particle::kGoodSTART|      
                               Particle::kGoodTRIGGER|    
                               Particle::kNoVETO|         
                               Particle::kGoodSTARTVETO|  
                               Particle::kGoodSTARTMETA
                              )) continue ; // bitwise add flags as suggested by Manuel 
                              

void DrawChannelCompatibility(double rMin = -5,double rMax=5)
{
  gROOT->ForceStyle();
  TFile *inf = TFile::Open("higgsCombineZ.ChannelCompatibilityCheck.mH120.root");
  RooFitResult *fit_nominal   = (RooFitResult *)inf->Get("fit_nominal");
  RooFitResult *fit_alternate = (RooFitResult *)inf->Get("fit_alternate");
  RooRealVar *rFit = (RooRealVar*)fit_nominal->floatParsFinal().find("r");
  
  TString prefix = TString::Format("_ChannelCompatibilityCheck_%s_","r");

  int nChann = 0;
  TIterator *iter = fit_alternate->floatParsFinal().createIterator();
  for (RooAbsArg *a = (RooAbsArg *) iter->Next(); a != 0; a = (RooAbsArg *) iter->Next()) {
    if (TString(a->GetName()).Index(prefix) == 0) nChann++;
  }
  TH2F *frame = new TH2F("frame",";best fit #sigma/#sigma_{SM};",1,rMin,rMax,nChann,0,nChann);

  iter->Reset(); 
  int iChann = 0; 
  TGraphAsymmErrors *points = new TGraphAsymmErrors(nChann);
  float chi2(0.0);
  for (RooAbsArg *a = (RooAbsArg *) iter->Next(); a != 0; a = (RooAbsArg *) iter->Next()) {
    if (TString(a->GetName()).Index(prefix) == 0) {
      RooRealVar *ri = (RooRealVar *) a;
      TString channel = a->GetName(); 
      channel.ReplaceAll(prefix,"");
      points->SetPoint(iChann,ri->getVal(),iChann+0.5);
      cout<<channel<<" "<<ri->getVal()<<" "<<ri->getAsymErrorLo()<<" +"<<ri->getAsymErrorHi()<<endl;
      chi2 += pow((ri->getVal()-rFit->getVal())/ri->getError(),2);
      points->SetPointError(iChann,-ri->getAsymErrorLo(),ri->getAsymErrorHi(),0,0);
      //points->SetPointError(iChann,ri->getAsymErrorHi(),ri->getAsymErrorHi(),0,0);
      iChann++;
      frame->GetYaxis()->SetBinLabel(iChann, channel);
    }
  }
  cout<<"Combined fit: "<<rFit->getVal()<<" "<<rFit->getAsymErrorLo()<<" +"<<rFit->getAsymErrorHi()<<endl;
  cout<<"chi2 = "<<chi2<<endl;
  points->SetLineColor(kRed);
  points->SetLineWidth(3);
  points->SetMarkerStyle(21);

  TCanvas *can = new TCanvas("ChannelCompatibility_Z","ChannelCompatibility_Z",900,600);
  frame->GetXaxis()->SetNdivisions(505);
  frame->GetXaxis()->SetTitleSize(0.06);
  frame->GetXaxis()->SetTitleOffset(0.9);
  frame->GetXaxis()->SetLabelSize(0.05);
  frame->GetYaxis()->SetLabelSize(0.1);
  frame->Draw(); 
  //gStyle->SetOptStat(0);
  TBox globalFitBand(rFit->getVal()+rFit->getAsymErrorLo(), 0, rFit->getVal()+rFit->getAsymErrorHi(), nChann);
  //TBox globalFitBand(rFit->getVal()-rFit->getAsymErrorHi(), 0, rFit->getVal()+rFit->getAsymErrorHi(), nChann);
  globalFitBand.SetFillStyle(3013);
  globalFitBand.SetFillColor(65);
  globalFitBand.SetLineStyle(0);
  globalFitBand.DrawClone();
  TLine globalFitLine(rFit->getVal(), 0, rFit->getVal(), nChann);
  globalFitLine.SetLineWidth(4);
  globalFitLine.SetLineColor(214);
  globalFitLine.DrawClone();
  points->Draw("P SAME");
  gPad->Update();

  TLine *ln0 = new TLine(1,gPad->GetFrame()->GetY1(),1,gPad->GetFrame()->GetY2());
  ln0->SetLineColor(kBlack);
  ln0->SetLineWidth(1);
  ln0->SetLineStyle(2);
  ln0->Draw("same");
}