]>
Commit | Line | Data |
---|---|---|
53fd478b | 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 | // ************************************************************************** | |
c28632f0 | 15 | |
c28632f0 | 16 | |
0422a446 | 17 | #include "AliRICHv1.h" //class header |
53fd478b | 18 | #include "AliRICHParam.h" |
19 | #include "AliRICHChamber.h" | |
237c933d | 20 | #include <TParticle.h> |
88cb7938 | 21 | #include <TRandom.h> |
88cb7938 | 22 | #include <TVirtualMC.h> |
d128c9d1 | 23 | #include <TPDGCode.h> |
0422a446 | 24 | #include <AliStack.h> //Hits2SDigits() |
c60862bf | 25 | #include <AliConst.h> |
26 | #include <AliPDG.h> | |
27 | #include <AliRun.h> | |
5d12ce38 | 28 | #include <AliMC.h> |
0422a446 | 29 | #include <AliRawDataHeader.h> //Digits2Raw() |
c28632f0 | 30 | |
d128c9d1 | 31 | ClassImp(AliRICHv1) |
3582c1f9 | 32 | //__________________________________________________________________________________________________ |
d128c9d1 | 33 | void AliRICHv1::StepManager() |
53fd478b | 34 | { |
e42a7b46 | 35 | // Full Step Manager. |
998b831f | 36 | // 3- Ckovs absorbed on Collection electrods |
37 | // 5- Ckovs absorbed on Cathode wires | |
38 | // 6- Ckovs absorbed on Anod wires | |
e42a7b46 | 39 | |
ed3ceb24 | 40 | Int_t copy; |
c60862bf | 41 | static Int_t iCurrentChamber; |
d6fb41ac | 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"); | |
e42a7b46 | 47 | //history of Cerenkovs |
48 | if(gMC->TrackPid()==kCerenkov){ | |
d6fb41ac | 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 | |
e42a7b46 | 53 | } |
54 | ||
3582c1f9 | 55 | //Treat photons |
56 | static TLorentzVector cerX4; | |
d6fb41ac | 57 | if((gMC->TrackPid()==kCerenkov||gMC->TrackPid()==kFeedback)&&gMC->CurrentVolID(copy)==idRPC){//photon in PC |
998b831f | 58 | if(gMC->Edep()>0){//photon in PC +DE |
e42a7b46 | 59 | if(IsLostByFresnel()){ |
998b831f | 60 | if(gMC->TrackPid()==kCerenkov) fCounters(7)++;// 7- Ckovs reflected from CsI |
e42a7b46 | 61 | gMC->StopTrack(); |
62 | return; | |
63 | } | |
998b831f | 64 | gMC->TrackPosition(cerX4); gMC->CurrentVolOffID(2,iCurrentChamber);//RICH-RPPF-RPC |
c60862bf | 65 | |
0fe8fa07 | 66 | HitAdd(iCurrentChamber,gAlice->GetMCApp()->GetCurrentTrackNumber(),cerX4.Vect(),cerX4.Vect());//HIT for PHOTON in conditions CF+CSI+DE |
998b831f | 67 | fCounters(8)++;//4- Ckovs converted to electron on CsI |
3582c1f9 | 68 | GenerateFeedbacks(iCurrentChamber); |
998b831f | 69 | }//photon in PC and DE >0 |
70 | }//photon in PC | |
ed3ceb24 | 71 | |
72 | //Treat charged particles | |
73 | static Float_t eloss; | |
3582c1f9 | 74 | static TLorentzVector mipInX4,mipOutX4; |
d6fb41ac | 75 | if(gMC->TrackCharge() && gMC->CurrentVolID(copy)==idRGAP){//MIP in GAP |
998b831f | 76 | gMC->CurrentVolOffID(1,iCurrentChamber);//RICH-RGAP |
3582c1f9 | 77 | if(gMC->IsTrackEntering()||gMC->IsNewTrack()) {//MIP in GAP entering or newly created |
ed3ceb24 | 78 | eloss=0; |
af3d25a6 | 79 | gMC->TrackPosition(mipInX4); |
3582c1f9 | 80 | }else if(gMC->IsTrackExiting()||gMC->IsTrackStop()||gMC->IsTrackDisappeared()){//MIP in GAP exiting or disappeared |
ed3ceb24 | 81 | eloss+=gMC->Edep();//take into account last step dEdX |
af3d25a6 | 82 | gMC->TrackPosition(mipOutX4); |
0fe8fa07 | 83 | HitAdd(iCurrentChamber,gAlice->GetMCApp()->GetCurrentTrackNumber(),mipInX4.Vect(),mipOutX4.Vect(),eloss);//HIT for MIP: MIP in GAP Exiting |
ed3ceb24 | 84 | GenerateFeedbacks(iCurrentChamber,eloss);//MIP+GAP+Exit |
85 | }else//MIP in GAP going inside | |
86 | eloss += gMC->Edep(); | |
87 | }//MIP in GAP | |
3582c1f9 | 88 | }//StepManager() |
89 | //__________________________________________________________________________________________________ | |
90 | Bool_t AliRICHv1::IsLostByFresnel() | |
ed3ceb24 | 91 | { |
e42a7b46 | 92 | // Calculate probability for the photon to be lost by Fresnel reflection. |
3582c1f9 | 93 | TLorentzVector p4; |
94 | Double_t mom[3],localMom[3]; | |
ab296411 | 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; | |
3582c1f9 | 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)){ | |
998b831f | 102 | AliDebug(1,"Photon lost"); |
3582c1f9 | 103 | return kTRUE; |
104 | }else | |
105 | return kFALSE; | |
106 | }//IsLostByFresnel() | |
107 | //__________________________________________________________________________________________________ | |
e42a7b46 | 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); | |
998b831f | 115 | TVector2 x2=C(iChamber)->Mrs2Pc(x4);//hit position on photocathode plane |
116 | TVector2 xspe=x2; | |
0fe8fa07 | 117 | Int_t sector=P()->Loc2Sec(xspe); if(sector==-1) return; //hit in dead zone, nothing to produce |
e42a7b46 | 118 | Int_t iTotQdc=P()->TotQdc(x2,eloss); |
998b831f | 119 | Int_t iNphotons=gMC->GetRandom()->Poisson(P()->C(iChamber)->AlphaFeedback(sector)*iTotQdc); |
120 | AliDebug(1,Form("N photons=%i",iNphotons)); | |
e42a7b46 | 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 | |
998b831f | 182 | AliDebug(1,"Stop."); |
e42a7b46 | 183 | }//GenerateFeedbacks() |
0422a446 | 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 | //__________________________________________________________________________________________________ |