[u/mrichter/AliRoot.git] / EVGEN / AliDimuCombinator.cxx

//
//
//
//
#include "AliDimuCombinator.h" 
#include "AliRun.h" 
#include "TRandom.h" 
//
ClassImp(AliDimuCombinator)
//
//                       Iterators
// 
 TParticle* AliDimuCombinator::FirstMuon()
     {
	 fimuon1=fimin1;
	 fmuon1 = (TParticle*) fPartArray->UncheckedAt(fimuon1);
	 while(Type(fmuon1)!=5 && Type(fmuon1)!=6) {
	     fimuon1++;
	     if (fimuon1 >= fimax1) {fmuon1=0; break;}
	     fmuon1 = (TParticle*) fPartArray->UncheckedAt(fimuon1);
	 }
	 return fmuon1;
     }
	     
 TParticle* AliDimuCombinator::FirstMuonSelected()
     {
	 TParticle * muon=FirstMuon();
	 while(muon!=0 && !Selected(muon)) {muon=NextMuon();}
	 return muon;
     }
	     

 TParticle* AliDimuCombinator::NextMuon()
     {
	 fimuon1++;
	 if (fimuon1>=fNParticle) {fmuon1 = 0; return fmuon1;}
	 
	 fmuon1 = (TParticle*) fPartArray->UncheckedAt(fimuon1);
	 while(Type(fmuon1)!=5 && Type(fmuon1)!=6) {
	     fimuon1++;
	     if (fimuon1>=fimax1) {fmuon1 = 0; break;}
	     fmuon1 = (TParticle*) fPartArray->UncheckedAt(fimuon1);
	 }
	 return fmuon1;
     }

TParticle* AliDimuCombinator::NextMuonSelected()
     {
	 TParticle * muon=NextMuon();
	 while(muon !=0 && !Selected(muon)) {muon=NextMuon();}
	 return muon;
     }


 void AliDimuCombinator::FirstPartner()
     {
	 if (fimin1==fimin2) {
	     fimuon2=fimuon1+1;
	 } else {
	     fimuon2=fimin2;
	 }
	 if (fimuon2 >= fimax2) {fmuon2=0; return;}
	 fmuon2 = (TParticle*) fPartArray->UncheckedAt(fimuon2);
	 while(Type(fmuon2)!=5 && Type(fmuon2)!=6) {
	     fimuon2++;
	     if (fimuon2 >= fimax2) {fmuon2=0; break;}
	     fmuon2 = (TParticle*) fPartArray->UncheckedAt(fimuon2);
	 }
     }
void AliDimuCombinator::FirstPartnerSelected()
{
	 FirstPartner();
	 while(fmuon2 !=0 && !Selected(fmuon2)) {NextPartner();}
}


 void AliDimuCombinator::NextPartner()
     {
	 fimuon2++;
	 if (fimuon2>=fimax2) {fmuon2 = 0; return;}

	 
	 fmuon2 = (TParticle*) fPartArray->UncheckedAt(fimuon2);

	 while(Type(fmuon2)!=5 && Type(fmuon2)!=6) {
	     fimuon2++;
	     if (fimuon2>=fimax2) {fmuon2 = 0; break;}
	     fmuon2 = (TParticle*) fPartArray->UncheckedAt(fimuon2);
	 }

     }

void AliDimuCombinator::NextPartnerSelected()
{
	 NextPartner();
	 while(fmuon2 !=0 && !Selected(fmuon2)) {NextPartner();}
}


 TParticle*  AliDimuCombinator::Partner()
     {
	 return fmuon2;
     }

void AliDimuCombinator::FirstMuonPair(TParticle* & muon1, TParticle* & muon2)
     {
	 FirstMuon();
	 FirstPartner();
	 muon1=fmuon1;
	 muon2=fmuon2;	 
     }
void AliDimuCombinator::NextMuonPair(TParticle* & muon1, TParticle* & muon2)
     {
	 NextPartner();
	 if (!Partner()) {
	     NextMuon();
	     FirstPartner();
	 }
	 muon1=fmuon1;
	 muon2=fmuon2;	 
     }
void AliDimuCombinator::FirstMuonPairSelected(TParticle* & muon1, TParticle* & muon2)
     {
	 FirstMuonSelected();
	 FirstPartnerSelected();
	 muon1=fmuon1;
	 muon2=fmuon2;	 
     }
void AliDimuCombinator::NextMuonPairSelected(TParticle* & muon1, TParticle* & muon2)
     {
	 NextPartnerSelected();
	 if (!Partner()) {
	     NextMuonSelected();
	     FirstPartnerSelected();
	 }
	 muon1=fmuon1;
	 muon2=fmuon2;	 
     }
void AliDimuCombinator::ResetRange()
{
    fimin1=fimin2=0;
    fimax1=fimax2=fNParticle;
}

void AliDimuCombinator::SetFirstRange(Int_t from, Int_t to)
{
    fimin1=from;
    fimax1=to;
    if (fimax1 > fNParticle) fimax1=fNParticle;
}

void AliDimuCombinator::SetSecondRange(Int_t from, Int_t to)
{
    fimin2=from;
    fimax2=to;
    if (fimax2 > fNParticle) fimax2=fNParticle;
}
//
//                       Selection
//

Bool_t AliDimuCombinator::Selected(TParticle* part)
{
// 
//
    if (part==0) {return 0;}
    
    if (part->Pt() > fPtMin && part->Eta()>fEtaMin && part->Eta()<fEtaMax) {
	return 1;
    } else {
	return 0;
    }
    
    
}

Bool_t AliDimuCombinator::Selected(TParticle* part1, TParticle* part2)
{
     return Selected(part1)*Selected(part2);
}
//
//                       Kinematics
//
Float_t AliDimuCombinator::Mass(TParticle* part1, TParticle* part2)
{
    Float_t px,py,pz,e;
    px=part1->Px()+part2->Px();
    py=part1->Py()+part2->Py();
    pz=part1->Pz()+part2->Pz();    
    e =part1->Energy()+part2->Energy();
    Float_t p=px*px+py*py+pz*pz;
    if (e*e < p) {
	return -1; 
    } else {
	return TMath::Sqrt(e*e-p);
    }
}

Float_t AliDimuCombinator::PT(TParticle* part1, TParticle* part2)
{
    Float_t px,py;
    px=part1->Px()+part2->Px();
    py=part1->Py()+part2->Py();
    return TMath::Sqrt(px*px+py*py);
}

Float_t AliDimuCombinator::Pz(TParticle* part1, TParticle* part2)
{
    return part1->Pz()+part2->Pz();
}

Float_t AliDimuCombinator::Y(TParticle* part1, TParticle* part2)
{
    Float_t pz,e;
    pz=part1->Pz()+part2->Pz();
    e =part1->Energy()+part2->Energy();
    return 0.5*TMath::Log((e+pz)/(e-pz));
}
//                  Response
//
void AliDimuCombinator::SmearGauss(Float_t width, Float_t & value)
{
    value+=gRandom->Gaus(0, width);
}
//              Weighting
// 

Float_t AliDimuCombinator::Decay_Prob(TParticle* part)
{
    Float_t d, h, theta, CTau;
    TParticle* parent = Parent(part);
    Int_t ipar=Type(parent);
    if (ipar==8 || ipar==9) {
	CTau=780.4;
    } else if (ipar==11 || ipar==12) {
	CTau=370.9;
    } else {
	CTau=0;
    }
    
    
    Float_t GammaBeta=(parent->P())/(parent->GetMass());
//
// this part is still very ALICE muon-arm specific
//
    theta=parent->Theta();
    h=90*TMath::Tan(theta);
    
    if (h<4) {
	d=4/TMath::Sin(theta);
    } else {
	d=90/TMath::Cos(theta);
    }
    
    if (CTau > 0) {
	return 1-TMath::Exp(-d/CTau/GammaBeta);
    } else {
	return 1;
    }
}

Float_t AliDimuCombinator::Weight(TParticle* part1, TParticle* part2)
{
    Float_t wgt=(part1->GetWeight())*(part2->GetWeight());
    
    if (Correlated(part1, part2)) {
	return wgt/(Parent(part1)->GetWeight())*fRate1;
    } else {
	return wgt*fRate1*fRate2;
    }
} 


Float_t AliDimuCombinator::Weight(TParticle* part)
{
    return (part->GetWeight())*(Parent(part)->GetWeight())*fRate1;
}
Bool_t  AliDimuCombinator::Correlated(TParticle* part1, TParticle* part2)
{
    if (Origin(part1) == Origin(part2)) {
	return kTRUE;
    } else {
	return kFALSE;
    }
}

TParticle* AliDimuCombinator::Parent(TParticle* part)
{
    return (TParticle*) (fPartArray->UncheckedAt(part->GetFirstMother()));
}

Int_t AliDimuCombinator::Origin(TParticle* part)
{
    Int_t iparent= part->GetFirstMother();
    if (iparent < 0) return iparent;
    Int_t ip;
    while(1) {
	ip=((TParticle*) fPartArray->UncheckedAt(iparent))->GetFirstMother();
	if (ip < 0) {
	    break;
	} else {
	    iparent=ip;
	}
    }
    return iparent;
}
Commit	Line	Data
	1	//
	2	//
	3	//
	4	//
	5	#include "AliDimuCombinator.h"
	6	#include "AliRun.h"
	7	#include "TRandom.h"
	8	//
	9	ClassImp(AliDimuCombinator)
	10	//
	11	// Iterators
	12	//
	13	TParticle* AliDimuCombinator::FirstMuon()
	14	{
	15	fimuon1=fimin1;
	16	fmuon1 = (TParticle*) fPartArray->UncheckedAt(fimuon1);
	17	while(Type(fmuon1)!=5 && Type(fmuon1)!=6) {
	18	fimuon1++;
	19	if (fimuon1 >= fimax1) {fmuon1=0; break;}
	20	fmuon1 = (TParticle*) fPartArray->UncheckedAt(fimuon1);
	21	}
	22	return fmuon1;
	23	}
	24
	25	TParticle* AliDimuCombinator::FirstMuonSelected()
	26	{
	27	TParticle * muon=FirstMuon();
	28	while(muon!=0 && !Selected(muon)) {muon=NextMuon();}
	29	return muon;
	30	}
	31
	32
	33	TParticle* AliDimuCombinator::NextMuon()
	34	{
	35	fimuon1++;
	36	if (fimuon1>=fNParticle) {fmuon1 = 0; return fmuon1;}
	37
	38	fmuon1 = (TParticle*) fPartArray->UncheckedAt(fimuon1);
	39	while(Type(fmuon1)!=5 && Type(fmuon1)!=6) {
	40	fimuon1++;
	41	if (fimuon1>=fimax1) {fmuon1 = 0; break;}
	42	fmuon1 = (TParticle*) fPartArray->UncheckedAt(fimuon1);
	43	}
	44	return fmuon1;
	45	}
	46
	47	TParticle* AliDimuCombinator::NextMuonSelected()
	48	{
	49	TParticle * muon=NextMuon();
	50	while(muon !=0 && !Selected(muon)) {muon=NextMuon();}
	51	return muon;
	52	}
	53
	54
	55	void AliDimuCombinator::FirstPartner()
	56	{
	57	if (fimin1==fimin2) {
	58	fimuon2=fimuon1+1;
	59	} else {
	60	fimuon2=fimin2;
	61	}
	62	if (fimuon2 >= fimax2) {fmuon2=0; return;}
	63	fmuon2 = (TParticle*) fPartArray->UncheckedAt(fimuon2);
	64	while(Type(fmuon2)!=5 && Type(fmuon2)!=6) {
	65	fimuon2++;
	66	if (fimuon2 >= fimax2) {fmuon2=0; break;}
	67	fmuon2 = (TParticle*) fPartArray->UncheckedAt(fimuon2);
	68	}
	69	}
	70	void AliDimuCombinator::FirstPartnerSelected()
	71	{
	72	FirstPartner();
	73	while(fmuon2 !=0 && !Selected(fmuon2)) {NextPartner();}
	74	}
	75
	76
	77	void AliDimuCombinator::NextPartner()
	78	{
	79	fimuon2++;
	80	if (fimuon2>=fimax2) {fmuon2 = 0; return;}
	81
	82
	83	fmuon2 = (TParticle*) fPartArray->UncheckedAt(fimuon2);
	84
	85	while(Type(fmuon2)!=5 && Type(fmuon2)!=6) {
	86	fimuon2++;
	87	if (fimuon2>=fimax2) {fmuon2 = 0; break;}
	88	fmuon2 = (TParticle*) fPartArray->UncheckedAt(fimuon2);
	89	}
	90
	91	}
	92
	93	void AliDimuCombinator::NextPartnerSelected()
	94	{
	95	NextPartner();
	96	while(fmuon2 !=0 && !Selected(fmuon2)) {NextPartner();}
	97	}
	98
	99
	100	TParticle* AliDimuCombinator::Partner()
	101	{
	102	return fmuon2;
	103	}
	104
	105	void AliDimuCombinator::FirstMuonPair(TParticle* & muon1, TParticle* & muon2)
	106	{
	107	FirstMuon();
	108	FirstPartner();
	109	muon1=fmuon1;
	110	muon2=fmuon2;
	111	}
	112	void AliDimuCombinator::NextMuonPair(TParticle* & muon1, TParticle* & muon2)
	113	{
	114	NextPartner();
	115	if (!Partner()) {
	116	NextMuon();
	117	FirstPartner();
	118	}
	119	muon1=fmuon1;
	120	muon2=fmuon2;
	121	}
	122	void AliDimuCombinator::FirstMuonPairSelected(TParticle* & muon1, TParticle* & muon2)
	123	{
	124	FirstMuonSelected();
	125	FirstPartnerSelected();
	126	muon1=fmuon1;
	127	muon2=fmuon2;
	128	}
	129	void AliDimuCombinator::NextMuonPairSelected(TParticle* & muon1, TParticle* & muon2)
	130	{
	131	NextPartnerSelected();
	132	if (!Partner()) {
	133	NextMuonSelected();
	134	FirstPartnerSelected();
	135	}
	136	muon1=fmuon1;
	137	muon2=fmuon2;
	138	}
	139	void AliDimuCombinator::ResetRange()
	140	{
	141	fimin1=fimin2=0;
	142	fimax1=fimax2=fNParticle;
	143	}
	144
	145	void AliDimuCombinator::SetFirstRange(Int_t from, Int_t to)
	146	{
	147	fimin1=from;
	148	fimax1=to;
	149	if (fimax1 > fNParticle) fimax1=fNParticle;
	150	}
	151
	152	void AliDimuCombinator::SetSecondRange(Int_t from, Int_t to)
	153	{
	154	fimin2=from;
	155	fimax2=to;
	156	if (fimax2 > fNParticle) fimax2=fNParticle;
	157	}
	158	//
	159	// Selection
	160	//
	161
	162	Bool_t AliDimuCombinator::Selected(TParticle* part)
	163	{
	164	//
	165	//
	166	if (part==0) {return 0;}
	167
	168	if (part->Pt() > fPtMin && part->Eta()>fEtaMin && part->Eta()<fEtaMax) {
	169	return 1;
	170	} else {
	171	return 0;
	172	}
	173
	174
	175	}
	176
	177	Bool_t AliDimuCombinator::Selected(TParticle* part1, TParticle* part2)
	178	{
	179	return Selected(part1)*Selected(part2);
	180	}
	181	//
	182	// Kinematics
	183	//
	184	Float_t AliDimuCombinator::Mass(TParticle* part1, TParticle* part2)
	185	{
	186	Float_t px,py,pz,e;
	187	px=part1->Px()+part2->Px();
	188	py=part1->Py()+part2->Py();
	189	pz=part1->Pz()+part2->Pz();
	190	e =part1->Energy()+part2->Energy();
	191	Float_t p=pxpx+pypy+pz*pz;
	192	if (e*e < p) {
	193	return -1;
	194	} else {
	195	return TMath::Sqrt(e*e-p);
	196	}
	197	}
	198
	199	Float_t AliDimuCombinator::PT(TParticle* part1, TParticle* part2)
	200	{
	201	Float_t px,py;
	202	px=part1->Px()+part2->Px();
	203	py=part1->Py()+part2->Py();
	204	return TMath::Sqrt(pxpx+pypy);
	205	}
	206
	207	Float_t AliDimuCombinator::Pz(TParticle* part1, TParticle* part2)
	208	{
	209	return part1->Pz()+part2->Pz();
	210	}
	211
	212	Float_t AliDimuCombinator::Y(TParticle* part1, TParticle* part2)
	213	{
	214	Float_t pz,e;
	215	pz=part1->Pz()+part2->Pz();
	216	e =part1->Energy()+part2->Energy();
	217	return 0.5*TMath::Log((e+pz)/(e-pz));
	218	}
	219	// Response
	220	//
	221	void AliDimuCombinator::SmearGauss(Float_t width, Float_t & value)
	222	{
	223	value+=gRandom->Gaus(0, width);
	224	}
	225	// Weighting
	226	//
	227
	228	Float_t AliDimuCombinator::Decay_Prob(TParticle* part)
	229	{
	230	Float_t d, h, theta, CTau;
	231	TParticle* parent = Parent(part);
	232	Int_t ipar=Type(parent);
	233	if (ipar==8 \|\| ipar==9) {
	234	CTau=780.4;
	235	} else if (ipar==11 \|\| ipar==12) {
	236	CTau=370.9;
	237	} else {
	238	CTau=0;
	239	}
	240
	241
	242	Float_t GammaBeta=(parent->P())/(parent->GetMass());
	243	//
	244	// this part is still very ALICE muon-arm specific
	245	//
	246	theta=parent->Theta();
	247	h=90*TMath::Tan(theta);
	248
	249	if (h<4) {
	250	d=4/TMath::Sin(theta);
	251	} else {
	252	d=90/TMath::Cos(theta);
	253	}
	254
	255	if (CTau > 0) {
	256	return 1-TMath::Exp(-d/CTau/GammaBeta);
	257	} else {
	258	return 1;
	259	}
	260	}
	261
	262	Float_t AliDimuCombinator::Weight(TParticle* part1, TParticle* part2)
	263	{
	264	Float_t wgt=(part1->GetWeight())*(part2->GetWeight());
	265
	266	if (Correlated(part1, part2)) {
	267	return wgt/(Parent(part1)->GetWeight())*fRate1;
	268	} else {
	269	return wgtfRate1fRate2;
	270	}
	271	}
	272
	273
	274	Float_t AliDimuCombinator::Weight(TParticle* part)
	275	{
	276	return (part->GetWeight())(Parent(part)->GetWeight())fRate1;
	277	}
	278	Bool_t AliDimuCombinator::Correlated(TParticle* part1, TParticle* part2)
	279	{
	280	if (Origin(part1) == Origin(part2)) {
	281	return kTRUE;
	282	} else {
	283	return kFALSE;
	284	}
	285	}
	286
	287	TParticle* AliDimuCombinator::Parent(TParticle* part)
	288	{
	289	return (TParticle*) (fPartArray->UncheckedAt(part->GetFirstMother()));
	290	}
	291
	292	Int_t AliDimuCombinator::Origin(TParticle* part)
	293	{
	294	Int_t iparent= part->GetFirstMother();
	295	if (iparent < 0) return iparent;
	296	Int_t ip;
	297	while(1) {
	298	ip=((TParticle*) fPartArray->UncheckedAt(iparent))->GetFirstMother();
	299	if (ip < 0) {
	300	break;
	301	} else {
	302	iparent=ip;
	303	}
	304	}
	305	return iparent;
	306	}
	307