[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 "TString.h"
#include "TFile.h"
#include "TLorentzVector.h"
#include "TList.h"
#include "TTree.h"
#include "TParticle.h"
#include "AliRun.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 *) 
{
  // does the job
  
  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 id_g = -1;
    GetGammaJet(particleList,gamma, id_g);
    GetLeadingCharge(particleList,charge, id_g);
    GetLeadingPi0(particleList,pi0);
    Info("Pi0Decay", "Gamma: %f %d", gamma.Energy(), id_g) ;
    Info("Pi0Decay", "Charge: %f", charge.Energy()) ;
    Info("Pi0Decay", "Pi0: %f", pi0.Energy()) ;
    //    GetLeadingGammaPair(particleList, gammapair, id_g, 
    //			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) 
{
  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) 
{
  TParticle *particle = 0x0;
  
  Int_t iPrimary=-1;
  Double_t pt_max = 0, pt_i = 0;
  TIter next(&particleList);
  while ( (particle = (TParticle*)next()) ) {
    iPrimary++;  
    Int_t ksCode = particle->GetStatusCode();
    
    if((ksCode == 1)&&(id != iPrimary)
       &&(particle->GetPDG(0)->Charge()!=0)){
      pt_i = particle->Pt(); 
      if(pt_i> pt_max){
	pt_max = pt_i;
 	particle->Momentum(charge);
      }//pt_max   
    }// kscode == 1
  }// while
}

//____________________________________________________________________________
void AliPHOSGammaJet::GetLeadingPi0(TList &particleList, TLorentzVector &pi0) 
{
  TParticle *particle = 0x0;
  
  Int_t iPrimary=-1;
  Double_t pt_max = 0, pt_i = 0;
  TIter next(&particleList);
  while ( (particle = (TParticle*)next()) ) {
    iPrimary++;  
    Int_t ksCode = particle->GetStatusCode();
    
    if((ksCode == 2))
      {
	pt_i = particle->Pt(); 
	if(pt_i> pt_max){
	  pt_max = pt_i;
	  particle->Momentum(pi0);
	}//pt_max   
      }// 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 ---
	26
	27	#include "TString.h"
	28	#include "TFile.h"
	29	#include "TLorentzVector.h"
	30	#include "TList.h"
	31	#include "TTree.h"
	32	#include "TParticle.h"
	33	#include "AliRun.h"
	34	#include "AliPHOSGammaJet.h"
	35	#include "AliPHOSGetter.h"
	36
	37	ClassImp(AliPHOSGammaJet)
	38
	39	//____________________________________________________________________________
	40	AliPHOSGammaJet::AliPHOSGammaJet() {
	41	// ctor
	42	}
	43
	44	//____________________________________________________________________________
	45	AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) {
	46	// ctor
	47	AliPHOSGetter::Instance(inputfilename) ;
	48	}
	49
	50	//____________________________________________________________________________
90cceaf6	51	AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) : TTask(gj) {
77414c90	52	}
	53
	54	//____________________________________________________________________________
	55	AliPHOSGammaJet::~AliPHOSGammaJet() {
	56	}
	57
	58	//____________________________________________________________________________
90cceaf6	59	void AliPHOSGammaJet::Exec(Option_t *)
77414c90	60	{
	61	// does the job
	62
	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;
	117	Int_t id_g = -1;
	118	GetGammaJet(particleList,gamma, id_g);
	119	GetLeadingCharge(particleList,charge, id_g);
90cceaf6	120	GetLeadingPi0(particleList,pi0);
77414c90	121	Info("Pi0Decay", "Gamma: %f %d", gamma.Energy(), id_g) ;
	122	Info("Pi0Decay", "Charge: %f", charge.Energy()) ;
	123	Info("Pi0Decay", "Pi0: %f", pi0.Energy()) ;
	124	// GetLeadingGammaPair(particleList, gammapair, id_g,
	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	{
	172	TParticle *particle = 0x0;
	173
	174	Int_t iPrimary=-1, id_motherg = -1;
	175
	176	TIter next(&particleList);
	177	while ( (particle = (TParticle*)next()) ) {
	178	iPrimary++;
	179	Int_t ksCode = particle->GetStatusCode();
	180	Int_t iMother= particle->GetMother(0);
	181
	182	if (ksCode == 21 && iMother == -1)
	183	if(particle->GetPdgCode()==22){
	184	id_motherg = iPrimary;
185	//cout<<"idmother "<<id_motherg<<endl;
186	}
187	if(ksCode == 1){
188
189	if(id_motherg == particle->GetMother(0)){
190	particle->Momentum(gamma);
191	id = iPrimary;
192	break;
193	}
194	}// kscode == 1
195	}// while
196	}
197
198	//____________________________________________________________________________
199	void AliPHOSGammaJet::GetLeadingCharge(TList &particleList, TLorentzVector &charge, Int_t & id)
200	{
201	TParticle *particle = 0x0;
202
203	Int_t iPrimary=-1;
204	Double_t pt_max = 0, pt_i = 0;
205	TIter next(&particleList);
206	while ( (particle = (TParticle*)next()) ) {
207	iPrimary++;
208	Int_t ksCode = particle->GetStatusCode();
209
210	if((ksCode == 1)&&(id != iPrimary)
211	&&(particle->GetPDG(0)->Charge()!=0)){
212	pt_i = particle->Pt();
213	if(pt_i> pt_max){
214	pt_max = pt_i;
215	particle->Momentum(charge);
216	}//pt_max
217	}// kscode == 1
218	}// while
219	}
220
221	//____________________________________________________________________________
90cceaf6	222	void AliPHOSGammaJet::GetLeadingPi0(TList &particleList, TLorentzVector &pi0)
77414c90	223	{
	224	TParticle *particle = 0x0;
	225
	226	Int_t iPrimary=-1;
	227	Double_t pt_max = 0, pt_i = 0;
	228	TIter next(&particleList);
	229	while ( (particle = (TParticle*)next()) ) {
	230	iPrimary++;
	231	Int_t ksCode = particle->GetStatusCode();
	232
	233	if((ksCode == 2))
77414c90	234	{
	235	pt_i = particle->Pt();
	236	if(pt_i> pt_max){
	237	pt_max = pt_i;
	238	particle->Momentum(pi0);
	239	}//pt_max
	240	}// kscode == 1
	241	}// while
	242	}
	243
	244	//____________________________________________________________________________
90cceaf6	245	// void AliPHOSGammaJet::GetLeadingGammaPair(TList &particleList, TLorentzVector &gammapair, Int_t & id,
	246	// Double_t & thetacut,Double_t & ratiocut1, Double_t & ratiocut2,
	247	// Double_t & invmasscut1,Double_t & invmasscut2)
	248	// {
77414c90	249	// TParticle *particle = 0x0;
	250
	251	// Int_t iPrimary=-1;
	252	// Double_t El = 0, E12 = 0;
	253	// TLorentzVector gamma_i,gamma_j;
	254	// TIter next(&particleList);
	255	// while ( (particle = (TParticle*)next()) ) {
	256	// iPrimary++;
	257	// Int_t ksCode = particle->GetStatusCode();
	258	// Int_t ksPdg = particle->GetPdgCode();
	259	// if((ksCode == 1) && (iPrimary != id) && (ksPdg == 22)){
	260	// particle->Momentum(gamma_i);
	261	// Int_t jPrimary=-1;
	262	// TIter next2(&particleList);
	263	// while ( (particle = (TParticle*)next2()) ) {
	264	// jPrimary++;
	265	// if(jPrimary>iPrimary){
	266	// ksCode = particle->GetStatusCode();
	267	// ksPdg = particle->GetPdgCode();
	268	// if((ksCode == 1) && (iPrimary != id) && (ksPdg == 22)){
	269	// particle->Momentum(gamma_j);
	270	// if(gamma_j.Angle(gamma_i.Vect())<thetacut){
	271	// Float_t invmass = (gamma_i+gamma_j).M();
	272	// h_invmass->Fill(Eg,invmass);
	273	// if((invmass>invmasscut1) && (invmass<invmasscut2)){
	274	// E12 = (gamma_i+gamma_j).Energy();
	275	// if(E12>El && (E12/Eg>ratiocut1) && (E12/Eg<ratiocut2)){
	276	// //cout<<E12<<" "<<E12/Eg<<endl;
	277	// El = E12;
	278	// id_i = iPrimary;
	279	// id_j = jPrimary;
	280	// gammapair = gamma_i+gamma_j;
	281	// }//E12>El && (E12/Eg>0.2 && E12/Eg<1.)
	282	// }//(invmass>0.125) && (invmass<0.145)
	283	// }//gamma_j.Angle(gamma_i.Vect())<0.04
	284	// }//(ksCode == 1)
	285	// }
	286	// }//while
	287	// // cout<<"jPrimary "<<jPrimary<<endl;
	288	// }// if kscode 1
	289	// }//while
	290	// // cout<<"iPrimary "<<iPrimary<<endl;
90cceaf6	291	// }