]> git.uio.no Git - u/mrichter/AliRoot.git/blame - RICH/AliRICHv1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New Stepmanager function.
[u/mrichter/AliRoot.git] / RICH / AliRICHv1.cxx
CommitLineData
c28632f0 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
88cb7938 16/* $Id$ */
c28632f0 17
88cb7938 18#include "Riostream.h"
19
c28632f0 20#include <TNode.h>
237c933d 21#include <TParticle.h>
88cb7938 22#include <TRandom.h>
23#include <TTUBE.h>
24#include <TVirtualMC.h>
d128c9d1 25#include <TPDGCode.h>
c28632f0 26
88cb7938 27#include "AliConst.h"
28#include "AliPDG.h"
29#include "AliRICHGeometry.h"
64e9b5aa 30#include "AliRICHResponse.h"
64e9b5aa 31#include "AliRICHResponseV0.h"
88cb7938 32#include "AliRICHSegmentationV1.h"
33#include "AliRICHv1.h"
c28632f0 34#include "AliRun.h"
c28632f0 35
d128c9d1 36ClassImp(AliRICHv1)
37//______________________________________________________________________________
c28632f0 38AliRICHv1::AliRICHv1(const char *name, const char *title)
d128c9d1 39 :AliRICH(name,title)
3ea9cb08 40{// Full version of RICH with hits and diagnostics
41 if(GetDebug())Info("named ctor","Start.");
42 if(GetDebug())Info("named ctor","Stop.");
43}//named ctor
44//______________________________________________________________________________
237c933d 45void AliRICHv1::Init()
3ea9cb08 46{//nothing to do here, all the work in ctor or CreateGeometry()
47 if(GetDebug())Info("Init","Start.");
48 if(GetDebug())Info("Init","Stop.");
49}//void AliRICHv1::Init()
d128c9d1 50//______________________________________________________________________________
51void AliRICHv1::StepManager()
52{//Full Step Manager
53
54 Int_t copy, id;
55 static Int_t idvol;
56 static Int_t vol[2];
57 Int_t ipart;
58 static Float_t hits[22];
59 static Float_t ckovData[19];
60 TLorentzVector position;
61 TLorentzVector momentum;
62 Float_t pos[3];
63 Float_t mom[4];
64 Float_t localPos[3];
65 Float_t localMom[4];
66 Float_t localTheta,localPhi;
67 Float_t theta,phi;
68 Float_t destep, step;
69 Double_t ranf[2];
942194a4 70 Int_t nPads=-1;
d128c9d1 71 Float_t coscerenkov;
72 static Float_t eloss, xhit, yhit, tlength;
73 const Float_t kBig=1.e10;
74
75 TClonesArray &lhits = *fHits;
642f15cf 76 TParticle *current = (TParticle*)(*gAlice->Particles())[gAlice->GetCurrentTrackNumber()];
d128c9d1 77
d128c9d1 78
79
80 id=gMC->CurrentVolID(copy);
81 idvol = copy-1;
82 Float_t cherenkovLoss=0;
d128c9d1 83
84 gMC->TrackPosition(position);
85 pos[0]=position(0);
86 pos[1]=position(1);
87 pos[2]=position(2);
d128c9d1 88 ckovData[1] = pos[0]; // X-position for hit
89 ckovData[2] = pos[1]; // Y-position for hit
90 ckovData[3] = pos[2]; // Z-position for hit
91 ckovData[6] = 0; // dummy track length
d128c9d1 92
93 /********************Store production parameters for Cerenkov photons************************/
853634d3 94
95 if (gMC->TrackPid() == 50000050) { //is it a Cerenkov photon?
d128c9d1 96
97 //if (gMC->VolId("GAP ")==gMC->CurrentVolID(copy))
98 //{
99 Float_t ckovEnergy = current->Energy();
100 //energy interval for tracking
101 if (ckovEnergy > 5.6e-09 && ckovEnergy < 7.8e-09 )
102 //if (ckovEnergy > 0)
103 {
104 if (gMC->IsTrackEntering()){ //is track entering?
105 //printf("Track entered (1)\n");
106 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
107 { //is it in freo?
108 if (gMC->IsNewTrack()){ //is it the first step?
109 //printf("I'm in!\n");
110 Int_t mother = current->GetFirstMother();
111
d128c9d1 112
113 ckovData[10] = mother;
642f15cf 114 ckovData[11] = gAlice->GetCurrentTrackNumber();
d128c9d1 115 ckovData[12] = 1; //Media where photon was produced 1->Freon, 2->Quarz
116 //printf("Produced in FREO\n");
117 fCkovNumber++;
118 fFreonProd=1;
d128c9d1 119 } //first step question
120 } //freo question
121
122 if (gMC->IsNewTrack()){ //is it first step?
123 if (gMC->VolId("QUAR")==gMC->CurrentVolID(copy)) //is it in quarz?
124 {
125 ckovData[12] = 2;
d128c9d1 126 } //quarz question
127 } //first step question
128
d128c9d1 129 } //track entering question
130
131 if (ckovData[12] == 1) //was it produced in Freon?
132 //if (fFreonProd == 1)
133 {
134 if (gMC->IsTrackEntering()){ //is track entering?
d128c9d1 135 if (gMC->VolId("META")==gMC->CurrentVolID(copy)) //is it in gap?
136 {
137 //printf("Got in META\n");
138 gMC->TrackMomentum(momentum);
139 mom[0]=momentum(0);
140 mom[1]=momentum(1);
141 mom[2]=momentum(2);
142 mom[3]=momentum(3);
143
144 gMC->Gmtod(mom,localMom,2);
145 Float_t cophi = TMath::Cos(TMath::ATan2(localMom[0], localMom[1]));
146 Float_t t = (1. - .025 / cophi) * (1. - .05 / cophi);
147 /**************** Photons lost in second grid have to be calculated by hand************/
148 gMC->GetRandom()->RndmArray(1,ranf);
149 if (ranf[0] > t) {
150 gMC->StopTrack();
151 ckovData[13] = 5;
642f15cf 152 AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
d128c9d1 153 }
154 /**********************************************************************************/
155 } //gap
156
d128c9d1 157 if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy)) //is it in csi?
158 {
159 //printf("Got in CSI\n");
160 gMC->TrackMomentum(momentum);
161 mom[0]=momentum(0);
162 mom[1]=momentum(1);
163 mom[2]=momentum(2);
164 mom[3]=momentum(3);
165
166 gMC->Gmtod(mom,localMom,2);
167 /********* Photons lost by Fresnel reflection have to be calculated by hand********/
168 /***********************Cerenkov phtons (always polarised)*************************/
169 Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
170 Double_t localRt = TMath::Sqrt(localTc);
171 localTheta = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])));
172 Double_t cotheta = TMath::Abs(cos(localTheta));
173 Float_t t = Fresnel(ckovEnergy*1e9,cotheta,1);
174 gMC->GetRandom()->RndmArray(1,ranf);
175 if (ranf[0] < t) {
176 gMC->StopTrack();
177 ckovData[13] = 6;
642f15cf 178 AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
d128c9d1 179
180 //printf("Added One (2)!\n");
181 //printf("Lost by Fresnel\n");
182 }
183 /**********************************************************************************/
184 }
185 } //track entering?
186
187
188 /********************Evaluation of losses************************/
189 /******************still in the old fashion**********************/
190
191 TArrayI procs;
192 Int_t i1 = gMC->StepProcesses(procs); //number of physics mechanisms acting on the particle
193 for (Int_t i = 0; i < i1; ++i) {
194 // Reflection loss
195 if (procs[i] == kPLightReflection) { //was it reflected
196 ckovData[13]=10;
197 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
198 ckovData[13]=1;
199 if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
200 ckovData[13]=2;
d128c9d1 201 } //reflection question
202
203 // Absorption loss
204 else if (procs[i] == kPLightAbsorption) { //was it absorbed?
205 //printf("Got in absorption\n");
206 ckovData[13]=20;
207 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
208 ckovData[13]=11;
209 if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
210 ckovData[13]=12;
211 if (gMC->CurrentVolID(copy) == gMC->VolId("META"))
212 ckovData[13]=13;
213 if (gMC->CurrentVolID(copy) == gMC->VolId("GAP "))
214 ckovData[13]=13;
215
216 if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC"))
217 ckovData[13]=15;
218
219 // CsI inefficiency
220 if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
221 ckovData[13]=16;
222 }
223 gMC->StopTrack();
642f15cf 224 AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
d128c9d1 225 } //absorption question
226
227
228 // Photon goes out of tracking scope
229 else if (procs[i] == kPStop) { //is it below energy treshold?
230 ckovData[13]=21;
231 gMC->StopTrack();
642f15cf 232 AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
d128c9d1 233 } // energy treshold question
234 } //number of mechanisms cycle
235 /**********************End of evaluation************************/
236 } //freon production question
237 } //energy interval question
d128c9d1 238 } //cerenkov photon question
239
240 /**************************************End of Production Parameters Storing*********************/
241
242
243 /*******************************Treat photons that hit the CsI (Ckovs and Feedbacks)************/
244
245 if (gMC->TrackPid() == 50000050 || gMC->TrackPid() == 50000051) {
d128c9d1 246
247 if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
248 {
d128c9d1 249 if (gMC->Edep() > 0.){
250 gMC->TrackPosition(position);
251 gMC->TrackMomentum(momentum);
252 pos[0]=position(0);
253 pos[1]=position(1);
254 pos[2]=position(2);
255 mom[0]=momentum(0);
256 mom[1]=momentum(1);
257 mom[2]=momentum(2);
258 mom[3]=momentum(3);
259 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
260 Double_t rt = TMath::Sqrt(tc);
261 theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
262 phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
263
264 gMC->CurrentVolOffID(2,copy);
265 vol[0]=copy;
266 idvol=vol[0]-1;
267
268
269 gMC->Gmtod(pos,localPos,1);
270
d128c9d1 271
272 gMC->Gmtod(mom,localMom,2);
273
d128c9d1 274
275 gMC->CurrentVolOffID(2,copy);
276 vol[0]=copy;
277 idvol=vol[0]-1;
278
d128c9d1 279
d128c9d1 280 ((AliRICHChamber*)fChambers->At(idvol))
281 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
282 if(idvol<kNCH) {
283 ckovData[0] = gMC->TrackPid(); // particle type
284 ckovData[1] = pos[0]; // X-position for hit
285 ckovData[2] = pos[1]; // Y-position for hit
286 ckovData[3] = pos[2]; // Z-position for hit
287 ckovData[4] = theta; // theta angle of incidence
288 ckovData[5] = phi; // phi angle of incidence
853634d3 289 ckovData[8] = (Float_t) fNsdigits; // first sdigit
d128c9d1 290 ckovData[9] = -1; // last pad hit
291 ckovData[13] = 4; // photon was detected
292 ckovData[14] = mom[0];
293 ckovData[15] = mom[1];
294 ckovData[16] = mom[2];
295
296 destep = gMC->Edep();
297 gMC->SetMaxStep(kBig);
298 cherenkovLoss += destep;
299 ckovData[7]=cherenkovLoss;
300
301 //nPads = Hits2SDigits(localPos[0],localPos[2],cherenkovLoss,idvol,kCerenkov);//for photons in CsI kir
302
853634d3 303 if (fNsdigits > (Int_t)ckovData[8]) {
d128c9d1 304 ckovData[8]= ckovData[8]+1;
853634d3 305 ckovData[9]= (Float_t) fNsdigits;
d128c9d1 306 }
307
d128c9d1 308
309 ckovData[17] = nPads;
853634d3 310 AliRICHhit *mipHit = (AliRICHhit*) (fHits->UncheckedAt(0));
d128c9d1 311 if (mipHit)
312 {
313 mom[0] = current->Px();
314 mom[1] = current->Py();
315 mom[2] = current->Pz();
316 Float_t mipPx = mipHit->MomX();
317 Float_t mipPy = mipHit->MomY();
318 Float_t mipPz = mipHit->MomZ();
319
320 Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
321 Float_t rt = TMath::Sqrt(r);
322 Float_t mipR = mipPx*mipPx + mipPy*mipPy + mipPz*mipPz;
323 Float_t mipRt = TMath::Sqrt(mipR);
324 if ((rt*mipRt) > 0)
325 {
326 coscerenkov = (mom[0]*mipPx + mom[1]*mipPy + mom[2]*mipPz)/(rt*mipRt);
327 }
328 else
329 {
330 coscerenkov = 0;
331 }
332 Float_t cherenkov = TMath::ACos(coscerenkov);
333 ckovData[18]=cherenkov;
334 }
335 //if (sector != -1)
336 //{
642f15cf 337 AddHit(gAlice->GetCurrentTrackNumber(),vol,ckovData);
338 AddCerenkov(gAlice->GetCurrentTrackNumber(),vol,ckovData);
d128c9d1 339 //printf("Added One (5)!\n");
340 //}
341 }
342 }
343 }
344 }
345
346 /***********************************************End of photon hits*********************************************/
347
348
349 /**********************************************Charged particles treatment*************************************/
350
853634d3 351 else if (gMC->TrackCharge()){//is MIP?
d128c9d1 352 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
353 {
354 gMC->TrackMomentum(momentum);
355 mom[0]=momentum(0);
356 mom[1]=momentum(1);
357 mom[2]=momentum(2);
358 mom[3]=momentum(3);
359 hits [19] = mom[0];
360 hits [20] = mom[1];
361 hits [21] = mom[2];
362 fFreonProd=1;
363 }
364
853634d3 365 if(gMC->VolId("GAP ")==gMC->CurrentVolID(copy)) {//is in GAP?
d128c9d1 366// Get current particle id (ipart), track position (pos) and momentum (mom)
367
368 gMC->CurrentVolOffID(3,copy);
369 vol[0]=copy;
370 idvol=vol[0]-1;
d128c9d1 371 gMC->TrackPosition(position);
372 gMC->TrackMomentum(momentum);
373 pos[0]=position(0);
374 pos[1]=position(1);
375 pos[2]=position(2);
376 mom[0]=momentum(0);
377 mom[1]=momentum(1);
378 mom[2]=momentum(2);
379 mom[3]=momentum(3);
380
381 gMC->Gmtod(pos,localPos,1);
d128c9d1 382 gMC->Gmtod(mom,localMom,2);
d128c9d1 383 ipart = gMC->TrackPid();
853634d3 384 destep = gMC->Edep();step = gMC->TrackStep();// momentum loss and steplength in last step
385 if(gMC->IsTrackEntering()){ // record hits when track enters ...
386
d128c9d1 387 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
388 Double_t rt = TMath::Sqrt(tc);
389 theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
390 phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
391
392
393 Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
394 Double_t localRt = TMath::Sqrt(localTc);
395 localTheta = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])))*kRaddeg;
396 localPhi = Float_t(TMath::ATan2(Double_t(localMom[2]),Double_t(localMom[0])))*kRaddeg;
397
398 hits[0] = Float_t(ipart); // particle type
399 hits[1] = localPos[0]; // X-position for hit
400 hits[2] = localPos[1]; // Y-position for hit
401 hits[3] = localPos[2]; // Z-position for hit
402 hits[4] = localTheta; // theta angle of incidence
403 hits[5] = localPhi; // phi angle of incidence
853634d3 404 hits[8] = (Float_t) fNsdigits; // first sdigit
d128c9d1 405 hits[9] = -1; // last pad hit
406 hits[13] = fFreonProd; // did id hit the freon?
407 hits[14] = mom[0];
408 hits[15] = mom[1];
409 hits[16] = mom[2];
410 hits[18] = 0; // dummy cerenkov angle
411
412 tlength = 0;
413 eloss = 0;
414 fFreonProd = 0;
415
416 Chamber(idvol).LocaltoGlobal(localPos,hits+1);
417
418
419 //To make chamber coordinates x-y had to pass localPos[0], localPos[2]
420 xhit = localPos[0];
421 yhit = localPos[2];
422 // Only if not trigger chamber
423 if(idvol<kNCH) {
d128c9d1 424 // Initialize hit position (cursor) in the segmentation model
d128c9d1 425 ((AliRICHChamber*)fChambers->At(idvol))
426 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
427 }
428 }
429
d128c9d1 430 // Calculate the charge induced on a pad (disintegration) in case
d128c9d1 431 // Mip left chamber ...
432 if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
433 gMC->SetMaxStep(kBig);
434 eloss += destep;
435 tlength += step;
436
437
438 // Only if not trigger chamber
439 if(idvol<kNCH) {
440 if (eloss > 0)
441 {
442 if(gMC->TrackPid() == kNeutron)
443 printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
444 //nPads = Hits2SDigits(xhit,yhit,eloss,idvol,kMip); //for MIP kir
445 hits[17] = nPads;
d128c9d1 446 }
447 }
448
449 hits[6]=tlength;
450 hits[7]=eloss;
853634d3 451 if (fNsdigits > (Int_t)hits[8]) {
d128c9d1 452 hits[8]= hits[8]+1;
853634d3 453 hits[9]= (Float_t) fNsdigits;
d128c9d1 454 }
455
853634d3 456 new(lhits[fNhits++]) AliRICHhit(fIshunt,gAlice->GetCurrentTrackNumber(),vol,hits);
d128c9d1 457 eloss = 0;
d128c9d1 458 // Check additional signal generation conditions
459 // defined by the segmentation
460 // model (boundary crossing conditions)
461 } else if
d128c9d1 462 (((AliRICHChamber*)fChambers->At(idvol))
463 ->SigGenCond(localPos[0], localPos[2], localPos[1]))
464 {
d128c9d1 465 ((AliRICHChamber*)fChambers->At(idvol))
466 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
467 if (eloss > 0)
468 {
469 if(gMC->TrackPid() == kNeutron)
470 printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
471 //nPads = Hits2SDigits(xhit,yhit,eloss,idvol,kMip);//for N kir
472 hits[17] = nPads;
473 //printf("Npads:%d",NPads);
474 }
475 xhit = localPos[0];
476 yhit = localPos[2];
477 eloss = destep;
478 tlength += step ;
479 //
480 // nothing special happened, add up energy loss
481 } else {
482 eloss += destep;
483 tlength += step ;
484 }
853634d3 485 }//is in GAP?
486 }//is MIP?
d128c9d1 487}//void AliRICHv1::StepManager()