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