PWGGA/GammaConv/AliAODConversionMother.cxx

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