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