RICH/AliRICHResponseV0.cxx

   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
  16 /*
  17   $Log$
  18 */
  19
  20 #include "AliRICHResponseV0.h"
  21 #include "AliRICHSegmentation.h"
  22 #include "AliRun.h"
  23
  24 #include <TMath.h>
  25 #include <TRandom.h>
  26 #include <TParticle.h>
  27 //___________________________________________
  28 ClassImp(AliRICHResponseV0)
  29
  30 Float_t AliRICHResponseV0::IntPH(Float_t eloss)
  31 {
  32     // Get number of electrons and return charge
  33
  34     Int_t nel;
  35     nel= Int_t(eloss/fEIonisation);
  36
  37     Float_t charge=0;
  38     if (nel == 0) nel=1;
  39     for (Int_t i=1;i<=nel;i++) {
  40         charge -= fChargeSlope*TMath::Log(gRandom->Rndm());
  41     }
  42     return charge;
  43 }
  44
  45 Float_t AliRICHResponseV0::IntPH()
  46 {
  47
  48 //  Get number of electrons and return charge, for a single photon
  49
  50     Float_t charge = -fChargeSlope*TMath::Log(gRandom->Rndm());
  51     return charge;
  52 }
  53
  54
  55
  56 // -------------------------------------------
  57 Float_t AliRICHResponseV0::IntXY(AliRICHSegmentation * segmentation)
  58 {
  59
  60     const Float_t kInversePitch = 1/fPitch;
  61     Float_t response;
  62 //
  63 //  Integration limits defined by segmentation model
  64 //
  65
  66     Float_t xi1, xi2, yi1, yi2;
  67     segmentation->IntegrationLimits(xi1,xi2,yi1,yi2);
  68
  69     xi1=xi1*kInversePitch;
  70     xi2=xi2*kInversePitch;
  71     yi1=yi1*kInversePitch;
  72     yi2=yi2*kInversePitch;
  73
  74     //printf("Integration Limits: %f-%f, %f-%f\n",xi1,xi2,yi1,yi2);
  75
  76     //printf("KInversePitch:%f\n",kInversePitch);
  77
  78     //
  79 // The Mathieson function
  80     Double_t ux1=fSqrtKx3*TMath::TanH(fKx2*xi1);
  81     Double_t ux2=fSqrtKx3*TMath::TanH(fKx2*xi2);
  82
  83     Double_t uy1=fSqrtKy3*TMath::TanH(fKy2*yi1);
  84     Double_t uy2=fSqrtKy3*TMath::TanH(fKy2*yi2);
  85
  86     //printf("Integration Data: %f-%f, %f-%f\n",ux1,ux2,uy1,uy2);
  87
  88     //printf("%f %f %f %f\n",fSqrtKx3,fKx2,fKy4,fKx4);
  89
  90     response=4.*fKx4*(TMath::ATan(ux2)-TMath::ATan(ux1))*fKy4*(TMath::ATan(uy2)-TMath::ATan(uy1));
  91
  92     //printf("Response:%f\n",response);
  93
  94     return response;
  95
  96 }
  97
  98 Int_t AliRICHResponseV0::FeedBackPhotons(Float_t *source, Float_t qtot)
  99 {
 100   //
 101   // Generate FeedBack photons
 102   //
 103   Int_t j, ipart, nt;
 104
 105   Int_t sNfeed=0;
 106
 107
 108   // Local variables
 109   Float_t cthf, ranf[2], phif, enfp = 0, sthf;
 110   Int_t i, ifeed;
 111   Float_t e1[3], e2[3], e3[3];
 112   Float_t vmod, uswop;
 113   Float_t fp, random;
 114   Float_t dir[3], phi;
 115   Int_t nfp;
 116   Float_t pol[3], mom[3];
 117   TLorentzVector position;
 118   //
 119   // Determine number of feedback photons
 120
 121   //  Get weight of current particle
 122   TParticle *current = (TParticle*)
 123     (*gAlice->Particles())[gAlice->CurrentTrack()];
 124
 125   ifeed = Int_t(current->GetWeight()/100+0.5);
 126   ipart = gMC->TrackPid();
 127   fp = fAlphaFeedback * qtot;
 128   nfp = gRandom->Poisson(fp);
 129
 130   // This call to fill the time of flight
 131   gMC->TrackPosition(position);
 132   //
 133   // Generate photons
 134   for (i = 0; i <nfp; i++) {
 135
 136     // Direction
 137     gMC->Rndm(ranf, 2);
 138     cthf = ranf[0] * 2 - 1.;
 139     if (cthf < 0)  continue;
 140     sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
 141     phif = ranf[1] * 2 * TMath::Pi();
 142     //
 143     gMC->Rndm(&random, 1);
 144     if (random <= .57) {
 145       enfp = 7.5e-9;
 146     } else if (random <= .7) {
 147       enfp = 6.4e-9;
 148     } else {
 149       enfp = 7.9e-9;
 150     }
 151
 152     dir[0] = sthf * TMath::Sin(phif);
 153     dir[1] = cthf;
 154     dir[2] = sthf * TMath::Cos(phif);
 155     gMC->Gdtom(dir, mom, 2);
 156     mom[0]*=enfp;
 157     mom[1]*=enfp;
 158     mom[2]*=enfp;
 159
 160     // Polarisation
 161     e1[0] = 0;
 162     e1[1] = -dir[2];
 163     e1[2] = dir[1];
 164
 165     e2[0] = -dir[1];
 166     e2[1] = dir[0];
 167     e2[2] = 0;
 168
 169     e3[0] = dir[1];
 170     e3[1] = 0;
 171     e3[2] = -dir[0];
 172
 173     vmod=0;
 174     for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
 175     if (!vmod) for(j=0;j<3;j++) {
 176       uswop=e1[j];
 177       e1[j]=e3[j];
 178       e3[j]=uswop;
 179     }
 180     vmod=0;
 181     for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
 182     if (!vmod) for(j=0;j<3;j++) {
 183       uswop=e2[j];
 184       e2[j]=e3[j];
 185       e3[j]=uswop;
 186     }
 187
 188     vmod=0;
 189     for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
 190     vmod=TMath::Sqrt(1/vmod);
 191     for(j=0;j<3;j++) e1[j]*=vmod;
 192
 193     vmod=0;
 194     for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
 195     vmod=TMath::Sqrt(1/vmod);
 196     for(j=0;j<3;j++) e2[j]*=vmod;
 197
 198     gMC->Rndm(ranf, 1);
 199     phi = ranf[0] * 2 * TMath::Pi();
 200     for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi);
 201     gMC->Gdtom(pol, pol, 2);
 202
 203     // Put photon on the stack and label it as feedback (51, 52)
 204     ++sNfeed;
 205
 206     gAlice->SetTrack(Int_t(1), gAlice->CurrentTrack(), Int_t(50000051),
 207                      mom,source,pol,position[3],
 208                      "Feedback", nt, 1.);
 209   }
 210   return(sNfeed);
 211 }
 212
 213
 214
 215