RICH/AliRICHv1.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 <TParticle.h>
  18 #include <TRandom.h>
  19 #include <TVirtualMC.h>
  20 #include <TPDGCode.h>
  21
  22 #include "AliRICHv1.h"
  23 #include "AliRICHParam.h"
  24 #include <AliConst.h>
  25 #include <AliPDG.h>
  26 #include <AliRun.h>
  27 #include <AliMC.h>
  28
  29 ClassImp(AliRICHv1)
  30 //______________________________________________________________________________
  31 AliRICHv1::AliRICHv1(const char *name, const char *title)
  32           :AliRICH(name,title)
  33 {// Full version of RICH with hits and diagnostics
  34   if(GetDebug())Info("named ctor","Start.");
  35   if(GetDebug())Info("named ctor","Stop.");
  36 }//named ctor
  37 //______________________________________________________________________________
  38 void AliRICHv1::Init()
  39 {//nothing to do here, all the work in ctor or CreateGeometry()
  40   if(GetDebug())Info("Init","Start.");
  41   if(GetDebug())Info("Init","Stop.");
  42 }//void AliRICHv1::Init()
  43 //______________________________________________________________________________
  44 void AliRICHv1::StepManager()
  45 {//Full Step Manager
  46
  47   Int_t          copy, id;
  48   static Int_t   iCurrentChamber;
  49   static Int_t   vol[2];
  50   Int_t          ipart;
  51   static Float_t hits[22];
  52   static Float_t ckovData[19];
  53   TLorentzVector x4,p4;
  54   Float_t        pos[3],mom[4],localPos[3],localMom[4];
  55   Float_t        theta,phi,localTheta,localPhi;
  56   Float_t        destep, step;
  57   Double_t        ranf[2];
  58   Int_t          nPads=kBad;
  59   Float_t        coscerenkov;
  60   static Float_t eloss,  tlength;
  61   const  Float_t kBig=1.e10;
  62
  63   TParticle *current = (TParticle*)(*gAlice->GetMCApp()->Particles())[gAlice->GetMCApp()->GetCurrentTrackNumber()];
  64
  65
  66
  67   id=gMC->CurrentVolID(copy); iCurrentChamber=copy;
  68   Float_t cherenkovLoss=0;
  69
  70   gMC->TrackPosition(x4);    pos[0]=x4(0);    pos[1]=x4(1);    pos[2]=x4(2);
  71   ckovData[1] = pos[0];   ckovData[2] = pos[1];   ckovData[3] = pos[2];
  72   ckovData[6] = 0;                      // dummy track length
  73
  74     /********************Store production parameters for Cerenkov photons************************/
  75
  76   if(gMC->TrackPid()==kCerenkov){//C
  77     Float_t ckovEnergy = current->Energy();
  78     if(ckovEnergy > 5.6e-09 && ckovEnergy < 7.8e-09 ){//C+E
  79       if(gMC->IsTrackEntering()){//C+E+enter
  80         if(gMC->VolId("FRE1")==gMC->CurrentVolID(copy)||gMC->VolId("FRE2")==gMC->CurrentVolID(copy)){//C+E+enter+FRE
  81           if(gMC->IsNewTrack()){//C+E+enter+FRE+new
  82             Int_t mother=current->GetFirstMother();
  83             ckovData[10]=mother;
  84             ckovData[11]=gAlice->GetMCApp()->GetCurrentTrackNumber();
  85             ckovData[12]=1;             //Media where photon was produced 1->Freon, 2->Quarz
  86             fCkovNumber++;
  87             fFreonProd=1;
  88           }//C+E+enter+FRE+new
  89         }//C+E+enter+FRE
  90         if(gMC->IsNewTrack()&&gMC->VolId("QUAR")==gMC->CurrentVolID(copy))  ckovData[12]=2;
  91       }//C+E+enter
  92       if(ckovData[12]==1){//C+E+produced in Freon
  93         if(gMC->IsTrackEntering()){                                     //is track entering?
  94                     if (gMC->VolId("META")==gMC->CurrentVolID(copy)){                //is it in gap?
  95                   gMC->TrackMomentum(p4);   mom[0]=p4(0);   mom[1]=p4(1);   mom[2]=p4(2);   mom[3]=p4(3);
  96
  97                         gMC->Gmtod(mom,localMom,2);
  98                         Float_t cophi = TMath::Cos(TMath::ATan2(localMom[0], localMom[1]));
  99                         Float_t t = (1. - .025 / cophi) * (1. - .05 /  cophi);
 100                         /**************** Photons lost in second grid have to be calculated by hand************/
 101                         gMC->GetRandom()->RndmArray(1,ranf);
 102                         if (ranf[0] > t) {
 103                           gMC->StopTrack();
 104                           ckovData[13] = 5;
 105                           //AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
 106                         }
 107                         /**********************************************************************************/
 108                       }    //gap
 109
 110                     if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy)){             //is it in csi?
 111
 112                         gMC->TrackMomentum(p4); mom[0]=p4(0);   mom[1]=p4(1);   mom[2]=p4(2);   mom[3]=p4(3);
 113
 114                         gMC->Gmtod(mom,localMom,2);
 115                         /********* Photons lost by Fresnel reflection have to be calculated by hand********/
 116                         /***********************Cerenkov phtons (always polarised)*************************/
 117                         Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
 118                         Double_t localRt = TMath::Sqrt(localTc);
 119                         localTheta   = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])));
 120                         Double_t cotheta = TMath::Abs(cos(localTheta));
 121                         Float_t t = Fresnel(ckovEnergy*1e9,cotheta,1);
 122                             gMC->GetRandom()->RndmArray(1,ranf);
 123                             if (ranf[0] < t) {
 124                               gMC->StopTrack();
 125                               ckovData[13] = 6;
 126                               //AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
 127
 128                               //printf("Added One (2)!\n");
 129                               //printf("Lost by Fresnel\n");
 130                             }
 131                             /**********************************************************************************/
 132                       }
 133                   } //track entering?
 134                   /********************Evaluation of losses************************/
 135                   /******************still in the old fashion**********************/
 136                   TArrayI procs;
 137                   Int_t i1 = gMC->StepProcesses(procs);            //number of physics mechanisms acting on the particle
 138                   for (Int_t i = 0; i < i1; ++i) {
 139                     //        Reflection loss
 140                     if (procs[i] == kPLightReflection) {        //was it reflected
 141                       ckovData[13]=10;
 142                       if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
 143                         ckovData[13]=1;
 144                       if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
 145                         ckovData[13]=2;
 146                     } //reflection question
 147
 148                     //        Absorption loss
 149                     else if (procs[i] == kPLightAbsorption) {              //was it absorbed?
 150                       //printf("Got in absorption\n");
 151                       ckovData[13]=20;
 152                       if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
 153                         ckovData[13]=11;
 154                       if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
 155                         ckovData[13]=12;
 156                       if (gMC->CurrentVolID(copy) == gMC->VolId("META"))
 157                         ckovData[13]=13;
 158                       if (gMC->CurrentVolID(copy) == gMC->VolId("GAP "))
 159                         ckovData[13]=13;
 160
 161                       if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC"))
 162                         ckovData[13]=15;
 163
 164                       //        CsI inefficiency
 165                       if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
 166                         ckovData[13]=16;
 167                       }
 168                       gMC->StopTrack();
 169                       //AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
 170                     } //absorption question
 171
 172
 173                     //        Photon goes out of tracking scope
 174                     else if (procs[i] == kPStop) {                 //is it below energy treshold?
 175                       ckovData[13]=21;
 176                       gMC->StopTrack();
 177                       //AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
 178                     }   // energy treshold question
 179                   }  //number of mechanisms cycle
 180                   /**********************End of evaluation************************/
 181       }//C+E+produced in Freon
 182     }//C+E
 183   }//C
 184 //*************************************End of Production Parameters Storing*********************
 185
 186   if(gMC->TrackPid()==kCerenkov||gMC->TrackPid()==kFeedback){//CF
 187     if(gMC->VolId("CSI ")==gMC->CurrentVolID(copy)){//CF+CSI
 188       if(gMC->Edep()>0.){//CF+CSI+DE
 189         gMC->TrackPosition(x4);   pos[0]=x4(0);   pos[1]=x4(1);   pos[2]=x4(2);
 190         gMC->TrackMomentum(p4);   mom[0]=p4(0);   mom[1]=p4(1);   mom[2]=p4(2);   mom[3]=p4(3);
 191
 192         Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
 193         Double_t rt = TMath::Sqrt(tc);
 194         theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
 195         phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
 196
 197         gMC->CurrentVolOffID(2,copy);   vol[0]=copy;    iCurrentChamber=vol[0];
 198
 199         gMC->Gmtod(pos,localPos,1);     gMC->Gmtod(mom,localMom,2);
 200
 201         Param()->SigGenInit(localPos[0],localPos[2]);
 202         ckovData[0]=gMC->TrackPid();
 203         ckovData[1]=pos[0];      ckovData[2]=pos[1];    ckovData[3]=pos[2];
 204         ckovData[4]=theta;       ckovData[5]=phi;       //theta-phi angles of incidence
 205         ckovData[8]=(Float_t) fNsdigits;      // first sdigit
 206         ckovData[9]=-1;                       // last pad hit
 207         ckovData[13]=4;                       // photon was detected
 208         ckovData[14]=mom[0];     ckovData[15]=mom[1];    ckovData[16]=mom[2];
 209
 210         destep = gMC->Edep();
 211         gMC->SetMaxStep(kBig);
 212         cherenkovLoss  += destep;
 213         ckovData[7]=cherenkovLoss;
 214
 215         GenerateFeedbacks(iCurrentChamber,cherenkovLoss);//CF+CSI+DE
 216
 217         if (fNsdigits > (Int_t)ckovData[8]) {
 218           ckovData[8]= ckovData[8]+1;
 219           ckovData[9]= (Float_t) fNsdigits;
 220         }
 221
 222         ckovData[17] = nPads;
 223         AliRICHhit *mipHit =  (AliRICHhit*) (fHits->UncheckedAt(0));
 224         if(mipHit){
 225           mom[0] = current->Px();   mom[1] = current->Py();   mom[2] = current->Pz();
 226           Float_t mipPx = mipHit->MomX();   Float_t mipPy = mipHit->MomY();   Float_t mipPz = mipHit->MomZ();
 227
 228           Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
 229           Float_t rt = TMath::Sqrt(r);
 230           Float_t mipR = mipPx*mipPx + mipPy*mipPy + mipPz*mipPz;
 231           Float_t mipRt = TMath::Sqrt(mipR);
 232           if((rt*mipRt) > 0)
 233             coscerenkov = (mom[0]*mipPx + mom[1]*mipPy + mom[2]*mipPz)/(rt*mipRt);
 234           else
 235             coscerenkov = 0;
 236           ckovData[18]=TMath::ACos(coscerenkov);//Cerenkov angle
 237         }
 238         AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),vol,ckovData);//PHOTON HIT CF+CSI+DE
 239         //AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
 240       }//CF+CSI+DE
 241     }//CF+CSI
 242   }/*CF*/else if(gMC->TrackCharge()){//MIP
 243     if(gMC->VolId("FRE1")==gMC->CurrentVolID(copy)||gMC->VolId("FRE2")==gMC->CurrentVolID(copy)){//MIP+FRE
 244       gMC->TrackMomentum(p4);       mom[0]=p4(0);   mom[1]=p4(1);   mom[2]=p4(2);   mom[3]=p4(3);
 245       hits[19]=mom[0];      hits [20] = mom[1];     hits [21] = mom[2];     fFreonProd=1;
 246     }//MIP+FRE
 247     if(gMC->VolId("GAP ")==gMC->CurrentVolID(copy)){//MIP+GAP
 248       gMC->CurrentVolOffID(3,copy);    vol[0]=copy;    iCurrentChamber=vol[0];
 249       gMC->TrackPosition(x4);   pos[0]=x4(0);   pos[1]=x4(1);   pos[2]=x4(2);
 250       gMC->TrackMomentum(p4);   mom[0]=p4(0);   mom[1]=p4(1);   mom[2]=p4(2);    mom[3]=p4(3);
 251       gMC->Gmtod(pos,localPos,1);           gMC->Gmtod(mom,localMom,2);
 252       ipart =gMC->TrackPid();
 253       destep = gMC->Edep();step   = gMC->TrackStep();// momentum loss and steplength in last step
 254       if(gMC->IsTrackEntering()){//MIP+GAP+Enter  record hit when mip enters ...
 255         Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
 256         Double_t rt = TMath::Sqrt(tc);
 257         theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
 258         phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
 259         Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
 260         Double_t localRt = TMath::Sqrt(localTc);
 261         localTheta   = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])))*kRaddeg;
 262         localPhi     = Float_t(TMath::ATan2(Double_t(localMom[2]),Double_t(localMom[0])))*kRaddeg;
 263         hits[0] = Float_t(ipart);         // particle type
 264         hits[1] = localPos[0];    hits[2] = localPos[1];     hits[3] = localPos[2];
 265         hits[4] = localTheta;     hits[5] = localPhi;              // theta-phi angles of incidence
 266         hits[8] = (Float_t) fNsdigits;    // first sdigit
 267         hits[9] = -1;                     // last pad hit
 268         hits[13] = fFreonProd;           // did id hit the freon?
 269         hits[14] = mom[0];        hits[15] = mom[1];        hits[16] = mom[2];
 270         hits[18] = 0;               // dummy cerenkov angle
 271         tlength = 0;        eloss   = 0;        fFreonProd = 0;
 272         C(iCurrentChamber)->LocaltoGlobal(localPos,hits+1);
 273
 274         Param()->SigGenInit(localPos[0], localPos[2]);
 275       }/*MIP+GAP+Enter*/else if(gMC->IsTrackExiting()||gMC->IsTrackStop()||gMC->IsTrackDisappeared()){//MIP+GAP+Exit
 276         gMC->SetMaxStep(kBig);
 277         eloss   += destep;
 278         tlength += step;
 279         if (eloss > 0) {
 280           if(gMC->TrackPid() == kNeutron) printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
 281           GenerateFeedbacks(iCurrentChamber,eloss);//MIP+GAP+Exit
 282           hits[17] = nPads;
 283         }
 284         hits[6]=tlength;        hits[7]=eloss;
 285         if(fNsdigits > (Int_t)hits[8]) {
 286           hits[8]= hits[8]+1;
 287           hits[9]= (Float_t) fNsdigits;
 288         }
 289         AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),vol,hits);//MIP HIT MIP+GAP+Exit
 290         eloss = 0;
 291       }/*MIP+GAP+Exit*/else if(Param()->SigGenCond(localPos[0], localPos[2])){//MIP+GAP+Spec
 292         Param()->SigGenInit(localPos[0], localPos[2]);
 293         if(eloss>0){
 294           if(gMC->TrackPid() == kNeutron) printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
 295           GenerateFeedbacks(iCurrentChamber,eloss);//MIP+GAP+Spec
 296           hits[17] = nPads;
 297         }
 298
 299         eloss    = destep;      tlength += step ;
 300       }/*MIP+GAP+Spec*/else{//MIP+GAP+nothing special
 301         eloss   += destep;
 302         tlength += step ;
 303       }//MIP+GAP+nothing special
 304     }//MIP+GAP
 305   }//MIP
 306 }//void AliRICHv1::StepManager()