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