RICH/AliRICHResponse.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 #include <TMath.h>
  18 #include <TParticle.h>
  19 #include <TRandom.h>
  20 #include <TVirtualMC.h>
  21
  22 #include "AliRICHResponse.h"
  23 #include "AliRun.h"
  24 #include "AliSegmentation.h"
  25 #include "AliMC.h"
  26
  27 ClassImp(AliRICHResponse)
  28 //__________________________________________________________________________________________________
  29 AliRICHResponse::AliRICHResponse()
  30 {
  31    SetSigmaIntegration(5.);
  32    SetChargeSlope(27.);
  33    SetChargeSpread(0.18, 0.18);
  34    SetMaxAdc(4096);
  35    SetAlphaFeedback(0.036);
  36    SetEIonisation(26.e-9);
  37    SetSqrtKx3(0.77459667);
  38    SetKx2(0.962);
  39    SetKx4(0.379);
  40    SetSqrtKy3(0.77459667);
  41    SetKy2(0.962);
  42    SetKy4(0.379);
  43    SetPitch(0.25);
  44    SetWireSag(1);                     // 1->On, 0->Off
  45    SetVoltage(2150);                  // Should only be 2000, 2050, 2100 or 2150
  46 }//AliRICHResponse::ctor()
  47 //__________________________________________________________________________________________________
  48 Float_t AliRICHResponse::IntPH(Float_t eloss, Float_t yhit)
  49 {
  50     // Get number of electrons and return charge
  51
  52     Int_t nel;
  53     nel= Int_t(eloss/fEIonisation);
  54
  55     Float_t charge=0;
  56     Double_t gain_var=1;
  57
  58     if (nel == 0) nel=1;
  59
  60     if (fWireSag)
  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           }
  66         if (fVoltage==2100)
  67             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;
  68         if (fVoltage==2050)
  69             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;
  70         if (fVoltage==2000)
  71             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;
  72
  73         gain_var = gain_var/100;
  74         //printf("Yhit:%f, Gain variation:%f\n",yhit,gain_var);
  75
  76         Float_t gain = (fChargeSlope + fChargeSlope*gain_var)*.9;
  77         //printf(" Yhit:%f, Gain variation:%f\n",yhit, gain);
  78
  79         for (Int_t i=1;i<=nel;i++) {
  80           charge -= gain*TMath::Log(gRandom->Rndm());
  81         }
  82       }
  83     else
  84       {
  85         for (Int_t i=1;i<=nel;i++) {
  86           charge -= fChargeSlope*TMath::Log(gRandom->Rndm());
  87         }
  88       }
  89
  90     return charge;
  91 }//InitPH()
  92 //__________________________________________________________________________________________________
  93 Float_t AliRICHResponse::IntPH(Float_t yhit)
  94 {//  Get number of electrons and return charge, for a single photon
  95
  96   Float_t charge=0;
  97   Double_t gain_var=1;
  98
  99    if (fWireSag)
 100       {
 101         if (fVoltage==2150)
 102           {
 103             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;
 104             //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;
 105           }
 106         if (fVoltage==2100)
 107             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;
 108         if (fVoltage==2050)
 109             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;
 110         if (fVoltage==2000)
 111             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;
 112
 113         gain_var = gain_var/100;
 114         //printf(" Yhit:%f, Gain variation:%f\n",yhit, gain_var);
 115
 116         Float_t gain = (fChargeSlope + fChargeSlope*gain_var)*.9;
 117
 118         charge -= gain*TMath::Log(gRandom->Rndm());
 119         //printf(" Yhit:%f, Gain variation:%f\n",yhit, gain);
 120       }
 121    else
 122      {
 123        charge -= fChargeSlope*TMath::Log(gRandom->Rndm());
 124      }
 125     return charge;
 126 }//IntPH()
 127 //__________________________________________________________________________________________________
 128 Float_t AliRICHResponse::IntXY(AliSegmentation * segmentation)
 129 {
 130
 131     const Float_t kInversePitch = 1/fPitch;
 132     Float_t response;
 133
 134     //  Integration limits defined by segmentation model
 135
 136     Float_t xi1, xi2, yi1, yi2;
 137     segmentation->IntegrationLimits(xi1,xi2,yi1,yi2);
 138
 139     xi1=xi1*kInversePitch;
 140     xi2=xi2*kInversePitch;
 141     yi1=yi1*kInversePitch;
 142     yi2=yi2*kInversePitch;
 143
 144     //printf("Integration Limits: %f-%f, %f-%f\n",xi1,xi2,yi1,yi2);
 145
 146     //printf("KInversePitch:%f\n",kInversePitch);
 147
 148     //
 149 // The Mathieson function
 150     Double_t ux1=fSqrtKx3*TMath::TanH(fKx2*xi1);
 151     Double_t ux2=fSqrtKx3*TMath::TanH(fKx2*xi2);
 152
 153     Double_t uy1=fSqrtKy3*TMath::TanH(fKy2*yi1);
 154     Double_t uy2=fSqrtKy3*TMath::TanH(fKy2*yi2);
 155
 156     //printf("Integration Data: %f-%f, %f-%f\n",ux1,ux2,uy1,uy2);
 157
 158     //printf("%f %f %f %f\n",fSqrtKx3,fKx2,fKy4,fKx4);
 159
 160     response=4.*fKx4*(TMath::ATan(ux2)-TMath::ATan(ux1))*fKy4*(TMath::ATan(uy2)-TMath::ATan(uy1));
 161
 162     //printf("Response:%f\n",response);
 163
 164     return response;
 165
 166 }//IntXY
 167 //__________________________________________________________________________________________________
 168 Int_t AliRICHResponse::FeedBackPhotons(Float_t *source, Float_t qtot)
 169 {  // Generate FeedBack photons
 170   Int_t j, ipart, nt;
 171
 172   Int_t sNfeed=0;
 173
 174
 175   // Local variables
 176   Double_t ranf[2];
 177   Float_t cthf, phif, enfp = 0, sthf;
 178   Int_t i, ifeed;
 179   Float_t e1[3], e2[3], e3[3];
 180   Float_t vmod, uswop;
 181   Float_t fp;
 182   Double_t random;
 183   Float_t dir[3], phi;
 184   Int_t nfp;
 185   Float_t pol[3], mom[4];
 186   TLorentzVector position;
 187   //
 188   // Determine number of feedback photons
 189
 190   //  Get weight of current particle
 191   TParticle *current = (TParticle*)
 192     (*gAlice->GetMCApp()->Particles())[gAlice->GetMCApp()->GetCurrentTrackNumber()];
 193
 194   ifeed = Int_t(current->GetWeight()/100+0.5);
 195   ipart = gMC->TrackPid();
 196   fp = fAlphaFeedback * qtot;
 197   nfp = gRandom->Poisson(fp);
 198
 199   // This call to fill the time of flight
 200   gMC->TrackPosition(position);
 201   //printf("Track position: %f %f %f %15.12f\n", position[0],position[1],position[2],position[3]);
 202   //
 203   // Generate photons
 204   for (i = 0; i <nfp; i++) {
 205
 206     // Direction
 207     gMC->GetRandom()->RndmArray(2,ranf);
 208     cthf = ranf[0] * 2 - 1.;
 209     if (cthf < 0)  continue;
 210     sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
 211     phif = ranf[1] * 2 * TMath::Pi();
 212     //
 213     //gMC->Rndm(&random, 1);
 214     gMC->GetRandom()->RndmArray(1, &random);
 215     if (random <= .57) {
 216       enfp = 7.5e-9;
 217     } else if (random <= .7) {
 218       enfp = 6.4e-9;
 219     } else {
 220       enfp = 7.9e-9;
 221     }
 222
 223     dir[0] = sthf * TMath::Sin(phif);
 224     dir[1] = cthf;
 225     dir[2] = sthf * TMath::Cos(phif);
 226     gMC->Gdtom(dir, mom, 2);
 227     mom[0]*=enfp;
 228     mom[1]*=enfp;
 229     mom[2]*=enfp;
 230     mom[3] = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2]);
 231     //printf("Dir %f %f %f\n",dir[0],dir[1],dir[2]);
 232     //printf("Momentum %15.12f %15.12f %15.12f\n",mom[0],mom[1],mom[2]);
 233     //printf("Energy %e\n", mom[3]);
 234
 235     // Polarisation
 236     e1[0] = 0;
 237     e1[1] = -dir[2];
 238     e1[2] = dir[1];
 239
 240     e2[0] = -dir[1];
 241     e2[1] = dir[0];
 242     e2[2] = 0;
 243
 244     e3[0] = dir[1];
 245     e3[1] = 0;
 246     e3[2] = -dir[0];
 247
 248     vmod=0;
 249     for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
 250     if (!vmod) for(j=0;j<3;j++) {
 251       uswop=e1[j];
 252       e1[j]=e3[j];
 253       e3[j]=uswop;
 254     }
 255     vmod=0;
 256     for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
 257     if (!vmod) for(j=0;j<3;j++) {
 258       uswop=e2[j];
 259       e2[j]=e3[j];
 260       e3[j]=uswop;
 261     }
 262
 263     vmod=0;
 264     for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
 265     vmod=TMath::Sqrt(1/vmod);
 266     for(j=0;j<3;j++) e1[j]*=vmod;
 267
 268     vmod=0;
 269     for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
 270     vmod=TMath::Sqrt(1/vmod);
 271     for(j=0;j<3;j++) e2[j]*=vmod;
 272
 273     //gMC->Rndm(ranf, 1);
 274     gMC->GetRandom()->RndmArray(1,ranf);
 275     phi = ranf[0] * 2 * TMath::Pi();
 276     for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi);
 277     gMC->Gdtom(pol, pol, 2);
 278
 279     // Put photon on the stack and label it as feedback (51, 52)
 280     ++sNfeed;
 281
 282     gAlice->GetMCApp()->PushTrack(Int_t(1), gAlice->GetMCApp()->GetCurrentTrackNumber(), Int_t(50000051),
 283                      mom[0],mom[1],mom[2],mom[3],source[0],source[1],source[2],position[3],pol[0],pol[1],pol[2],
 284                      kPFeedBackPhoton, nt, 1.);
 285
 286   }
 287   return(sNfeed);
 288 }//FeedBackPhotons()