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 "AliRICHv1.h"     //class header
  18 #include "AliRICHParam.h"
  19 #include "AliRICHChamber.h"
  20 #include <TParticle.h>
  21 #include <TRandom.h>
  22 #include <TVirtualMC.h>
  23 #include <TPDGCode.h>
  24 #include <AliStack.h>      //Hits2SDigits()
  25 #include <AliConst.h>
  26 #include <AliPDG.h>
  27 #include <AliRun.h>
  28 #include <AliMC.h>
  29 #include <AliRawDataHeader.h> //Digits2Raw()
  30
  31 ClassImp(AliRICHv1)
  32 //__________________________________________________________________________________________________
  33 void AliRICHv1::StepManager()
  34 {
  35 // Full Step Manager.
  36 // 3- Ckovs absorbed on Collection electrods
  37 // 5- Ckovs absorbed on Cathode wires
  38 // 6- Ckovs absorbed on Anod wires
  39
  40   Int_t          copy;
  41   static Int_t   iCurrentChamber;
  42   static Int_t idRRAD = gMC->VolId("RRAD");
  43   static Int_t idRRWI = gMC->VolId("RRWI");
  44   static Int_t idRICH = gMC->VolId("RICH");
  45   static Int_t idRPC  = gMC->VolId("RPC ");
  46   static Int_t idRGAP = gMC->VolId("RGAP");
  47 //history of Cerenkovs
  48   if(gMC->TrackPid()==kCerenkov){
  49     if( gMC->IsNewTrack()   && gMC->CurrentVolID(copy)==idRRAD) fCounters(0)++;// 0- Ckovs produced in radiator
  50     if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==idRRAD) fCounters(1)++;// 1- Ckovs absorbed in radiator
  51     if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==idRRWI) fCounters(2)++;// 2- Ckovs absorbed in radiator window
  52     if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==idRICH) fCounters(4)++;// 4- Ckovs absorbed in CH4
  53   }
  54
  55 //Treat photons
  56   static TLorentzVector cerX4;
  57   if((gMC->TrackPid()==kCerenkov||gMC->TrackPid()==kFeedback)&&gMC->CurrentVolID(copy)==idRPC){//photon in PC
  58     if(gMC->Edep()>0){//photon in PC +DE
  59       if(IsLostByFresnel()){
  60         if(gMC->TrackPid()==kCerenkov) fCounters(7)++;// 7- Ckovs reflected from CsI
  61         gMC->StopTrack();
  62         return;
  63       }
  64       gMC->TrackPosition(cerX4); gMC->CurrentVolOffID(2,iCurrentChamber);//RICH-RPPF-RPC
  65
  66       HitAdd(iCurrentChamber,gAlice->GetMCApp()->GetCurrentTrackNumber(),cerX4.Vect(),cerX4.Vect());//HIT for PHOTON in conditions CF+CSI+DE
  67       fCounters(8)++;//4- Ckovs converted to electron on CsI
  68       GenerateFeedbacks(iCurrentChamber);
  69     }//photon in PC and DE >0
  70   }//photon in PC
  71
  72 //Treat charged particles
  73   static Float_t eloss;
  74   static TLorentzVector mipInX4,mipOutX4;
  75   if(gMC->TrackCharge() && gMC->CurrentVolID(copy)==idRGAP){//MIP in GAP
  76     gMC->CurrentVolOffID(1,iCurrentChamber);//RICH-RGAP
  77     if(gMC->IsTrackEntering()||gMC->IsNewTrack()) {//MIP in GAP entering or newly created
  78       eloss=0;
  79       gMC->TrackPosition(mipInX4);
  80     }else if(gMC->IsTrackExiting()||gMC->IsTrackStop()||gMC->IsTrackDisappeared()){//MIP in GAP exiting or disappeared
  81       eloss+=gMC->Edep();//take into account last step dEdX
  82       gMC->TrackPosition(mipOutX4);
  83       HitAdd(iCurrentChamber,gAlice->GetMCApp()->GetCurrentTrackNumber(),mipInX4.Vect(),mipOutX4.Vect(),eloss);//HIT for MIP: MIP in GAP Exiting
  84       GenerateFeedbacks(iCurrentChamber,eloss);//MIP+GAP+Exit
  85     }else//MIP in GAP going inside
  86       eloss   += gMC->Edep();
  87   }//MIP in GAP
  88 }//StepManager()
  89 //__________________________________________________________________________________________________
  90 Bool_t AliRICHv1::IsLostByFresnel()
  91 {
  92 // Calculate probability for the photon to be lost by Fresnel reflection.
  93   TLorentzVector p4;
  94   Double_t mom[3],localMom[3];
  95   gMC->TrackMomentum(p4);   mom[0]=p4(1);   mom[1]=p4(2);   mom[2]=p4(3);
  96   localMom[0]=0; localMom[1]=0; localMom[2]=0;
  97   gMC->Gmtod(mom,localMom,2);
  98   Double_t localTc    = localMom[0]*localMom[0]+localMom[2]*localMom[2];
  99   Double_t localTheta = TMath::ATan2(TMath::Sqrt(localTc),localMom[1]);
 100   Double_t cotheta = TMath::Abs(TMath::Cos(localTheta));
 101   if(gMC->GetRandom()->Rndm() < Fresnel(p4.E()*1e9,cotheta,1)){
 102     AliDebug(1,"Photon lost");
 103     return kTRUE;
 104   }else
 105     return kFALSE;
 106 }//IsLostByFresnel()
 107 //__________________________________________________________________________________________________
 108 void AliRICHv1::GenerateFeedbacks(Int_t iChamber,Float_t eloss)
 109 {
 110 // Generate FeedBack photons for the current particle. To be invoked from StepManager().
 111 // eloss=0 means photon so only pulse height distribution is to be analysed. This one is done in AliRICHParam::TotQdc()
 112
 113   TLorentzVector x4;
 114   gMC->TrackPosition(x4);
 115   TVector2 x2=C(iChamber)->Mrs2Pc(x4);//hit position on photocathode plane
 116   TVector2 xspe=x2;
 117   Int_t sector=P()->Loc2Sec(xspe);  if(sector==-1) return; //hit in dead zone, nothing to produce
 118   Int_t iTotQdc=P()->TotQdc(x2,eloss);
 119   Int_t iNphotons=gMC->GetRandom()->Poisson(P()->C(iChamber)->AlphaFeedback(sector)*iTotQdc);
 120   AliDebug(1,Form("N photons=%i",iNphotons));
 121   Int_t j;
 122   Float_t cthf, phif, enfp = 0, sthf, e1[3], e2[3], e3[3], vmod, uswop,dir[3], phi,pol[3], mom[4];
 123 //Generate photons
 124   for(Int_t i=0;i<iNphotons;i++){//feedbacks loop
 125     Double_t ranf[2];
 126     gMC->GetRandom()->RndmArray(2,ranf);    //Sample direction
 127     cthf=ranf[0]*2-1.0;
 128     if(cthf<0) continue;
 129     sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
 130     phif = ranf[1] * 2 * TMath::Pi();
 131
 132     if(Double_t randomNumber=gMC->GetRandom()->Rndm()<=0.57)
 133       enfp = 7.5e-9;
 134     else if(randomNumber<=0.7)
 135       enfp = 6.4e-9;
 136     else
 137       enfp = 7.9e-9;
 138
 139
 140     dir[0] = sthf * TMath::Sin(phif);    dir[1] = cthf;    dir[2] = sthf * TMath::Cos(phif);
 141     gMC->Gdtom(dir, mom, 2);
 142     mom[0]*=enfp;    mom[1]*=enfp;    mom[2]*=enfp;
 143     mom[3] = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2]);
 144
 145     // Polarisation
 146     e1[0]=      0;    e1[1]=-dir[2];    e1[2]= dir[1];
 147     e2[0]=-dir[1];    e2[1]= dir[0];    e2[2]=      0;
 148     e3[0]= dir[1];    e3[1]=      0;    e3[2]=-dir[0];
 149
 150     vmod=0;
 151     for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
 152     if (!vmod) for(j=0;j<3;j++) {
 153       uswop=e1[j];
 154       e1[j]=e3[j];
 155       e3[j]=uswop;
 156     }
 157     vmod=0;
 158     for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
 159     if (!vmod) for(j=0;j<3;j++) {
 160       uswop=e2[j];
 161       e2[j]=e3[j];
 162       e3[j]=uswop;
 163     }
 164
 165     vmod=0;  for(j=0;j<3;j++) vmod+=e1[j]*e1[j];  vmod=TMath::Sqrt(1/vmod);  for(j=0;j<3;j++) e1[j]*=vmod;
 166     vmod=0;  for(j=0;j<3;j++) vmod+=e2[j]*e2[j];  vmod=TMath::Sqrt(1/vmod);  for(j=0;j<3;j++) e2[j]*=vmod;
 167
 168     phi = gMC->GetRandom()->Rndm()* 2 * TMath::Pi();
 169     for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi);
 170     gMC->Gdtom(pol, pol, 2);
 171     Int_t outputNtracksStored;
 172     gAlice->GetMCApp()->PushTrack(1,                 //do not transport
 173                      gAlice->GetMCApp()->GetCurrentTrackNumber(),//parent track
 174                      kFeedback,                      //PID
 175                      mom[0],mom[1],mom[2],mom[3],    //track momentum
 176                      x4.X(),x4.Y(),x4.Z(),x4.T(),    //track origin
 177                      pol[0],pol[1],pol[2],           //polarization
 178                      kPFeedBackPhoton,
 179                      outputNtracksStored,
 180                      1.0);
 181   }//feedbacks loop
 182   AliDebug(1,"Stop.");
 183 }//GenerateFeedbacks()
 184 //__________________________________________________________________________________________________
 185 void AliRICHv1::Hits2SDigits()
 186 {
 187 // Create a list of sdigits corresponding to list of hits. Every hit generates one or more sdigits.
 188   AliDebug(1,"Start.");
 189   for(Int_t iEventN=0;iEventN<GetLoader()->GetRunLoader()->GetAliRun()->GetEventsPerRun();iEventN++){//events loop
 190     GetLoader()->GetRunLoader()->GetEvent(iEventN);//get next event
 191
 192     if(!GetLoader()->TreeH()) GetLoader()->LoadHits();    GetLoader()->GetRunLoader()->LoadHeader();
 193     if(!GetLoader()->GetRunLoader()->TreeK())             GetLoader()->GetRunLoader()->LoadKinematics();//from
 194     if(!GetLoader()->TreeS()) GetLoader()->MakeTree("S"); MakeBranch("S");//to
 195
 196     for(Int_t iPrimN=0;iPrimN<GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop
 197       GetLoader()->TreeH()->GetEntry(iPrimN);
 198       for(Int_t iHitN=0;iHitN<Hits()->GetEntries();iHitN++){//hits loop
 199         AliRICHHit *pHit=(AliRICHHit*)Hits()->At(iHitN);//get current hit
 200         TVector2 x2 = C(pHit->C())->Mrs2Anod(0.5*(pHit->InX3()+pHit->OutX3()));//hit position in the anod plane
 201         Int_t iTotQdc=P()->TotQdc(x2,pHit->Eloss());//total charge produced by hit, 0 if hit in dead zone
 202         if(iTotQdc==0) continue;
 203         //
 204         //need to quantize the anod....
 205         TVector padHit=AliRICHParam::Loc2Pad(x2);
 206         TVector2 padHitXY=AliRICHParam::Pad2Loc(padHit);
 207         TVector2 anod;
 208         if((x2.Y()-padHitXY.Y())>0) anod.Set(x2.X(),padHitXY.Y()+AliRICHParam::PitchAnod()/2);
 209         else anod.Set(x2.X(),padHitXY.Y()-AliRICHParam::PitchAnod()/2);
 210         //end to quantize anod
 211         //
 212         TVector area=P()->Loc2Area(anod);//determine affected pads, dead zones analysed inside
 213         AliDebug(1,Form("hitanod(%6.2f,%6.2f)->area(%3.0f,%3.0f)-(%3.0f,%3.0f) QDC=%4i",anod.X(),anod.Y(),area[0],area[1],area[2],area[3],iTotQdc));
 214         TVector pad(2);
 215         for(pad[1]=area[1];pad[1]<=area[3];pad[1]++)//affected pads loop
 216           for(pad[0]=area[0];pad[0]<=area[2];pad[0]++){
 217             Double_t padQdc=iTotQdc*P()->FracQdc(anod,pad);
 218             AliDebug(1,Form("current pad(%3.0f,%3.0f) with QDC  =%6.2f",pad[0],pad[1],padQdc));
 219             if(padQdc>0.1) SDigitAdd(pHit->C(),pad,padQdc,GetLoader()->GetRunLoader()->Stack()->Particle(pHit->GetTrack())->GetPdgCode(),pHit->GetTrack());
 220           }//affected pads loop
 221       }//hits loop
 222     }//prims loop
 223     GetLoader()->TreeS()->Fill();
 224     GetLoader()->WriteSDigits("OVERWRITE");
 225     SDigitsReset();
 226   }//events loop
 227   GetLoader()->UnloadHits(); GetLoader()->GetRunLoader()->UnloadHeader(); GetLoader()->GetRunLoader()->UnloadKinematics();
 228   GetLoader()->UnloadSDigits();
 229   AliDebug(1,"Stop.");
 230 }//Hits2SDigits()
 231 //__________________________________________________________________________________________________
 232 void AliRICHv1::Digits2Raw()
 233 {
 234 //Creates raw data files in DDL format. Invoked by AliSimulation
 235 //loop over events is done outside in AliSimulation
 236 //Arguments: none
 237 //  Returns: none
 238   AliDebug(1,"Start.");
 239   GetLoader()->LoadDigits();
 240   GetLoader()->TreeD()->GetEntry(0);
 241
 242   ofstream file[AliRICHDigit::kNddls];   //output streams
 243   Int_t    cnt[AliRICHDigit::kNddls];        //data words counters for DDLs
 244   AliRawDataHeader header; //empty DDL header
 245   UInt_t w32=0;            //32 bits data word
 246
 247   for(Int_t i=0;i<AliRICHDigit::kNddls;i++){//open all 14 DDL in parallel
 248     file[i].open(Form("RICH_%4i.ddl",AliRICHDigit::kDdlOffset+i));  //first is 0x700
 249     file[i].write((char*)&header,sizeof(header));                //write dummy header as place holder, actual will be written later when total size of DDL is known
 250     cnt[i]=0; //reset counters
 251   }
 252
 253   for(Int_t iChN=1;iChN<=kNchambers;iChN++){ //digits are stored on chamber by chamber basis
 254     for(Int_t iDigN=0;iDigN<Digits(iChN)->GetEntriesFast();iDigN++){//digits loop for a given chamber
 255       AliRICHDigit *pDig=(AliRICHDigit*)Digits(iChN)->At(iDigN);
 256       Int_t ddl=pDig->Dig2Raw(w32);  //ddl is 0..13
 257       file[ddl].write((char*)&w32,sizeof(w32));  cnt[ddl]++;//write formated digit to the propriate file (as decided in Dig2Raw) and increment corresponding counter
 258     }//digits loop for a given chamber
 259   }//chambers loop
 260   for(Int_t i=0;i<AliRICHDigit::kNddls;i++){
 261     header.fSize=sizeof(header)+cnt[i]*sizeof(w32); //now calculate total number of bytes for each DDL file
 262     header.SetAttribute(0);
 263     file[i].seekp(0); file[i].write((char*)&header,sizeof(header));//rewrite DDL header with fSize field properly set
 264     file[i].close();                                               //close DDL file
 265   }
 266   GetLoader()->UnloadDigits();
 267   AliDebug(1,"Stop.");
 268 }//Digits2Raw()
 269 //__________________________________________________________________________________________________