//Declaration
Float_t GetR(Float_t phi1, Float_t phi2, Float_t eta1, Float_t eta2);

void pythia_multiT_test(const Char_t *out_name="pythia_dphi_histsNew.root",
		    const Int_t n_evt=100) {
  gSystem->Load("libPhysics");
  gSystem->Load("libEG");
  gSystem->Load("$OPTSTAR/alt/lib/libPythia6");
  gSystem->Load("libEGPythia6");

  const Float_t pt_trig_min=5; 
  const Float_t pt_assoc_min=2;
  const Float_t mJetR=0.3;
  const Int_t n_trig = 0;
  const Int_t n_assoc = 0;
  TDatabasePDG *pdg = TDatabasePDG::Instance(); // TDataBasePDG contains info on all particle properties

  TPythia6 *pythia=TPythia6::Instance();

  pythia->SetMRPY(1,clock());  // set random seed
  pythia->SetMRPY(2,0);

  //pythia->SetCKIN(3,2); //99.9);   // minimum hard pt
  //pythia->SetCKIN(4,100.1);  // maximum hard pt

  pythia->Initialize("CMS","p","p",200);

  TFile fout(out_name,"RECREATE");
   
  TH1D *dphi_trig_assoc = new TH1D("dphi_trig_assoc","dphi_trig_assoc",24,-0.5*TMath::Pi(),1.5*TMath::Pi());
  TH1D *dphi_trig_assoc3 = new TH1D("dphi_trig_assoc3","dphi_trig_assoc",24,-0.5*TMath::Pi(),1.5*TMath::Pi());
  TH1D *deta_assoc = new TH1D("deta_assoc","deta_assoc;#eta;dN/d#eta",30,-1.5,1.5);
  TH1D *dphi_assoc = new TH1D("dphi_assoc","dphi_assoc;#phi;dN/d#phi",25,-2*TMath::Pi(),2*TMath::Pi());
  TH1D *dpt_jet = new TH1D("dpt_jet","dpt_jet;p_{T,jet};dN/dp_{T,jet}",24,4.0,30.0);
  TH1D *mPtSeeds = new TH1D("PtSeeds","dN/dPtJet;p_{T,jet}(Gev)", 25,0.0,25.0);//Single Particle Triggers
  TH1D *Triggers = new TH1D("dNdPtJet","dN/dPtJet;p_{T,jet}(Gev)", 25,0.0,25.0);//All Cluster Triggers  
  TH1D *cluster_tracks = new TH1D("cluster_tracks","counts;tracks", 20,0.0,10.0);//Tracks in each cluster
  TH3D *dPhiJetPtAssocJetPt = new TH3D("dPhiJetPtAssocJetPt","dPhiJetPtAssoc;p_{T,assoc} (GeV);#Delta#phi;p_{T,jet} ",25,0.0,25.0,24,-0.5*TMath::Pi(),1.5*TMath::Pi(),24,6.0,30.0);
  TH3D *mdPhiSeedAssocSingleJetPt = new TH3D("dPhiSeedAssocSingleJetPt","dPhiJetPtAssoc;p_{T,assoc} (GeV);#Delta#phi;p_{T,jet} ",25,0.0,25.0,24,-0.5*TMath::Pi(),1.5*TMath::Pi(),26,4.0,30.0);
  TH3D *mdPhiSeedAssocSingleClusterJetPt = new TH3D("dPhiSeedAssocSingleClusterJetPt","dPhiJetPtAssocSingleCluster;p_{T,assoc} (GeV);#Delta#phi;p_{T,jet} ",25,0.0,25.0,24,-0.5*TMath::Pi(),1.5*TMath::Pi(),26,4.0,30.0);//Clusters with only 1 particle  JetPt, associated Pt, delta phi     
  TH3D * mdPhiSeedAssocMultipleClusterJetPt = new TH3D("dPhiSeedAssocMultipleClusterJetPt","dPhiJetPtAssocMultipleCluster;p_{T,assoc} (GeV);#Delta#phi;p_{T,jet} ",25,0.0,25.0,24,-0.5*TMath::Pi(),1.5*TMath::Pi(),26,4.0,30.0);//Clusters with more than 1 particle JetPt, associated Pt, delta phi   
  TH1D *mPtMultiSeeds = new TH1D("PtMultiSeeds","dN/dPtJet;p_{T,jet}(Gev)", 25,0.0,25.0);//Multiple Particle Cluster Triggers 
  TH1D *mPtSingleSeeds = new TH1D("PtSingleSeeds","dN/dPtJet;p_{T,jet}(Gev)", 25,0.0,25.0);//Single Particle Cluster Triggers 
   
  
  for (Int_t i_evt=0; i_evt<n_evt; i_evt++) {
    if (i_evt % 1000 == 0)
      cout << "Event " << i_evt << endl;
    pythia->GenerateEvent();
    /*
      TObjArray *plist=pythia->GetListOfParticles();
      for (Int_t i_part=2; i_part<8; i_part++) {
      part=(TMCParticle*) plist->At(i_part);
      cout << i_part << " " << part->GetKS() << " " << part->GetKF() 
      << " " << part->GetParent() << " ";
      Float_t pt=sqrt(part->GetPx()*part->GetPx()+part->GetPy()*part->GetPy());
      Float_t ptot=sqrt(pt*pt+part->GetPz()*part->GetPz());

      cout << ptot << " " << pt << endl; 
      }
    */
    
    TObjArray *plist=pythia->GetListOfParticles();
    TIter part_it(plist);
    while ((part=(TMCParticle*)part_it())) {
      // Check trigger particle
      if (part->GetKS()<10 &&  // select only final state particles
          sqrt(part->GetPx()*part->GetPx()+part->GetPy()*part->GetPy()) > pt_trig_min &&
          pdg->GetParticle(part->GetKF())->Charge() != 0) {  // select only charged particles
	Int_t n_trk = 1;
	Float_t phi_trig = TMath::ATan2(part->GetPy(), -part->GetPx());
	Float_t pt_jet = sqrt(part->GetPx()*part->GetPx()+part->GetPy()*part->GetPy());
	Float_t pt_seed = sqrt(part->GetPx()*part->GetPx()+part->GetPy()*part->GetPy());
        Float_t theta_trig = TMath::ATan2(part->GetPz(), -pt_seed);
	Float_t eta_trig = - TMath::Log(TMath::Tan(0.5*theta_trig));
	    if (fabs(eta_trig) > 1.0)continue;
	n_trig++;
    	    
	TIter part_it2(plist); 
   	// look for associated particles
	while ((part2=(TMCParticle*)part_it2())) {
	  if (part != part2 &&
	      part2->GetKS()<10 &&  // select only final state particles
	      sqrt(part2->GetPx()*part2->GetPx()+part2->GetPy()*part2->GetPy()) > pt_assoc_min &&
	      pdg->GetParticle(part2->GetKF())->Charge() != 0) {  // select only charged particles
	    n_assoc++;
	    Float_t phi_assoc = TMath::ATan2(part2->GetPy(), -part2->GetPx());
	    dphi_assoc->Fill(phi_assoc);		
	    Float_t dphi = phi_trig - phi_assoc;
	    if (dphi < -0.5*TMath::Pi())
	      dphi += 2*TMath::Pi();
	    if (dphi > 1.5*TMath::Pi())
	      dphi -= 2*TMath::Pi();
	    Float_t pt_assoc = sqrt(part2->GetPx()*part2->GetPx()+part2->GetPy()*part2->GetPy());
	    Float_t theta_assoc = TMath::ATan2(part2->GetPz(), -pt_assoc);
	    Float_t eta_assoc = - TMath::Log(TMath::Tan(0.5*theta_assoc)); 
	    if (fabs(eta_assoc) > 1.0)continue;
	    deta_assoc->Fill(eta_assoc);		
	    dphi_trig_assoc->Fill(dphi);
	    //mdPhiSeedAssocSingleJetPt->Fill(sqrt(part2->GetPx()*part2->GetPx()+part2->GetPy()*part2->GetPy()),dphi,pt_jet);
	    mdPhiSeedAssocSingleJetPt->Fill(sqrt(part2->GetPx()*part2->GetPx()+part2->GetPy()*part2->GetPy()),dphi,pt_seed);
	    //mPtSeeds->Fill(pt_jet);
	    mPtSeeds->Fill(pt_seed);
	    if (pt_assoc > pt_seed)continue;
	    if (GetR(phi_trig,phi_assoc,eta_trig,eta_assoc) < mJetR) {
	    pt_jet += pt_assoc;
	    n_trk++;
	    }
	  }
	}

	TIter part_it3(plist);
        while ((part3=(TMCParticle*)part_it3())) {
	  if (part != part3 && 
	      part3->GetKS()<10 &&  // select only final state particles
	      sqrt(part3->GetPx()*part3->GetPx()+part3->GetPy()*part3->GetPy()) > pt_assoc_min &&
	      pdg->GetParticle(part3->GetKF())->Charge() != 0) {  // select only charged particles
	    Float_t phi_assoc3 = TMath::ATan2(part3->GetPy(), -part3->GetPx());
	    Float_t dphi3 = phi_trig - phi_assoc3;
	    if (dphi3 < -0.5*TMath::Pi())
		dphi3 += 2*TMath::Pi();
	    if (dphi3 > 1.5*TMath::Pi())
		dphi3 -= 2*TMath::Pi();
	    Float_t pt_assoc3 = sqrt(part3->GetPx()*part3->GetPx()+part3->GetPy()*part3->GetPy());
	    Float_t theta_assoc3 = TMath::ATan2(part3->GetPz(), -pt_assoc3);
	    Float_t eta_assoc3 = - TMath::Log(TMath::Tan(0.5*theta_assoc3)); 
	    if (fabs(eta_assoc3) > 1.0)continue;
	   
	    cluster_tracks->Fill(n_trk);
	    
	    dPhiJetPtAssocJetPt->Fill(sqrt(part3->GetPx()*part3->GetPx()+part3->GetPy()*part3->GetPy()),dphi3,pt_jet);
	    Triggers->Fill(pt_jet);

	    if(n_trk == 1)
	    {
		mdPhiSeedAssocSingleClusterJetPt->Fill(sqrt(part3->GetPx()*part3->GetPx()+part3->GetPy()*part3->GetPy()),dphi3,pt_jet);
		mPtSingleSeeds->Fill(pt_jet);
	    }
	    
	    if(n_trk > 1)
	    {
	   	 mdPhiSeedAssocMultipleClusterJetPt->Fill(sqrt(part3->GetPx()*part3->GetPx()+part3->GetPy()*part3->GetPy()),dphi3,pt_jet);
		 mPtMultiSeeds->Fill(pt_jet);
 	    }

	    dphi_trig_assoc3->Fill(dphi3);
	  
	  }
	}
    }
      
    }//end of seed loop
  }
  
  fout.Write();
    cout << "number of triggers  " << n_trig << endl;
    cout << "number of associated  " << n_assoc << endl;
}

Float_t GetR(Float_t phi1, Float_t phi2, Float_t eta1, Float_t eta2){
 Float_t dphi=phi1-phi2;
   if (dphi < -TMath::Pi())
	    dphi += 2*TMath::Pi();
	      if (dphi > TMath::Pi())
	    dphi -= 2*TMath::Pi();
	      Float_t deta=eta2-eta1;
	        return sqrt(dphi*dphi+deta*deta);
}