[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliAODConversionMother.cxx

#include "AliAODConversionMother.h"
#include "AliKFConversionMother.h"
#include "AliKFParticle.h"


// Author D. Lohner (Daniel.Lohner@cern.ch)

using namespace std;

ClassImp(AliAODConversionMother)

AliAODConversionMother::AliAODConversionMother() :
AliAODConversionParticle(),
    fMCLabel(-1),
    fChi2(-1),
    fOpeningAngle(-1),
    fAlpha(-1),
    fWeight(1),
    fdcaBetweenPhotons(1),
    fdcaZPrimVtx(100),
    fdcaRPrimVtx(100),
    fQuality(0),
    fTrueMeson(0)
   
{
	fLabel[0] = -1;
	fLabel[1] = -1;
	fLabel[2] = 0;
   
   fProductionVtx[0]=0;
   fProductionVtx[1]=0;
   fProductionVtx[2]=0;
  
}

AliAODConversionMother::AliAODConversionMother(AliKFConversionMother *kf):
AliAODConversionParticle(),
fMCLabel(kf->GetMCLabel()),
fChi2(kf->GetChi2()),
fOpeningAngle(kf->GetOpeningAngle()),
fAlpha(kf->GetAlpha()),
fWeight(1),
fdcaBetweenPhotons(100),
fdcaZPrimVtx(100),
fdcaRPrimVtx(100),
fQuality(0),
fTrueMeson(0)
{
    // Set 4momentu
    SetPxPyPzE(kf->GetPx(),kf->GetPy(),kf->GetPz(),kf->GetE());

   fProductionVtx[0]=0;
   fProductionVtx[1]=0;
   fProductionVtx[2]=0;
  
     //Set Decay Photon Labels
    fLabel[0]=kf->GetGammaLabel(0);
    fLabel[1]=kf->GetGammaLabel(1);
    fLabel[2]=0;
}

AliAODConversionMother::AliAODConversionMother(AliAODConversionPhoton *y1,AliAODConversionPhoton *y2):
AliAODConversionParticle(),
fMCLabel(-1),
fChi2(-1),
fOpeningAngle(-1),
fAlpha(-1),
fWeight(1),
fdcaBetweenPhotons(1),
fdcaZPrimVtx(100),
fdcaRPrimVtx(100),
fQuality(0),
fTrueMeson(0)
{
    // Set 4momentum
    SetPxPyPzE(y1->Px()+y2->Px(),y1->Py()+y2->Py(),y1->Pz()+y2->Pz(),y1->E()+y2->E());

    // Calculate Opening Angle
    TVector3 v1(y1->Px(),y1->Py(),y1->Pz());
    TVector3 v2(y2->Px(),y2->Py(),y2->Pz());
    fOpeningAngle=v1.Angle(v2);
    fdcaBetweenPhotons = CalculateDistanceBetweenPhotons(y1,y2,fProductionVtx);
    DetermineMesonQuality(y1,y2);
    // Calculate Alpha
    if((y1->E()+y2->E()) != 0){
		fAlpha=TMath::Abs((y1->E()-y2->E())/(y1->E()+y2->E()));
    }

    // Set Chi2 to the mean chi2 of gammas
 //   fChi2=0.5*(y1->GetChi2perNDF()+y2->GetChi2perNDF());

    //Set Decay Photon Labels
    fLabel[0]=-1;
    fLabel[1]=-1;
    fLabel[2]=0;
}

AliAODConversionMother::AliAODConversionMother(AliAODConversionMother *meson,AliAODConversionPhoton *gamma):
AliAODConversionParticle(),
fMCLabel(-1),
fChi2(-1),
fOpeningAngle(-1),
fAlpha(-1),
fWeight(1),
fdcaBetweenPhotons(1),
fdcaZPrimVtx(100),
fdcaRPrimVtx(100),
fQuality(0),
fTrueMeson(0)
{
    // Set 4momentum
    SetPxPyPzE(meson->Px()+gamma->Px(),meson->Py()+gamma->Py(),meson->Pz()+gamma->Pz(),meson->E()+gamma->E());

    // Calculate Opening Angle
    TVector3 v1(meson->Px(),meson->Py(),meson->Pz());
    TVector3 v2(gamma->Px(),gamma->Py(),gamma->Pz());
    fOpeningAngle=v1.Angle(v2);
     
	fProductionVtx[0]=0;
	fProductionVtx[1]=0;
	fProductionVtx[2]=0;

    // Calculate Alpha
    if((meson->E()+gamma->E()) != 0){
		fAlpha=TMath::Abs((meson->E()-gamma->E())/(meson->E()+gamma->E()));
    }

    // Set Chi2 to the mean chi2 of gammas
	// fChi2=0.5*(y1->GetChi2perNDF()+y2->GetChi2perNDF());

    //Set Decay Photon Labels
    fLabel[0]=-1;
    fLabel[1]=-1;
    fLabel[2]=0;
}


AliAODConversionMother::~AliAODConversionMother() {
    // empty standard destructor
}

TParticle *AliAODConversionMother::GetMCParticle(AliStack *fMCStack){
    if(!fMCStack){AliError("MC Stack not defined");return 0x0;}

    if(fMCLabel>-1){
		return fMCStack->Particle(fMCLabel);
    }
    return 0x0;
}

Bool_t AliAODConversionMother::IsTrueMeson(AliStack *fMCStack,Int_t pdgcode){
    TParticle *part=GetMCParticle(fMCStack);

    if(part){
	// Check if it is a true photon
	if(part->GetPdgCode()==pdgcode){
	    return kTRUE;
	}
    }
    return kFALSE;
}

Float_t AliAODConversionMother::CalculateDistanceBetweenPhotons(AliAODConversionPhoton* y1, AliAODConversionPhoton* y2 , Double_t prodPoint[3]){

   TVector3 a(y1->GetConversionX(),y1->GetConversionY(),y1->GetConversionZ());
   TVector3 b(y1->GetPx(),y1->GetPy(),y1->GetPz());
   TVector3 c(y2->GetConversionX(),y2->GetConversionY(),y2->GetConversionZ());
   TVector3 d(y2->GetPx(),y2->GetPy(),y2->GetPz());
   
   TVector3 n = b.Cross(d);
   TVector3 nn = n.Unit();
   
   Double_t dist = 0;
   if (n.Mag() == 0){
      TVector3 e = a-c;
      if (d.Mag() != 0){
         dist = TMath::Abs((e.Cross(d)).Mag())/TMath::Abs(d.Mag());
      }
      prodPoint[0] = 0;
      prodPoint[1] = 0;
      prodPoint[2] = 0;
   } else {   
      dist = TMath::Abs(n.Dot(c-a))/TMath::Abs(n.Mag());
      Double_t lambda = (b.Dot(d) * (a-c).Dot(d) - d.Dot(d) * (a-c).Dot(b))/(b.Dot(b) * d.Dot(d) - TMath::Power(b.Dot(d),2));
      Double_t mu = ((a-c).Dot(d) * b.Dot(b) - (a-c).Dot(b) * b.Dot(d) )/(b.Dot(b) * d.Dot(d) - TMath::Power(b.Dot(d),2));
      
      TVector3 S1 = a + lambda* b;
      TVector3 S2 = c + mu* d;
      TVector3 Prod = S1 + 0.5*dist*(S2-S1).Unit();
      prodPoint[0] = Prod(0);
      prodPoint[1] = Prod(1);
      prodPoint[2] = Prod(2);
      
   }
   if (dist > 1000) dist = 999.;
   return dist;
}

///________________________________________________________________________
void AliAODConversionMother::CalculateDistanceOfClossetApproachToPrimVtx(const AliVVertex* primVertex){

   Double_t primCo[3] = {primVertex->GetX(),primVertex->GetY(),primVertex->GetZ()};

   Double_t absoluteP = TMath::Sqrt(TMath::Power(Px(),2) + TMath::Power(Py(),2) + TMath::Power(Pz(),2));   
   Double_t p[3] = {Px()/absoluteP,Py()/absoluteP,Pz()/absoluteP};
   Double_t CP[3];
   
   CP[0] =  fProductionVtx[0] - primCo[0];
   CP[1] =  fProductionVtx[1] - primCo[1];
   CP[2] =  fProductionVtx[2] - primCo[2];
   
   Double_t Lambda = - (CP[0]*p[0]+CP[1]*p[1]+CP[2]*p[2])/(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
   
   Double_t S[3];
   S[0] = fProductionVtx[0] + p[0]*Lambda;
   S[1] = fProductionVtx[1] + p[1]*Lambda;
   S[2] = fProductionVtx[2] + p[2]*Lambda;
   
   fdcaRPrimVtx = TMath::Sqrt( TMath::Power(primCo[0]-S[0],2) + TMath::Power(primCo[1]-S[1],2));
   fdcaZPrimVtx = primCo[2]-S[2];
   
   
//    cout << "DCA z: " << dca[1] << "\t DCA r: " << dca[0] << "\t DCA 3d: " << TMath::Sqrt(dca[1]*dca[1] + dca[0]*dca[0]) << endl;
   
   
   return;
}

///________________________________________________________________________
void AliAODConversionMother::DetermineMesonQuality(AliAODConversionPhoton* y1, AliAODConversionPhoton* y2){
   UChar_t photonQA1 = y1->GetPhotonQuality();
   UChar_t photonQA2 = y2->GetPhotonQuality();
   
   if (photonQA1 == 0 || photonQA2 == 0){
      fQuality = 0;
      return;
   } 
   if (photonQA1 == 1 && photonQA2 == 1){
      fQuality = 1;
      return;
   }               
   if (photonQA1 == 2 && photonQA2 == 2){
      fQuality = 4;
      return;
   }               
   if (photonQA1 == 3 && photonQA2 == 3){
      fQuality = 6;
      return;
   }         
   if (photonQA1 == 1){
       if (photonQA2 == 2){
        fQuality = 2;  
        return;
       }
       if (photonQA2 == 3){
        fQuality = 3;  
        return;
       }
   }
   if (photonQA2 == 1){
       if (photonQA1 == 2){
        fQuality = 2;  
        return;
       }
       if (photonQA1 == 3){
        fQuality = 3;  
        return;
       }
   }
   if ((photonQA1 == 2 && photonQA2 == 3)|| (photonQA1 == 3 && photonQA2 == 2)){
        fQuality = 5;  
        return;
   }
   
}
Commit	Line	Data
2eedd4ed	1	#include "AliAODConversionMother.h"
	2	#include "AliKFConversionMother.h"
	3	#include "AliKFParticle.h"
	4
	5
	6	// Author D. Lohner (Daniel.Lohner@cern.ch)
	7
	8	using namespace std;
	9
	10	ClassImp(AliAODConversionMother)
	11
	12	AliAODConversionMother::AliAODConversionMother() :
	13	AliAODConversionParticle(),
92efd725	14	fMCLabel(-1),
	15	fChi2(-1),
	16	fOpeningAngle(-1),
69b05e50	17	fAlpha(-1),
4803eb1f	18	fWeight(1),
	19	fdcaBetweenPhotons(1),
	20	fdcaZPrimVtx(100),
	21	fdcaRPrimVtx(100),
eed82f0e	22	fQuality(0),
eed82f0e	23	fTrueMeson(0)
4803eb1f	24
2eedd4ed	25	{
fd420a4f	26	fLabel[0] = -1;
fd420a4f	27	fLabel[1] = -1;
92efd725	28	fLabel[2] = 0;
4803eb1f	29
	30	fProductionVtx[0]=0;
	31	fProductionVtx[1]=0;
	32	fProductionVtx[2]=0;
	33
2eedd4ed	34	}
	35
	36	AliAODConversionMother::AliAODConversionMother(AliKFConversionMother *kf):
	37	AliAODConversionParticle(),
92efd725	38	fMCLabel(kf->GetMCLabel()),
2eedd4ed	39	fChi2(kf->GetChi2()),
2eedd4ed	40	fOpeningAngle(kf->GetOpeningAngle()),
69b05e50	41	fAlpha(kf->GetAlpha()),
4803eb1f	42	fWeight(1),
	43	fdcaBetweenPhotons(100),
	44	fdcaZPrimVtx(100),
	45	fdcaRPrimVtx(100),
eed82f0e	46	fQuality(0),
eed82f0e	47	fTrueMeson(0)
2eedd4ed	48	{
	49	// Set 4momentu
	50	SetPxPyPzE(kf->GetPx(),kf->GetPy(),kf->GetPz(),kf->GetE());
	51
4803eb1f	52	fProductionVtx[0]=0;
	53	fProductionVtx[1]=0;
	54	fProductionVtx[2]=0;
	55
2eedd4ed	56	//Set Decay Photon Labels
	57	fLabel[0]=kf->GetGammaLabel(0);
	58	fLabel[1]=kf->GetGammaLabel(1);
92efd725	59	fLabel[2]=0;
2eedd4ed	60	}
	61
	62	AliAODConversionMother::AliAODConversionMother(AliAODConversionPhoton y1,AliAODConversionPhoton y2):
	63	AliAODConversionParticle(),
92efd725	64	fMCLabel(-1),
2eedd4ed	65	fChi2(-1),
2eedd4ed	66	fOpeningAngle(-1),
69b05e50	67	fAlpha(-1),
4803eb1f	68	fWeight(1),
	69	fdcaBetweenPhotons(1),
	70	fdcaZPrimVtx(100),
	71	fdcaRPrimVtx(100),
eed82f0e	72	fQuality(0),
eed82f0e	73	fTrueMeson(0)
2eedd4ed	74	{
	75	// Set 4momentum
	76	SetPxPyPzE(y1->Px()+y2->Px(),y1->Py()+y2->Py(),y1->Pz()+y2->Pz(),y1->E()+y2->E());
	77
	78	// Calculate Opening Angle
	79	TVector3 v1(y1->Px(),y1->Py(),y1->Pz());
	80	TVector3 v2(y2->Px(),y2->Py(),y2->Pz());
	81	fOpeningAngle=v1.Angle(v2);
4803eb1f	82	fdcaBetweenPhotons = CalculateDistanceBetweenPhotons(y1,y2,fProductionVtx);
4803eb1f	83	DetermineMesonQuality(y1,y2);
2eedd4ed	84	// Calculate Alpha
2eedd4ed	85	if((y1->E()+y2->E()) != 0){
eed82f0e	86	fAlpha=TMath::Abs((y1->E()-y2->E())/(y1->E()+y2->E()));
2eedd4ed	87	}
	88
	89	// Set Chi2 to the mean chi2 of gammas
3b77b2d1	90	// fChi2=0.5*(y1->GetChi2perNDF()+y2->GetChi2perNDF());
2eedd4ed	91
	92	//Set Decay Photon Labels
	93	fLabel[0]=-1;
	94	fLabel[1]=-1;
a280ac15	95	fLabel[2]=0;
2eedd4ed	96	}
2eedd4ed	97
cb287740	98	AliAODConversionMother::AliAODConversionMother(AliAODConversionMother meson,AliAODConversionPhoton gamma):
	99	AliAODConversionParticle(),
	100	fMCLabel(-1),
	101	fChi2(-1),
	102	fOpeningAngle(-1),
	103	fAlpha(-1),
	104	fWeight(1),
	105	fdcaBetweenPhotons(1),
	106	fdcaZPrimVtx(100),
	107	fdcaRPrimVtx(100),
eed82f0e	108	fQuality(0),
eed82f0e	109	fTrueMeson(0)
cb287740	110	{
	111	// Set 4momentum
	112	SetPxPyPzE(meson->Px()+gamma->Px(),meson->Py()+gamma->Py(),meson->Pz()+gamma->Pz(),meson->E()+gamma->E());
	113
	114	// Calculate Opening Angle
	115	TVector3 v1(meson->Px(),meson->Py(),meson->Pz());
	116	TVector3 v2(gamma->Px(),gamma->Py(),gamma->Pz());
	117	fOpeningAngle=v1.Angle(v2);
	118
	119	fProductionVtx[0]=0;
	120	fProductionVtx[1]=0;
	121	fProductionVtx[2]=0;
	122
	123	// Calculate Alpha
	124	if((meson->E()+gamma->E()) != 0){
	125	fAlpha=TMath::Abs((meson->E()-gamma->E())/(meson->E()+gamma->E()));
	126	}
	127
	128	// Set Chi2 to the mean chi2 of gammas
eed82f0e	129	// fChi2=0.5*(y1->GetChi2perNDF()+y2->GetChi2perNDF());
cb287740	130
	131	//Set Decay Photon Labels
	132	fLabel[0]=-1;
	133	fLabel[1]=-1;
	134	fLabel[2]=0;
	135	}
	136
	137
2eedd4ed	138	AliAODConversionMother::~AliAODConversionMother() {
92efd725	139	// empty standard destructor
	140	}
	141
	142	TParticle AliAODConversionMother::GetMCParticle(AliStack fMCStack){
	143	if(!fMCStack){AliError("MC Stack not defined");return 0x0;}
2eedd4ed	144
92efd725	145	if(fMCLabel>-1){
	146	return fMCStack->Particle(fMCLabel);
	147	}
	148	return 0x0;
	149	}
	150
	151	Bool_t AliAODConversionMother::IsTrueMeson(AliStack *fMCStack,Int_t pdgcode){
	152	TParticle *part=GetMCParticle(fMCStack);
	153
	154	if(part){
	155	// Check if it is a true photon
	156	if(part->GetPdgCode()==pdgcode){
	157	return kTRUE;
	158	}
	159	}
	160	return kFALSE;
2eedd4ed	161	}
2eedd4ed	162
4803eb1f	163	Float_t AliAODConversionMother::CalculateDistanceBetweenPhotons(AliAODConversionPhoton* y1, AliAODConversionPhoton* y2 , Double_t prodPoint[3]){
	164
	165	TVector3 a(y1->GetConversionX(),y1->GetConversionY(),y1->GetConversionZ());
	166	TVector3 b(y1->GetPx(),y1->GetPy(),y1->GetPz());
	167	TVector3 c(y2->GetConversionX(),y2->GetConversionY(),y2->GetConversionZ());
	168	TVector3 d(y2->GetPx(),y2->GetPy(),y2->GetPz());
	169
	170	TVector3 n = b.Cross(d);
	171	TVector3 nn = n.Unit();
	172
	173	Double_t dist = 0;
	174	if (n.Mag() == 0){
	175	TVector3 e = a-c;
	176	if (d.Mag() != 0){
	177	dist = TMath::Abs((e.Cross(d)).Mag())/TMath::Abs(d.Mag());
	178	}
	179	prodPoint[0] = 0;
	180	prodPoint[1] = 0;
	181	prodPoint[2] = 0;
	182	} else {
	183	dist = TMath::Abs(n.Dot(c-a))/TMath::Abs(n.Mag());
	184	Double_t lambda = (b.Dot(d) * (a-c).Dot(d) - d.Dot(d) * (a-c).Dot(b))/(b.Dot(b) * d.Dot(d) - TMath::Power(b.Dot(d),2));
	185	Double_t mu = ((a-c).Dot(d) * b.Dot(b) - (a-c).Dot(b) * b.Dot(d) )/(b.Dot(b) * d.Dot(d) - TMath::Power(b.Dot(d),2));
	186
	187	TVector3 S1 = a + lambda* b;
	188	TVector3 S2 = c + mu* d;
	189	TVector3 Prod = S1 + 0.5dist(S2-S1).Unit();
	190	prodPoint[0] = Prod(0);
	191	prodPoint[1] = Prod(1);
	192	prodPoint[2] = Prod(2);
	193
	194	}
c15c23b8	195	if (dist > 1000) dist = 999.;
4803eb1f	196	return dist;
	197	}
	198
	199	///________________________________________________________________________
	200	void AliAODConversionMother::CalculateDistanceOfClossetApproachToPrimVtx(const AliVVertex* primVertex){
	201
	202	Double_t primCo[3] = {primVertex->GetX(),primVertex->GetY(),primVertex->GetZ()};
	203
	204	Double_t absoluteP = TMath::Sqrt(TMath::Power(Px(),2) + TMath::Power(Py(),2) + TMath::Power(Pz(),2));
	205	Double_t p[3] = {Px()/absoluteP,Py()/absoluteP,Pz()/absoluteP};
	206	Double_t CP[3];
	207
	208	CP[0] = fProductionVtx[0] - primCo[0];
	209	CP[1] = fProductionVtx[1] - primCo[1];
	210	CP[2] = fProductionVtx[2] - primCo[2];
	211
	212	Double_t Lambda = - (CP[0]p[0]+CP[1]p[1]+CP[2]p[2])/(p[0]p[0]+p[1]p[1]+p[2]p[2]);
	213
	214	Double_t S[3];
	215	S[0] = fProductionVtx[0] + p[0]*Lambda;
	216	S[1] = fProductionVtx[1] + p[1]*Lambda;
	217	S[2] = fProductionVtx[2] + p[2]*Lambda;
	218
	219	fdcaRPrimVtx = TMath::Sqrt( TMath::Power(primCo[0]-S[0],2) + TMath::Power(primCo[1]-S[1],2));
	220	fdcaZPrimVtx = primCo[2]-S[2];
	221
	222
	223	// cout << "DCA z: " << dca[1] << "\t DCA r: " << dca[0] << "\t DCA 3d: " << TMath::Sqrt(dca[1]dca[1] + dca[0]dca[0]) << endl;
	224
	225
	226	return;
	227	}
	228
	229	///________________________________________________________________________
	230	void AliAODConversionMother::DetermineMesonQuality(AliAODConversionPhoton* y1, AliAODConversionPhoton* y2){
	231	UChar_t photonQA1 = y1->GetPhotonQuality();
	232	UChar_t photonQA2 = y2->GetPhotonQuality();
	233
	234	if (photonQA1 == 0 \|\| photonQA2 == 0){
	235	fQuality = 0;
	236	return;
	237	}
	238	if (photonQA1 == 1 && photonQA2 == 1){
	239	fQuality = 1;
	240	return;
	241	}
	242	if (photonQA1 == 2 && photonQA2 == 2){
	243	fQuality = 4;
	244	return;
	245	}
	246	if (photonQA1 == 3 && photonQA2 == 3){
	247	fQuality = 6;
	248	return;
	249	}
	250	if (photonQA1 == 1){
	251	if (photonQA2 == 2){
	252	fQuality = 2;
	253	return;
	254	}
	255	if (photonQA2 == 3){
	256	fQuality = 3;
	257	return;
	258	}
	259	}
260	if (photonQA2 == 1){
261	if (photonQA1 == 2){
262	fQuality = 2;
263	return;
264	}
265	if (photonQA1 == 3){
266	fQuality = 3;
267	return;
268	}
269	}
270	if ((photonQA1 == 2 && photonQA2 == 3)\|\| (photonQA1 == 3 && photonQA2 == 2)){
271	fQuality = 5;
272	return;
273	}
274
275	}