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