[u/mrichter/AliRoot.git] / PHOS / AliPHOSGammaJet.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
/* $Id$ */

//_________________________________________________________________________
// Class for the analysis of gamma-jet correlations 
//   
// 
//*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN) 
//////////////////////////////////////////////////////////////////////////////


// --- ROOT system ---
#include "TRandom.h"

#include "TLorentzVector.h"
#include "TList.h"
#include "TParticle.h"
#include "AliPHOSGammaJet.h" 
#include "AliPHOSGetter.h" 

ClassImp(AliPHOSGammaJet)

//____________________________________________________________________________
AliPHOSGammaJet::AliPHOSGammaJet() {
  // ctor
}

//____________________________________________________________________________
AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) {
  // ctor 
  AliPHOSGetter::Instance(inputfilename) ; 
}

//____________________________________________________________________________
AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) : TTask(gj) {
}

//____________________________________________________________________________
AliPHOSGammaJet::~AliPHOSGammaJet() {
}

//____________________________________________________________________________
void AliPHOSGammaJet::Exec(Option_t * opt) 
{
  // does the job
  
  if (! opt) 
    return ; 

  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
  Int_t nEvent1 = gime->MaxEvent() ;   
  Int_t iEvent = 0 ; 
  for ( iEvent = 0 ; iEvent < nEvent1 ; iEvent++) {
    if (iEvent <= 100 || iEvent%10 == 0)
      Info("Exec", "Event %d", iEvent) ;     
    TParticle * particle = 0 ;
    //-----------------Fill list with particles--------------------
    TLorentzVector pPi0, pGamma1, pGamma2 ;
    Double_t angle = 0 ;
    TList particleList ;
    Int_t n = -1; 
    gime->Event(iEvent, "X") ; 
    Int_t  nparticles = gime->NPrimaries() ; 
    Int_t iParticle=0 ;
    for (iParticle=0 ; iParticle < nparticles ; iParticle++) {
      //Keep original partons
      particle = gime->Primary(iParticle) ; 
      if((particle->GetStatusCode()== 21)){
	particleList.Add(particle);
	n++;
      }
      //Keep Stable particles within eta range
      Float_t etacut = 0. ; 
      if((particle->GetStatusCode() == 1)&&
	 (TMath::Abs(particle->Eta())<etacut)){
	// Keep particles different from Pi0
	if(particle->GetPdgCode() != 111){
	  particleList.Add(particle);
	  n++;
	} 
	//Decay Pi0 and keep it with different status name
	//Keep decay photons
	if(particle->GetPdgCode() == 111) {
	  //cout<<"Pi0 "<<endl;
	  n += 3 ; 
	  particle->Momentum(pPi0);
	  Pi0Decay(particle->GetMass(),pPi0,pGamma1,pGamma2,angle);	
	  TParticle * photon1 = new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(),
					      pGamma1.Pz(),pGamma1.E(),0,0,0,0);
	  TParticle * photon2 = new TParticle(22,1,0,0,0,0,pGamma2.Px(),pGamma2.Py(),
					      pGamma2.Pz(),pGamma2.E(),0,0,0,0);
	  photon1->SetWeight(1);
	  photon2->SetWeight(2);
	  particle->SetStatusCode(2);
	  particleList.Add(particle);
	  particleList.Add(photon1);
	  particleList.Add(photon2);
	  //photon1->Print();
	  //photon2->Print();  
	}//if pi0
      }//final particle etacut
    }//for (iParticle<nParticle)
    TLorentzVector gamma,charge,pi0, gammapair;
    Int_t idg = -1;
    GetGammaJet(particleList,gamma, idg);
    GetLeadingCharge(particleList,charge, idg);
    GetLeadingPi0(particleList,pi0);
    Info("Pi0Decay", "Gamma: %f %d", gamma.Energy(), idg) ;
    Info("Pi0Decay", "Charge: %f", charge.Energy()) ;
    Info("Pi0Decay", "Pi0: %f", pi0.Energy()) ;
    //    GetLeadingGammaPair(particleList, gammapair, idg, 
    //			0.04,0.2, 1.0,0.13,0.14);
    Info("Pi0Decay", "Pair: %f", gammapair.Energy()) ;
  }//loop: events
}    

//____________________________________________________________________________
void AliPHOSGammaJet::Pi0Decay(Double_t mPi0, TLorentzVector &p0, 
			       TLorentzVector &p1, TLorentzVector &p2, Double_t &angle) {
  // Perform isotropic decay pi0 -> 2 photons
  // p0 is pi0 4-momentum (input)
  // p1 and p2 are photon 4-momenta (output)
  //  cout<<"Boost vector"<<endl;
  TVector3 b = p0.BoostVector();
  //cout<<"Parameters"<<endl;
  //Double_t mPi0   = p0.M();
  Double_t phi    = TMath::TwoPi() * gRandom->Rndm();
  Double_t cosThe = 2 * gRandom->Rndm() - 1;
  Double_t cosPhi = TMath::Cos(phi);
  Double_t sinPhi = TMath::Sin(phi);
  Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe);
  Double_t ePi0   = mPi0/2.;
  //cout<<"ePi0 "<<ePi0<<endl;
  //cout<<"Components"<<endl;
  p1.SetPx(+ePi0*cosPhi*sinThe);
  p1.SetPy(+ePi0*sinPhi*sinThe);
  p1.SetPz(+ePi0*cosThe);
  p1.SetE(ePi0);
  //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
  //cout<<"p1 Mass: "<<p1.Px()*p1.Px()+p1.Py()*p1.Py()+p1.Pz()*p1.Pz()-p1.E()*p1.E()<<endl;
  p2.SetPx(-ePi0*cosPhi*sinThe);
  p2.SetPy(-ePi0*sinPhi*sinThe);
  p2.SetPz(-ePi0*cosThe);
  p2.SetE(ePi0);
  //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
  //cout<<"p2 Mass: "<<p2.Px()*p2.Px()+p2.Py()*p2.Py()+p2.Pz()*p2.Pz()-p2.E()*p2.E()<<endl;
  //cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl;
  p1.Boost(b);
  //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
  p2.Boost(b);
  //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
  //cout<<"angle"<<endl;
  angle = p1.Angle(p2.Vect());
  //cout<<angle<<endl;
}
//____________________________________________________________________________
void AliPHOSGammaJet::GetGammaJet(TList &particleList, TLorentzVector &gamma, Int_t & id) 
{
  // Get the lists of jet particles and gamma
  TParticle *particle = 0x0;
  
  Int_t iPrimary=-1, id_motherg = -1;
  
  TIter next(&particleList);
  while ( (particle = (TParticle*)next()) ) {
    iPrimary++;  
    Int_t ksCode = particle->GetStatusCode();
    Int_t iMother= particle->GetMother(0);
    
    if (ksCode == 21 && iMother == -1)
      if(particle->GetPdgCode()==22){
	id_motherg = iPrimary;
	//cout<<"idmother "<<id_motherg<<endl;
      }
    if(ksCode == 1){
      
      if(id_motherg == particle->GetMother(0)){
 	particle->Momentum(gamma);
 	id = iPrimary;
 	break;
      }
    }// kscode == 1
  }// while
}
 
//____________________________________________________________________________
void AliPHOSGammaJet::GetLeadingCharge(TList &particleList, TLorentzVector &charge, Int_t & id) 
{
  // Gets the leading particle from the list of charged particles
  TParticle *particle = 0x0;
  
  Int_t iPrimary=-1;
  Double_t ptmax = 0, pti = 0;
  TIter next(&particleList);
  while ( (particle = (TParticle*)next()) ) {
    iPrimary++;  
    Int_t ksCode = particle->GetStatusCode();
    
    if((ksCode == 1)&&(id != iPrimary)
       &&(particle->GetPDG(0)->Charge()!=0)){
      pti = particle->Pt(); 
      if(pti> ptmax){
	ptmax = pti;
 	particle->Momentum(charge);
      }//ptmax   
    }// kscode == 1
  }// while
}

//____________________________________________________________________________
void AliPHOSGammaJet::GetLeadingPi0(TList &particleList, TLorentzVector &pi0) 
{
  // Gets the leading pi0 from the list of particles
  TParticle *particle = 0x0;
  
  Int_t iPrimary=-1;
  Double_t ptmax = 0, pti = 0;
  TIter next(&particleList);
  while ( (particle = (TParticle*)next()) ) {
    iPrimary++;  
    Int_t ksCode = particle->GetStatusCode();
    
    if((ksCode == 2))
      {
	pti = particle->Pt(); 
	if(pti> ptmax){
	  ptmax = pti;
	  particle->Momentum(pi0);
	}//ptmax   
      }// kscode == 1
  }// while
}

//____________________________________________________________________________
//  void AliPHOSGammaJet::GetLeadingGammaPair(TList &particleList, TLorentzVector &gammapair, Int_t & id, 
//  			 Double_t & thetacut,Double_t & ratiocut1, Double_t & ratiocut2,
//  			 Double_t & invmasscut1,Double_t & invmasscut2) 
//  {
//   TParticle *particle = 0x0;
  
//   Int_t  iPrimary=-1;
//   Double_t El = 0, E12 = 0;
//   TLorentzVector gamma_i,gamma_j;
//   TIter next(&particleList);
//   while ( (particle = (TParticle*)next()) ) {
//     iPrimary++;	  
//     Int_t ksCode = particle->GetStatusCode();
//     Int_t ksPdg = particle->GetPdgCode();
//     if((ksCode == 1) && (iPrimary != id) && (ksPdg == 22)){
//       particle->Momentum(gamma_i);
//       Int_t jPrimary=-1;
//       TIter next2(&particleList);
//       while ( (particle = (TParticle*)next2()) ) {
// 	jPrimary++;
// 	if(jPrimary>iPrimary){
// 	  ksCode = particle->GetStatusCode();
// 	  ksPdg = particle->GetPdgCode();
// 	  if((ksCode == 1) && (iPrimary != id) && (ksPdg == 22)){
// 	    particle->Momentum(gamma_j);
// 	    if(gamma_j.Angle(gamma_i.Vect())<thetacut){
// 	      Float_t invmass = (gamma_i+gamma_j).M();
// 	      h_invmass->Fill(Eg,invmass);
// 	      if((invmass>invmasscut1) && (invmass<invmasscut2)){
// 		E12 =  (gamma_i+gamma_j).Energy(); 
// 		if(E12>El && (E12/Eg>ratiocut1) && (E12/Eg<ratiocut2)){
// 		  //cout<<E12<<" "<<E12/Eg<<endl;
// 		  El = E12;
// 		  id_i = iPrimary;
// 		  id_j = jPrimary;
// 		  gammapair = gamma_i+gamma_j;
// 		}//E12>El && (E12/Eg>0.2 && E12/Eg<1.)
// 	      }//(invmass>0.125) && (invmass<0.145)
// 	    }//gamma_j.Angle(gamma_i.Vect())<0.04
// 	  }//(ksCode == 1)
// 	}
//       }//while
//       //	    cout<<"jPrimary "<<jPrimary<<endl;
//     }// if kscode 1
//   }//while
//   //	cout<<"iPrimary "<<iPrimary<<endl;
//  }
Commit	Line	Data
77414c90	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15	/* $Id$ */
	16
	17	//_________________________________________________________________________
	18	// Class for the analysis of gamma-jet correlations
	19	//
	20	//
	21	//*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN)
	22	//////////////////////////////////////////////////////////////////////////////
	23
	24
	25	// --- ROOT system ---
e957fea8	26	#include "TRandom.h"
77414c90	27
77414c90	28	#include "TLorentzVector.h"
77414c90	29	#include "TList.h"
77414c90	30	#include "TParticle.h"
77414c90	31	#include "AliPHOSGammaJet.h"
	32	#include "AliPHOSGetter.h"
	33
	34	ClassImp(AliPHOSGammaJet)
	35
	36	//____________________________________________________________________________
	37	AliPHOSGammaJet::AliPHOSGammaJet() {
	38	// ctor
	39	}
	40
	41	//____________________________________________________________________________
	42	AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) {
	43	// ctor
	44	AliPHOSGetter::Instance(inputfilename) ;
	45	}
	46
	47	//____________________________________________________________________________
90cceaf6	48	AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) : TTask(gj) {
77414c90	49	}
	50
	51	//____________________________________________________________________________
	52	AliPHOSGammaJet::~AliPHOSGammaJet() {
	53	}
	54
	55	//____________________________________________________________________________
e957fea8	56	void AliPHOSGammaJet::Exec(Option_t * opt)
77414c90	57	{
	58	// does the job
	59
e957fea8	60	if (! opt)
	61	return ;
	62
77414c90	63	AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
	64	Int_t nEvent1 = gime->MaxEvent() ;
	65	Int_t iEvent = 0 ;
	66	for ( iEvent = 0 ; iEvent < nEvent1 ; iEvent++) {
	67	if (iEvent <= 100 \|\| iEvent%10 == 0)
	68	Info("Exec", "Event %d", iEvent) ;
	69	TParticle * particle = 0 ;
	70	//-----------------Fill list with particles--------------------
	71	TLorentzVector pPi0, pGamma1, pGamma2 ;
	72	Double_t angle = 0 ;
	73	TList particleList ;
	74	Int_t n = -1;
	75	gime->Event(iEvent, "X") ;
	76	Int_t nparticles = gime->NPrimaries() ;
	77	Int_t iParticle=0 ;
	78	for (iParticle=0 ; iParticle < nparticles ; iParticle++) {
	79	//Keep original partons
	80	particle = gime->Primary(iParticle) ;
	81	if((particle->GetStatusCode()== 21)){
	82	particleList.Add(particle);
	83	n++;
	84	}
	85	//Keep Stable particles within eta range
	86	Float_t etacut = 0. ;
	87	if((particle->GetStatusCode() == 1)&&
	88	(TMath::Abs(particle->Eta())<etacut)){
	89	// Keep particles different from Pi0
	90	if(particle->GetPdgCode() != 111){
	91	particleList.Add(particle);
	92	n++;
	93	}
	94	//Decay Pi0 and keep it with different status name
	95	//Keep decay photons
	96	if(particle->GetPdgCode() == 111) {
	97	//cout<<"Pi0 "<<endl;
	98	n += 3 ;
	99	particle->Momentum(pPi0);
	100	Pi0Decay(particle->GetMass(),pPi0,pGamma1,pGamma2,angle);
	101	TParticle * photon1 = new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(),
	102	pGamma1.Pz(),pGamma1.E(),0,0,0,0);
	103	TParticle * photon2 = new TParticle(22,1,0,0,0,0,pGamma2.Px(),pGamma2.Py(),
	104	pGamma2.Pz(),pGamma2.E(),0,0,0,0);
	105	photon1->SetWeight(1);
	106	photon2->SetWeight(2);
	107	particle->SetStatusCode(2);
	108	particleList.Add(particle);
	109	particleList.Add(photon1);
	110	particleList.Add(photon2);
	111	//photon1->Print();
	112	//photon2->Print();
	113	}//if pi0
	114	}//final particle etacut
	115	}//for (iParticle<nParticle)
	116	TLorentzVector gamma,charge,pi0, gammapair;
e957fea8	117	Int_t idg = -1;
	118	GetGammaJet(particleList,gamma, idg);
	119	GetLeadingCharge(particleList,charge, idg);
90cceaf6	120	GetLeadingPi0(particleList,pi0);
e957fea8	121	Info("Pi0Decay", "Gamma: %f %d", gamma.Energy(), idg) ;
77414c90	122	Info("Pi0Decay", "Charge: %f", charge.Energy()) ;
77414c90	123	Info("Pi0Decay", "Pi0: %f", pi0.Energy()) ;
e957fea8	124	// GetLeadingGammaPair(particleList, gammapair, idg,
77414c90	125	// 0.04,0.2, 1.0,0.13,0.14);
	126	Info("Pi0Decay", "Pair: %f", gammapair.Energy()) ;
	127	}//loop: events
	128	}
	129
	130	//____________________________________________________________________________
	131	void AliPHOSGammaJet::Pi0Decay(Double_t mPi0, TLorentzVector &p0,
	132	TLorentzVector &p1, TLorentzVector &p2, Double_t &angle) {
	133	// Perform isotropic decay pi0 -> 2 photons
	134	// p0 is pi0 4-momentum (input)
	135	// p1 and p2 are photon 4-momenta (output)
	136	// cout<<"Boost vector"<<endl;
	137	TVector3 b = p0.BoostVector();
	138	//cout<<"Parameters"<<endl;
	139	//Double_t mPi0 = p0.M();
	140	Double_t phi = TMath::TwoPi() * gRandom->Rndm();
	141	Double_t cosThe = 2 * gRandom->Rndm() - 1;
	142	Double_t cosPhi = TMath::Cos(phi);
	143	Double_t sinPhi = TMath::Sin(phi);
	144	Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe);
	145	Double_t ePi0 = mPi0/2.;
	146	//cout<<"ePi0 "<<ePi0<<endl;
	147	//cout<<"Components"<<endl;
	148	p1.SetPx(+ePi0cosPhisinThe);
	149	p1.SetPy(+ePi0sinPhisinThe);
	150	p1.SetPz(+ePi0*cosThe);
	151	p1.SetE(ePi0);
	152	//cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
	153	//cout<<"p1 Mass: "<<p1.Px()p1.Px()+p1.Py()p1.Py()+p1.Pz()p1.Pz()-p1.E()p1.E()<<endl;
	154	p2.SetPx(-ePi0cosPhisinThe);
	155	p2.SetPy(-ePi0sinPhisinThe);
	156	p2.SetPz(-ePi0*cosThe);
	157	p2.SetE(ePi0);
	158	//cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
	159	//cout<<"p2 Mass: "<<p2.Px()p2.Px()+p2.Py()p2.Py()+p2.Pz()p2.Pz()-p2.E()p2.E()<<endl;
	160	//cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl;
	161	p1.Boost(b);
	162	//cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
	163	p2.Boost(b);
	164	//cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
	165	//cout<<"angle"<<endl;
	166	angle = p1.Angle(p2.Vect());
	167	//cout<<angle<<endl;
	168	}
	169	//____________________________________________________________________________
	170	void AliPHOSGammaJet::GetGammaJet(TList &particleList, TLorentzVector &gamma, Int_t & id)
	171	{
e957fea8	172	// Get the lists of jet particles and gamma
77414c90	173	TParticle *particle = 0x0;
	174
	175	Int_t iPrimary=-1, id_motherg = -1;
	176
	177	TIter next(&particleList);
	178	while ( (particle = (TParticle*)next()) ) {
	179	iPrimary++;
	180	Int_t ksCode = particle->GetStatusCode();
	181	Int_t iMother= particle->GetMother(0);
	182
	183	if (ksCode == 21 && iMother == -1)
	184	if(particle->GetPdgCode()==22){
	185	id_motherg = iPrimary;
	186	//cout<<"idmother "<<id_motherg<<endl;
	187	}
	188	if(ksCode == 1){
	189
	190	if(id_motherg == particle->GetMother(0)){
	191	particle->Momentum(gamma);
	192	id = iPrimary;
	193	break;
	194	}
	195	}// kscode == 1
	196	}// while
	197	}
	198
	199	//____________________________________________________________________________
	200	void AliPHOSGammaJet::GetLeadingCharge(TList &particleList, TLorentzVector &charge, Int_t & id)
	201	{
e957fea8	202	// Gets the leading particle from the list of charged particles
77414c90	203	TParticle *particle = 0x0;
	204
	205	Int_t iPrimary=-1;
e957fea8	206	Double_t ptmax = 0, pti = 0;
77414c90	207	TIter next(&particleList);
	208	while ( (particle = (TParticle*)next()) ) {
	209	iPrimary++;
	210	Int_t ksCode = particle->GetStatusCode();
	211
	212	if((ksCode == 1)&&(id != iPrimary)
	213	&&(particle->GetPDG(0)->Charge()!=0)){
e957fea8	214	pti = particle->Pt();
	215	if(pti> ptmax){
	216	ptmax = pti;
77414c90	217	particle->Momentum(charge);
e957fea8	218	}//ptmax
77414c90	219	}// kscode == 1
	220	}// while
	221	}
	222
	223	//____________________________________________________________________________
90cceaf6	224	void AliPHOSGammaJet::GetLeadingPi0(TList &particleList, TLorentzVector &pi0)
77414c90	225	{
e957fea8	226	// Gets the leading pi0 from the list of particles
77414c90	227	TParticle *particle = 0x0;
	228
	229	Int_t iPrimary=-1;
e957fea8	230	Double_t ptmax = 0, pti = 0;
77414c90	231	TIter next(&particleList);
	232	while ( (particle = (TParticle*)next()) ) {
	233	iPrimary++;
	234	Int_t ksCode = particle->GetStatusCode();
	235
	236	if((ksCode == 2))
77414c90	237	{
e957fea8	238	pti = particle->Pt();
	239	if(pti> ptmax){
	240	ptmax = pti;
77414c90	241	particle->Momentum(pi0);
e957fea8	242	}//ptmax
77414c90	243	}// kscode == 1
	244	}// while
	245	}
	246
	247	//____________________________________________________________________________
90cceaf6	248	// void AliPHOSGammaJet::GetLeadingGammaPair(TList &particleList, TLorentzVector &gammapair, Int_t & id,
	249	// Double_t & thetacut,Double_t & ratiocut1, Double_t & ratiocut2,
	250	// Double_t & invmasscut1,Double_t & invmasscut2)
	251	// {
77414c90	252	// TParticle *particle = 0x0;
	253
	254	// Int_t iPrimary=-1;
	255	// Double_t El = 0, E12 = 0;
	256	// TLorentzVector gamma_i,gamma_j;
	257	// TIter next(&particleList);
	258	// while ( (particle = (TParticle*)next()) ) {
	259	// iPrimary++;
	260	// Int_t ksCode = particle->GetStatusCode();
	261	// Int_t ksPdg = particle->GetPdgCode();
	262	// if((ksCode == 1) && (iPrimary != id) && (ksPdg == 22)){
	263	// particle->Momentum(gamma_i);
	264	// Int_t jPrimary=-1;
	265	// TIter next2(&particleList);
	266	// while ( (particle = (TParticle*)next2()) ) {
	267	// jPrimary++;
	268	// if(jPrimary>iPrimary){
	269	// ksCode = particle->GetStatusCode();
	270	// ksPdg = particle->GetPdgCode();
	271	// if((ksCode == 1) && (iPrimary != id) && (ksPdg == 22)){
	272	// particle->Momentum(gamma_j);
	273	// if(gamma_j.Angle(gamma_i.Vect())<thetacut){
	274	// Float_t invmass = (gamma_i+gamma_j).M();
	275	// h_invmass->Fill(Eg,invmass);
	276	// if((invmass>invmasscut1) && (invmass<invmasscut2)){
	277	// E12 = (gamma_i+gamma_j).Energy();
	278	// if(E12>El && (E12/Eg>ratiocut1) && (E12/Eg<ratiocut2)){
	279	// //cout<<E12<<" "<<E12/Eg<<endl;
	280	// El = E12;
	281	// id_i = iPrimary;
	282	// id_j = jPrimary;
	283	// gammapair = gamma_i+gamma_j;
	284	// }//E12>El && (E12/Eg>0.2 && E12/Eg<1.)
	285	// }//(invmass>0.125) && (invmass<0.145)
	286	// }//gamma_j.Angle(gamma_i.Vect())<0.04
	287	// }//(ksCode == 1)
	288	// }
	289	// }//while
	290	// // cout<<"jPrimary "<<jPrimary<<endl;
	291	// }// if kscode 1
	292	// }//while
	293	// // cout<<"iPrimary "<<iPrimary<<endl;
90cceaf6	294	// }