1 // **************************************************************************
2 // * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 // * Author: The ALICE Off-line Project. *
5 // * Contributors are mentioned in the code where appropriate. *
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 // **************************************************************************
17 #include "AliRICHv1.h" //class header
18 #include "AliRICHParam.h"
19 #include "AliRICHChamber.h"
20 #include <TParticle.h>
22 #include <TVirtualMC.h>
24 #include <AliStack.h> //Hits2SDigits()
29 #include <AliRawDataHeader.h> //Digits2Raw()
32 //__________________________________________________________________________________________________
33 void AliRICHv1::StepManager()
36 // 3- Ckovs absorbed on Collection electrods
37 // 5- Ckovs absorbed on Cathode wires
38 // 6- Ckovs absorbed on Anod wires
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
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
64 gMC->TrackPosition(cerX4); gMC->CurrentVolOffID(2,iCurrentChamber);//RICH-RPPF-RPC
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
72 //Treat charged particles
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
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
89 //__________________________________________________________________________________________________
90 Bool_t AliRICHv1::IsLostByFresnel()
92 // Calculate probability for the photon to be lost by Fresnel reflection.
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");
107 //__________________________________________________________________________________________________
108 void AliRICHv1::GenerateFeedbacks(Int_t iChamber,Float_t eloss)
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()
114 gMC->TrackPosition(x4);
115 TVector2 x2=C(iChamber)->Mrs2Pc(x4);//hit position on photocathode plane
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));
122 Float_t cthf, phif, enfp = 0, sthf, e1[3], e2[3], e3[3], vmod, uswop,dir[3], phi,pol[3], mom[4];
124 for(Int_t i=0;i<iNphotons;i++){//feedbacks loop
126 gMC->GetRandom()->RndmArray(2,ranf); //Sample direction
129 sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
130 phif = ranf[1] * 2 * TMath::Pi();
132 if(Double_t randomNumber=gMC->GetRandom()->Rndm()<=0.57)
134 else if(randomNumber<=0.7)
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]);
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];
151 for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
152 if (!vmod) for(j=0;j<3;j++) {
158 for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
159 if (!vmod) for(j=0;j<3;j++) {
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;
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
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
183 }//GenerateFeedbacks()
184 //__________________________________________________________________________________________________
185 void AliRICHv1::Hits2SDigits()
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
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
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;
204 //need to quantize the anod....
205 TVector padHit=AliRICHParam::Loc2Pad(x2);
206 TVector2 padHitXY=AliRICHParam::Pad2Loc(padHit);
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
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));
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
223 GetLoader()->TreeS()->Fill();
224 GetLoader()->WriteSDigits("OVERWRITE");
227 GetLoader()->UnloadHits(); GetLoader()->GetRunLoader()->UnloadHeader(); GetLoader()->GetRunLoader()->UnloadKinematics();
228 GetLoader()->UnloadSDigits();
231 //__________________________________________________________________________________________________
232 void AliRICHv1::Digits2Raw()
234 //Creates raw data files in DDL format. Invoked by AliSimulation
235 //loop over events is done outside in AliSimulation
238 AliDebug(1,"Start.");
239 GetLoader()->LoadDigits();
240 GetLoader()->TreeD()->GetEntry(0);
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
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
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
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
266 GetLoader()->UnloadDigits();
269 //__________________________________________________________________________________________________