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.5  2001/02/23 17:25:08  jbarbosa
  19   Re-definition of IntPH() to accomodate for wire sag effect.
  20
  21   Revision 1.4  2000/12/01 17:37:44  morsch
  22   Replace in argument of SetTrack(..) string constant by  kPFeedBackPhoton.
  23
  24   Revision 1.3  2000/10/03 21:44:09  morsch
  25   Use AliSegmentation and AliHit abstract base classes.
  26
  27   Revision 1.2  2000/10/02 21:28:12  fca
  28   Removal of useless dependecies via forward declarations
  29
  30   Revision 1.1  2000/06/12 15:29:37  jbarbosa
  31   Cleaned up version.
  32
  33 */
  34
  35 #include "AliRICHResponseV0.h"
  36 #include "AliSegmentation.h"
  37 #include "AliRun.h"
  38 #include "AliMC.h"
  39
  40 #include <TMath.h>
  41 #include <TRandom.h>
  42 #include <TParticle.h>
  43 //___________________________________________
  44 ClassImp(AliRICHResponseV0)
  45
  46 Float_t AliRICHResponseV0::IntPH(Float_t eloss, Float_t yhit)
  47 {
  48     // Get number of electrons and return charge
  49
  50     Int_t nel;
  51     nel= Int_t(eloss/fEIonisation);
  52
  53     Float_t charge=0;
  54     Double_t gain_var=1;
  55
  56     if (nel == 0) nel=1;
  57
  58     if (fWireSag)
  59       {
  60         //printf("Voltage:%d, Yhit:%f\n",fVoltage, yhit);
  61
  62         if (fVoltage==2150)
  63           {
  64             gain_var = 9e-6*TMath::Power(yhit,4) + 2e-7*TMath::Power(yhit,3) - 0.0316*TMath::Power(yhit,2) - 3e-4*yhit + 25.367;
  65             //gain_var = 9e-5*TMath::Power(yhit,4) + 2e-6*TMath::Power(yhit,3) - 0.316*TMath::Power(yhit,2) - 3e-3*yhit + 253.67;
  66           }
  67         if (fVoltage==2100)
  68             gain_var = 8e-6*TMath::Power(yhit,4) + 2e-7*TMath::Power(yhit,3) - 0.0283*TMath::Power(yhit,2) - 2e-4*yhit + 23.015;
  69         if (fVoltage==2050)
  70             gain_var = 7e-6*TMath::Power(yhit,4) + 1e-7*TMath::Power(yhit,3) - 0.0254*TMath::Power(yhit,2) - 2e-4*yhit + 20.888;
  71         if (fVoltage==2000)
  72             gain_var = 6e-6*TMath::Power(yhit,4) + 8e-8*TMath::Power(yhit,3) - 0.0227*TMath::Power(yhit,2) - 1e-4*yhit + 18.961;
  73
  74         gain_var = gain_var/100;
  75         //printf("Yhit:%f, Gain variation:%f\n",yhit,gain_var);
  76
  77         Float_t gain = (fChargeSlope + fChargeSlope*gain_var)*.9;
  78         //printf(" Yhit:%f, Gain variation:%f\n",yhit, gain);
  79
  80         for (Int_t i=1;i<=nel;i++) {
  81           charge -= gain*TMath::Log(gRandom->Rndm());
  82         }
  83       }
  84     else
  85       {
  86         for (Int_t i=1;i<=nel;i++) {
  87           charge -= fChargeSlope*TMath::Log(gRandom->Rndm());
  88         }
  89       }
  90
  91     return charge;
  92 }
  93
  94 Float_t AliRICHResponseV0::IntPH(Float_t yhit)
  95 {
  96
  97 //  Get number of electrons and return charge, for a single photon
  98
  99   Float_t charge=0;
 100   Double_t gain_var=1;
 101
 102    if (fWireSag)
 103       {
 104         if (fVoltage==2150)
 105           {
 106             gain_var = 9e-6*TMath::Power(yhit,4) + 2e-7*TMath::Power(yhit,3) - 0.0316*TMath::Power(yhit,2) - 3e-4*yhit + 25.367;
 107             //gain_var = 9e-5*TMath::Power(yhit,4) + 2e-6*TMath::Power(yhit,3) - 0.316*TMath::Power(yhit,2) - 3e-3*yhit + 253.67;
 108           }
 109         if (fVoltage==2100)
 110             gain_var = 8e-6*TMath::Power(yhit,4) + 2e-7*TMath::Power(yhit,3) - 0.0283*TMath::Power(yhit,2) - 2e-4*yhit + 23.015;
 111         if (fVoltage==2050)
 112             gain_var = 7e-6*TMath::Power(yhit,4) + 1e-7*TMath::Power(yhit,3) - 0.0254*TMath::Power(yhit,2) - 2e-4*yhit + 20.888;
 113         if (fVoltage==2000)
 114             gain_var = 6e-6*TMath::Power(yhit,4) + 8e-8*TMath::Power(yhit,3) - 0.0227*TMath::Power(yhit,2) - 1e-4*yhit + 18.961;
 115
 116         gain_var = gain_var/100;
 117         //printf(" Yhit:%f, Gain variation:%f\n",yhit, gain_var);
 118
 119         Float_t gain = (fChargeSlope + fChargeSlope*gain_var)*.9;
 120
 121         charge -= gain*TMath::Log(gRandom->Rndm());
 122         //printf(" Yhit:%f, Gain variation:%f\n",yhit, gain);
 123       }
 124    else
 125      {
 126        charge -= fChargeSlope*TMath::Log(gRandom->Rndm());
 127      }
 128     return charge;
 129 }
 130
 131
 132
 133 // -------------------------------------------
 134 Float_t AliRICHResponseV0::IntXY(AliSegmentation * segmentation)
 135 {
 136
 137     const Float_t kInversePitch = 1/fPitch;
 138     Float_t response;
 139 //
 140 //  Integration limits defined by segmentation model
 141 //
 142
 143     Float_t xi1, xi2, yi1, yi2;
 144     segmentation->IntegrationLimits(xi1,xi2,yi1,yi2);
 145
 146     xi1=xi1*kInversePitch;
 147     xi2=xi2*kInversePitch;
 148     yi1=yi1*kInversePitch;
 149     yi2=yi2*kInversePitch;
 150
 151     //printf("Integration Limits: %f-%f, %f-%f\n",xi1,xi2,yi1,yi2);
 152
 153     //printf("KInversePitch:%f\n",kInversePitch);
 154
 155     //
 156 // The Mathieson function
 157     Double_t ux1=fSqrtKx3*TMath::TanH(fKx2*xi1);
 158     Double_t ux2=fSqrtKx3*TMath::TanH(fKx2*xi2);
 159
 160     Double_t uy1=fSqrtKy3*TMath::TanH(fKy2*yi1);
 161     Double_t uy2=fSqrtKy3*TMath::TanH(fKy2*yi2);
 162
 163     //printf("Integration Data: %f-%f, %f-%f\n",ux1,ux2,uy1,uy2);
 164
 165     //printf("%f %f %f %f\n",fSqrtKx3,fKx2,fKy4,fKx4);
 166
 167     response=4.*fKx4*(TMath::ATan(ux2)-TMath::ATan(ux1))*fKy4*(TMath::ATan(uy2)-TMath::ATan(uy1));
 168
 169     //printf("Response:%f\n",response);
 170
 171     return response;
 172
 173 }
 174
 175 Int_t AliRICHResponseV0::FeedBackPhotons(Float_t *source, Float_t qtot)
 176 {
 177   //
 178   // Generate FeedBack photons
 179   //
 180   Int_t j, ipart, nt;
 181
 182   Int_t sNfeed=0;
 183
 184
 185   // Local variables
 186   Float_t cthf, ranf[2], phif, enfp = 0, sthf;
 187   Int_t i, ifeed;
 188   Float_t e1[3], e2[3], e3[3];
 189   Float_t vmod, uswop;
 190   Float_t fp, random;
 191   Float_t dir[3], phi;
 192   Int_t nfp;
 193   Float_t pol[3], mom[3];
 194   TLorentzVector position;
 195   //
 196   // Determine number of feedback photons
 197
 198   //  Get weight of current particle
 199   TParticle *current = (TParticle*)
 200     (*gAlice->Particles())[gAlice->CurrentTrack()];
 201
 202   ifeed = Int_t(current->GetWeight()/100+0.5);
 203   ipart = gMC->TrackPid();
 204   fp = fAlphaFeedback * qtot;
 205   nfp = gRandom->Poisson(fp);
 206
 207   // This call to fill the time of flight
 208   gMC->TrackPosition(position);
 209   //
 210   // Generate photons
 211   for (i = 0; i <nfp; i++) {
 212
 213     // Direction
 214     gMC->Rndm(ranf, 2);
 215     cthf = ranf[0] * 2 - 1.;
 216     if (cthf < 0)  continue;
 217     sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
 218     phif = ranf[1] * 2 * TMath::Pi();
 219     //
 220     gMC->Rndm(&random, 1);
 221     if (random <= .57) {
 222       enfp = 7.5e-9;
 223     } else if (random <= .7) {
 224       enfp = 6.4e-9;
 225     } else {
 226       enfp = 7.9e-9;
 227     }
 228
 229     dir[0] = sthf * TMath::Sin(phif);
 230     dir[1] = cthf;
 231     dir[2] = sthf * TMath::Cos(phif);
 232     gMC->Gdtom(dir, mom, 2);
 233     mom[0]*=enfp;
 234     mom[1]*=enfp;
 235     mom[2]*=enfp;
 236
 237     // Polarisation
 238     e1[0] = 0;
 239     e1[1] = -dir[2];
 240     e1[2] = dir[1];
 241
 242     e2[0] = -dir[1];
 243     e2[1] = dir[0];
 244     e2[2] = 0;
 245
 246     e3[0] = dir[1];
 247     e3[1] = 0;
 248     e3[2] = -dir[0];
 249
 250     vmod=0;
 251     for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
 252     if (!vmod) for(j=0;j<3;j++) {
 253       uswop=e1[j];
 254       e1[j]=e3[j];
 255       e3[j]=uswop;
 256     }
 257     vmod=0;
 258     for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
 259     if (!vmod) for(j=0;j<3;j++) {
 260       uswop=e2[j];
 261       e2[j]=e3[j];
 262       e3[j]=uswop;
 263     }
 264
 265     vmod=0;
 266     for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
 267     vmod=TMath::Sqrt(1/vmod);
 268     for(j=0;j<3;j++) e1[j]*=vmod;
 269
 270     vmod=0;
 271     for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
 272     vmod=TMath::Sqrt(1/vmod);
 273     for(j=0;j<3;j++) e2[j]*=vmod;
 274
 275     gMC->Rndm(ranf, 1);
 276     phi = ranf[0] * 2 * TMath::Pi();
 277     for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi);
 278     gMC->Gdtom(pol, pol, 2);
 279
 280     // Put photon on the stack and label it as feedback (51, 52)
 281     ++sNfeed;
 282
 283     gAlice->SetTrack(Int_t(1), gAlice->CurrentTrack(), Int_t(50000051),
 284                      mom,source,pol,position[3],
 285                      kPFeedBackPhoton, nt, 1.);
 286   }
 287   return(sNfeed);
 288 }
 289
 290
 291
 292