]> 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
Bug fix (I.Hrivnacova)
[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"
237c933d 30#include "AliRICHHit.h"
64e9b5aa 31#include "AliRICHResponse.h"
64e9b5aa 32#include "AliRICHResponseV0.h"
88cb7938 33#include "AliRICHSegmentationV1.h"
34#include "AliRICHv1.h"
c28632f0 35#include "AliRun.h"
c28632f0 36
d128c9d1 37ClassImp(AliRICHv1)
38//______________________________________________________________________________
c28632f0 39AliRICHv1::AliRICHv1(const char *name, const char *title)
d128c9d1 40 :AliRICH(name,title)
c28632f0 41{
c28632f0 42
237c933d 43// Full version of RICH with hits and diagnostics
44
64e9b5aa 45 // Version 0
46// Default Segmentation, no hits
03f3e21a 47 AliRICHSegmentationV1* segmentation = new AliRICHSegmentationV1;
64e9b5aa 48//
49// Segmentation parameters
c4384528 50 segmentation->SetPadSize(0.84,0.80);
51 segmentation->SetDAnod(0.84/2);
64e9b5aa 52//
53// Geometry parameters
54 AliRICHGeometry* geometry = new AliRICHGeometry;
55 geometry->SetGapThickness(8);
56 geometry->SetProximityGapThickness(.4);
2463d46d 57 geometry->SetQuartzLength(133);
58 geometry->SetQuartzWidth(127.9);
64e9b5aa 59 geometry->SetQuartzThickness(.5);
2463d46d 60 geometry->SetOuterFreonLength(133);
61 geometry->SetOuterFreonWidth(41.3);
62 geometry->SetInnerFreonLength(133);
63 geometry->SetInnerFreonWidth(41.3);
64 geometry->SetFreonThickness(1.5);
64e9b5aa 65//
66// Response parameters
c4384528 67 AliRICHResponseV0* response = new AliRICHResponseV0;
68 response->SetSigmaIntegration(5.);
69 response->SetChargeSlope(27.);
70 response->SetChargeSpread(0.18, 0.18);
71 response->SetMaxAdc(4096);
72 response->SetAlphaFeedback(0.036);
73 response->SetEIonisation(26.e-9);
74 response->SetSqrtKx3(0.77459667);
75 response->SetKx2(0.962);
76 response->SetKx4(0.379);
77 response->SetSqrtKy3(0.77459667);
78 response->SetKy2(0.962);
79 response->SetKy4(0.379);
80 response->SetPitch(0.25);
03f3e21a 81 response->SetWireSag(1); // 1->On, 0->Off
c4384528 82 response->SetVoltage(2150); // Should only be 2000, 2050, 2100 or 2150
83
64e9b5aa 84//
85// AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH");
86
237c933d 87 fCkovNumber=0;
88 fFreonProd=0;
64e9b5aa 89 Int_t i=0;
c28632f0 90
64e9b5aa 91 fChambers = new TObjArray(kNCH);
92 for (i=0; i<kNCH; i++) {
93
2682e810 94 //PH (*fChambers)[i] = new AliRICHChamber();
95 fChambers->AddAt(new AliRICHChamber(), i);
64e9b5aa 96
c28632f0 97 }
64e9b5aa 98
99 for (i=0; i<kNCH; i++) {
100 SetGeometryModel(i,geometry);
c4384528 101 SetSegmentationModel(i, segmentation);
102 SetResponseModel(i, response);
64e9b5aa 103 }
104
105
c28632f0 106}
107
237c933d 108void AliRICHv1::Init()
109{
110
9e1a0ddb 111 if(fDebug) {
112 printf("%s: *********************************** RICH_INIT ***********************************\n",ClassName());
113 printf("%s: * *\n",ClassName());
114 printf("%s: * AliRICHv1 Full version started *\n",ClassName());
115 printf("%s: * *\n",ClassName());
116 }
237c933d 117
118
a2f7eaf6 119 AliSegmentation* segmentation;
237c933d 120 AliRICHGeometry* geometry;
121 AliRICHResponse* response;
122
123
124 //
125 // Initialize Tracking Chambers
126 //
2fc6393f 127 for (Int_t i=0; i<kNCH; i++) {
237c933d 128 //printf ("i:%d",i);
2682e810 129 //PH ( (AliRICHChamber*) (*fChambers)[i])->Init(i);
130 ( (AliRICHChamber*)fChambers->At(i))->Init(i);
237c933d 131 }
132
133 //
134 // Set the chamber (sensitive region) GEANT identifier
135
2682e810 136 //PH ((AliRICHChamber*)(*fChambers)[0])->SetGid(1);
137 //PH ((AliRICHChamber*)(*fChambers)[1])->SetGid(2);
138 //PH ((AliRICHChamber*)(*fChambers)[2])->SetGid(3);
139 //PH ((AliRICHChamber*)(*fChambers)[3])->SetGid(4);
140 //PH ((AliRICHChamber*)(*fChambers)[4])->SetGid(5);
141 //PH ((AliRICHChamber*)(*fChambers)[5])->SetGid(6);
142 //PH ((AliRICHChamber*)(*fChambers)[6])->SetGid(7);
143
144 ((AliRICHChamber*)fChambers->At(0))->SetGid(1);
145 ((AliRICHChamber*)fChambers->At(1))->SetGid(2);
146 ((AliRICHChamber*)fChambers->At(2))->SetGid(3);
147 ((AliRICHChamber*)fChambers->At(3))->SetGid(4);
148 ((AliRICHChamber*)fChambers->At(4))->SetGid(5);
149 ((AliRICHChamber*)fChambers->At(5))->SetGid(6);
150 ((AliRICHChamber*)fChambers->At(6))->SetGid(7);
151
237c933d 152
237c933d 153 segmentation=Chamber(0).GetSegmentationModel(0);
154 geometry=Chamber(0).GetGeometryModel();
155 response=Chamber(0).GetResponseModel();
156
449e6743 157 Float_t offset = 490 + 1.276 - geometry->GetGapThickness()/2; //distance from center of mother volume to methane
158 Float_t deltaphi = 19.5; //phi angle between center of chambers - z direction
159 Float_t deltatheta = 20; //theta angle between center of chambers - x direction
160 Float_t cosphi = TMath::Cos(deltaphi*TMath::Pi()/180);
161 Float_t sinphi = TMath::Sin(deltaphi*TMath::Pi()/180);
162 Float_t costheta = TMath::Cos(deltatheta*TMath::Pi()/180);
163 Float_t sintheta = TMath::Sin(deltatheta*TMath::Pi()/180);
164
165 Float_t pos1[3]={0. , offset*cosphi , offset*sinphi};
166 Float_t pos2[3]={offset*sintheta , offset*costheta , 0. };
167 Float_t pos3[3]={0. , offset , 0.};
168 Float_t pos4[3]={-offset*sintheta , offset*costheta , 0.};
169 Float_t pos5[3]={offset*sinphi , offset*costheta*cosphi, -offset*sinphi};
170 Float_t pos6[3]={0. , offset*cosphi , -offset*sinphi};
171 Float_t pos7[3]={ -offset*sinphi , offset*costheta*cosphi, -offset*sinphi};
172
173 Chamber(0).SetChamberTransform(pos1[0],pos1[1],pos1[2],new TRotMatrix("rot993","rot993",90., 0. , 90. - deltaphi, 90. , deltaphi, -90. ));
174 Chamber(1).SetChamberTransform(pos2[0],pos2[1],pos2[2],new TRotMatrix("rot994","rot994",90., -deltatheta , 90. , 90.- deltatheta , 0. , 0. ));
175 Chamber(2).SetChamberTransform(pos3[0],pos3[1],pos3[2],new TRotMatrix("rot995","rot995",90., 0. , 90. , 90. , 0. , 0. ));
176 Chamber(3).SetChamberTransform(pos4[0],pos4[1],pos4[2],new TRotMatrix("rot996","rot996",90., deltatheta , 90. , 90 + deltatheta , 0. , 0. ));
177 Chamber(4).SetChamberTransform(pos5[0],pos5[1],pos5[2],new TRotMatrix("rot997","rot997",90., 360. - deltatheta, 108.2 , 90.- deltatheta ,18.2 , 90 - deltatheta));
178 Chamber(5).SetChamberTransform(pos6[0],pos6[1],pos6[2],new TRotMatrix("rot998","rot998",90., 0. , 90 + deltaphi , 90. , deltaphi, 90. ));
179 Chamber(6).SetChamberTransform(pos7[0],pos7[1],pos7[2],new TRotMatrix("rot999","rot999",90., deltatheta , 108.2 , 90.+ deltatheta ,18.2 , 90 + deltatheta));
237c933d 180
d128c9d1 181 if(GetDebug()) {
9e1a0ddb 182 printf("%s: * Pads : %3dx%3d *\n",
183 ClassName(),segmentation->Npx(),segmentation->Npy());
184 printf("%s: * Pad size : %5.2f x%5.2f mm2 *\n",
185 ClassName(),segmentation->Dpx(),segmentation->Dpy());
186 printf("%s: * Gap Thickness : %5.1f cm *\n",
187 ClassName(),geometry->GetGapThickness());
188 printf("%s: * Radiator Width : %5.1f cm *\n",
189 ClassName(),geometry->GetQuartzWidth());
190 printf("%s: * Radiator Length : %5.1f cm *\n",
191 ClassName(),geometry->GetQuartzLength());
192 printf("%s: * Freon Thickness : %5.1f cm *\n",
193 ClassName(),geometry->GetFreonThickness());
194 printf("%s: * Charge Slope : %5.1f ADC *\n",
195 ClassName(),response->ChargeSlope());
196 printf("%s: * Feedback Prob. : %5.2f %% *\n",
197 ClassName(),response->AlphaFeedback()*100);
9e1a0ddb 198 printf("%s: * *\n",
199 ClassName());
200 printf("%s: *********************************************************************************\n",
201 ClassName());
202 }
237c933d 203}
c28632f0 204
d128c9d1 205//______________________________________________________________________________
206void AliRICHv1::StepManager()
207{//Full Step Manager
208
209 Int_t copy, id;
210 static Int_t idvol;
211 static Int_t vol[2];
212 Int_t ipart;
213 static Float_t hits[22];
214 static Float_t ckovData[19];
215 TLorentzVector position;
216 TLorentzVector momentum;
217 Float_t pos[3];
218 Float_t mom[4];
219 Float_t localPos[3];
220 Float_t localMom[4];
221 Float_t localTheta,localPhi;
222 Float_t theta,phi;
223 Float_t destep, step;
224 Double_t ranf[2];
225 Int_t nPads;
226 Float_t coscerenkov;
227 static Float_t eloss, xhit, yhit, tlength;
228 const Float_t kBig=1.e10;
229
230 TClonesArray &lhits = *fHits;
231 TParticle *current = (TParticle*)(*gAlice->Particles())[gAlice->CurrentTrack()];
232
233 //if (current->Energy()>1)
234 //{
235
236 // Only gas gap inside chamber
237 // Tag chambers and record hits when track enters
238
239
240 id=gMC->CurrentVolID(copy);
241 idvol = copy-1;
242 Float_t cherenkovLoss=0;
243 //gAlice->KeepTrack(gAlice->CurrentTrack());
244
245 gMC->TrackPosition(position);
246 pos[0]=position(0);
247 pos[1]=position(1);
248 pos[2]=position(2);
249 //bzero((char *)ckovData,sizeof(ckovData)*19);
250 ckovData[1] = pos[0]; // X-position for hit
251 ckovData[2] = pos[1]; // Y-position for hit
252 ckovData[3] = pos[2]; // Z-position for hit
253 ckovData[6] = 0; // dummy track length
254 //ckovData[11] = gAlice->CurrentTrack();
255
256 //printf("\n+++++++++++\nTrack: %d\n++++++++++++\n",gAlice->CurrentTrack());
257
258 //AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH");
259
260 /********************Store production parameters for Cerenkov photons************************/
261//is it a Cerenkov photon?
262 if (gMC->TrackPid() == 50000050) {
263
264 //if (gMC->VolId("GAP ")==gMC->CurrentVolID(copy))
265 //{
266 Float_t ckovEnergy = current->Energy();
267 //energy interval for tracking
268 if (ckovEnergy > 5.6e-09 && ckovEnergy < 7.8e-09 )
269 //if (ckovEnergy > 0)
270 {
271 if (gMC->IsTrackEntering()){ //is track entering?
272 //printf("Track entered (1)\n");
273 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
274 { //is it in freo?
275 if (gMC->IsNewTrack()){ //is it the first step?
276 //printf("I'm in!\n");
277 Int_t mother = current->GetFirstMother();
278
279 //printf("Second Mother:%d\n",current->GetSecondMother());
280
281 ckovData[10] = mother;
282 ckovData[11] = gAlice->CurrentTrack();
283 ckovData[12] = 1; //Media where photon was produced 1->Freon, 2->Quarz
284 //printf("Produced in FREO\n");
285 fCkovNumber++;
286 fFreonProd=1;
287 //printf("Index: %d\n",fCkovNumber);
288 } //first step question
289 } //freo question
290
291 if (gMC->IsNewTrack()){ //is it first step?
292 if (gMC->VolId("QUAR")==gMC->CurrentVolID(copy)) //is it in quarz?
293 {
294 ckovData[12] = 2;
295 //printf("Produced in QUAR\n");
296 } //quarz question
297 } //first step question
298
299 //printf("Before %d\n",fFreonProd);
300 } //track entering question
301
302 if (ckovData[12] == 1) //was it produced in Freon?
303 //if (fFreonProd == 1)
304 {
305 if (gMC->IsTrackEntering()){ //is track entering?
306 //printf("Track entered (2)\n");
307 //printf("Current volume (should be META): %s\n",gMC->CurrentVolName());
308 //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("META"),gMC->CurrentVolID(copy));
309 if (gMC->VolId("META")==gMC->CurrentVolID(copy)) //is it in gap?
310 {
311 //printf("Got in META\n");
312 gMC->TrackMomentum(momentum);
313 mom[0]=momentum(0);
314 mom[1]=momentum(1);
315 mom[2]=momentum(2);
316 mom[3]=momentum(3);
317
318 gMC->Gmtod(mom,localMom,2);
319 Float_t cophi = TMath::Cos(TMath::ATan2(localMom[0], localMom[1]));
320 Float_t t = (1. - .025 / cophi) * (1. - .05 / cophi);
321 /**************** Photons lost in second grid have to be calculated by hand************/
322 gMC->GetRandom()->RndmArray(1,ranf);
323 if (ranf[0] > t) {
324 gMC->StopTrack();
325 ckovData[13] = 5;
326 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
327 //printf("Added One (1)!\n");
328 //printf("Lost one in grid\n");
329 }
330 /**********************************************************************************/
331 } //gap
332
333 //printf("Current volume (should be CSI) (1): %s\n",gMC->CurrentVolName());
334 //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
335 if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy)) //is it in csi?
336 {
337 //printf("Got in CSI\n");
338 gMC->TrackMomentum(momentum);
339 mom[0]=momentum(0);
340 mom[1]=momentum(1);
341 mom[2]=momentum(2);
342 mom[3]=momentum(3);
343
344 gMC->Gmtod(mom,localMom,2);
345 /********* Photons lost by Fresnel reflection have to be calculated by hand********/
346 /***********************Cerenkov phtons (always polarised)*************************/
347 Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
348 Double_t localRt = TMath::Sqrt(localTc);
349 localTheta = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])));
350 Double_t cotheta = TMath::Abs(cos(localTheta));
351 Float_t t = Fresnel(ckovEnergy*1e9,cotheta,1);
352 gMC->GetRandom()->RndmArray(1,ranf);
353 if (ranf[0] < t) {
354 gMC->StopTrack();
355 ckovData[13] = 6;
356 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
357
358 //printf("Added One (2)!\n");
359 //printf("Lost by Fresnel\n");
360 }
361 /**********************************************************************************/
362 }
363 } //track entering?
364
365
366 /********************Evaluation of losses************************/
367 /******************still in the old fashion**********************/
368
369 TArrayI procs;
370 Int_t i1 = gMC->StepProcesses(procs); //number of physics mechanisms acting on the particle
371 for (Int_t i = 0; i < i1; ++i) {
372 // Reflection loss
373 if (procs[i] == kPLightReflection) { //was it reflected
374 ckovData[13]=10;
375 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
376 ckovData[13]=1;
377 if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
378 ckovData[13]=2;
379 //gMC->StopTrack();
380 //AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
381 } //reflection question
382
383 // Absorption loss
384 else if (procs[i] == kPLightAbsorption) { //was it absorbed?
385 //printf("Got in absorption\n");
386 ckovData[13]=20;
387 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
388 ckovData[13]=11;
389 if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
390 ckovData[13]=12;
391 if (gMC->CurrentVolID(copy) == gMC->VolId("META"))
392 ckovData[13]=13;
393 if (gMC->CurrentVolID(copy) == gMC->VolId("GAP "))
394 ckovData[13]=13;
395
396 if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC"))
397 ckovData[13]=15;
398
399 // CsI inefficiency
400 if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
401 ckovData[13]=16;
402 }
403 gMC->StopTrack();
404 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
405 //printf("Added One (3)!\n");
406 //printf("Added cerenkov %d\n",fCkovNumber);
407 } //absorption question
408
409
410 // Photon goes out of tracking scope
411 else if (procs[i] == kPStop) { //is it below energy treshold?
412 ckovData[13]=21;
413 gMC->StopTrack();
414 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
415 //printf("Added One (4)!\n");
416 } // energy treshold question
417 } //number of mechanisms cycle
418 /**********************End of evaluation************************/
419 } //freon production question
420 } //energy interval question
421 //}//inside the proximity gap question
422 } //cerenkov photon question
423
424 /**************************************End of Production Parameters Storing*********************/
425
426
427 /*******************************Treat photons that hit the CsI (Ckovs and Feedbacks)************/
428
429 if (gMC->TrackPid() == 50000050 || gMC->TrackPid() == 50000051) {
430 //printf("Cerenkov\n");
431
432 //if (gMC->TrackPid() == 50000051)
433 //printf("Tracking a feedback\n");
434
435 if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
436 {
437 //printf("Current volume (should be CSI) (2): %s\n",gMC->CurrentVolName());
438 //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
439 //printf("Got in CSI\n");
440 //printf("Tracking a %d\n",gMC->TrackPid());
441 if (gMC->Edep() > 0.){
442 gMC->TrackPosition(position);
443 gMC->TrackMomentum(momentum);
444 pos[0]=position(0);
445 pos[1]=position(1);
446 pos[2]=position(2);
447 mom[0]=momentum(0);
448 mom[1]=momentum(1);
449 mom[2]=momentum(2);
450 mom[3]=momentum(3);
451 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
452 Double_t rt = TMath::Sqrt(tc);
453 theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
454 phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
455
456 gMC->CurrentVolOffID(2,copy);
457 vol[0]=copy;
458 idvol=vol[0]-1;
459
460
461 gMC->Gmtod(pos,localPos,1);
462
463 //Chamber(idvol).GlobaltoLocal(pos,localPos);
464
465 gMC->Gmtod(mom,localMom,2);
466
467 //Chamber(idvol).GlobaltoLocal(mom,localMom);
468
469 gMC->CurrentVolOffID(2,copy);
470 vol[0]=copy;
471 idvol=vol[0]-1;
472
473 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
474 //->Sector(localPos[0], localPos[2]);
475 //printf("Sector:%d\n",sector);
476
477 /*if (gMC->TrackPid() == 50000051){
478 fFeedbacks++;
479 printf("Feedbacks:%d\n",fFeedbacks);
480 }*/
481
482 //PH ((AliRICHChamber*) (*fChambers)[idvol])
483 ((AliRICHChamber*)fChambers->At(idvol))
484 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
485 if(idvol<kNCH) {
486 ckovData[0] = gMC->TrackPid(); // particle type
487 ckovData[1] = pos[0]; // X-position for hit
488 ckovData[2] = pos[1]; // Y-position for hit
489 ckovData[3] = pos[2]; // Z-position for hit
490 ckovData[4] = theta; // theta angle of incidence
491 ckovData[5] = phi; // phi angle of incidence
492 ckovData[8] = (Float_t) fNSDigits; // first sdigit
493 ckovData[9] = -1; // last pad hit
494 ckovData[13] = 4; // photon was detected
495 ckovData[14] = mom[0];
496 ckovData[15] = mom[1];
497 ckovData[16] = mom[2];
498
499 destep = gMC->Edep();
500 gMC->SetMaxStep(kBig);
501 cherenkovLoss += destep;
502 ckovData[7]=cherenkovLoss;
503
504 //nPads = Hits2SDigits(localPos[0],localPos[2],cherenkovLoss,idvol,kCerenkov);//for photons in CsI kir
505
506 if (fNSDigits > (Int_t)ckovData[8]) {
507 ckovData[8]= ckovData[8]+1;
508 ckovData[9]= (Float_t) fNSDigits;
509 }
510
511 //printf("Cerenkov loss: %f\n", cherenkovLoss);
512
513 ckovData[17] = nPads;
514 //printf("nPads:%d",nPads);
515
516 //TClonesArray *Hits = RICH->Hits();
517 AliRICHHit *mipHit = (AliRICHHit*) (fHits->UncheckedAt(0));
518 if (mipHit)
519 {
520 mom[0] = current->Px();
521 mom[1] = current->Py();
522 mom[2] = current->Pz();
523 Float_t mipPx = mipHit->MomX();
524 Float_t mipPy = mipHit->MomY();
525 Float_t mipPz = mipHit->MomZ();
526
527 Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
528 Float_t rt = TMath::Sqrt(r);
529 Float_t mipR = mipPx*mipPx + mipPy*mipPy + mipPz*mipPz;
530 Float_t mipRt = TMath::Sqrt(mipR);
531 if ((rt*mipRt) > 0)
532 {
533 coscerenkov = (mom[0]*mipPx + mom[1]*mipPy + mom[2]*mipPz)/(rt*mipRt);
534 }
535 else
536 {
537 coscerenkov = 0;
538 }
539 Float_t cherenkov = TMath::ACos(coscerenkov);
540 ckovData[18]=cherenkov;
541 }
542 //if (sector != -1)
543 //{
544 AddHit(gAlice->CurrentTrack(),vol,ckovData);
545 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
546 //printf("Added One (5)!\n");
547 //}
548 }
549 }
550 }
551 }
552
553 /***********************************************End of photon hits*********************************************/
554
555
556 /**********************************************Charged particles treatment*************************************/
557
558 else if (gMC->TrackCharge())
559 //else if (1 == 1)
560 {
561//If MIP
562 /*if (gMC->IsTrackEntering())
563 {
564 hits[13]=20;//is track entering?
565 }*/
566 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
567 {
568 gMC->TrackMomentum(momentum);
569 mom[0]=momentum(0);
570 mom[1]=momentum(1);
571 mom[2]=momentum(2);
572 mom[3]=momentum(3);
573 hits [19] = mom[0];
574 hits [20] = mom[1];
575 hits [21] = mom[2];
576 fFreonProd=1;
577 }
578
579 if (gMC->VolId("GAP ")== gMC->CurrentVolID(copy)) {
580// Get current particle id (ipart), track position (pos) and momentum (mom)
581
582 gMC->CurrentVolOffID(3,copy);
583 vol[0]=copy;
584 idvol=vol[0]-1;
585
586 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
587 //->Sector(localPos[0], localPos[2]);
588 //printf("Sector:%d\n",sector);
589
590 gMC->TrackPosition(position);
591 gMC->TrackMomentum(momentum);
592 pos[0]=position(0);
593 pos[1]=position(1);
594 pos[2]=position(2);
595 mom[0]=momentum(0);
596 mom[1]=momentum(1);
597 mom[2]=momentum(2);
598 mom[3]=momentum(3);
599
600 gMC->Gmtod(pos,localPos,1);
601
602 //Chamber(idvol).GlobaltoLocal(pos,localPos);
603
604 gMC->Gmtod(mom,localMom,2);
605
606 //Chamber(idvol).GlobaltoLocal(mom,localMom);
607
608 ipart = gMC->TrackPid();
609 //
610 // momentum loss and steplength in last step
611 destep = gMC->Edep();
612 step = gMC->TrackStep();
613
614 //
615 // record hits when track enters ...
616 if( gMC->IsTrackEntering()) {
617// gMC->SetMaxStep(fMaxStepGas);
618 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
619 Double_t rt = TMath::Sqrt(tc);
620 theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
621 phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
622
623
624 Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
625 Double_t localRt = TMath::Sqrt(localTc);
626 localTheta = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])))*kRaddeg;
627 localPhi = Float_t(TMath::ATan2(Double_t(localMom[2]),Double_t(localMom[0])))*kRaddeg;
628
629 hits[0] = Float_t(ipart); // particle type
630 hits[1] = localPos[0]; // X-position for hit
631 hits[2] = localPos[1]; // Y-position for hit
632 hits[3] = localPos[2]; // Z-position for hit
633 hits[4] = localTheta; // theta angle of incidence
634 hits[5] = localPhi; // phi angle of incidence
635 hits[8] = (Float_t) fNSDigits; // first sdigit
636 hits[9] = -1; // last pad hit
637 hits[13] = fFreonProd; // did id hit the freon?
638 hits[14] = mom[0];
639 hits[15] = mom[1];
640 hits[16] = mom[2];
641 hits[18] = 0; // dummy cerenkov angle
642
643 tlength = 0;
644 eloss = 0;
645 fFreonProd = 0;
646
647 Chamber(idvol).LocaltoGlobal(localPos,hits+1);
648
649
650 //To make chamber coordinates x-y had to pass localPos[0], localPos[2]
651 xhit = localPos[0];
652 yhit = localPos[2];
653 // Only if not trigger chamber
654 if(idvol<kNCH) {
655 //
656 // Initialize hit position (cursor) in the segmentation model
657 //PH ((AliRICHChamber*) (*fChambers)[idvol])
658 ((AliRICHChamber*)fChambers->At(idvol))
659 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
660 }
661 }
662
663 //
664 // Calculate the charge induced on a pad (disintegration) in case
665 //
666 // Mip left chamber ...
667 if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
668 gMC->SetMaxStep(kBig);
669 eloss += destep;
670 tlength += step;
671
672
673 // Only if not trigger chamber
674 if(idvol<kNCH) {
675 if (eloss > 0)
676 {
677 if(gMC->TrackPid() == kNeutron)
678 printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
679 //nPads = Hits2SDigits(xhit,yhit,eloss,idvol,kMip); //for MIP kir
680 hits[17] = nPads;
681 //printf("nPads:%d",nPads);
682 }
683 }
684
685 hits[6]=tlength;
686 hits[7]=eloss;
687 if (fNSDigits > (Int_t)hits[8]) {
688 hits[8]= hits[8]+1;
689 hits[9]= (Float_t) fNSDigits;
690 }
691
692 //if(sector !=-1)
693 new(lhits[fNhits++]) AliRICHHit(fIshunt,gAlice->CurrentTrack(),vol,hits);
694 eloss = 0;
695 //
696 // Check additional signal generation conditions
697 // defined by the segmentation
698 // model (boundary crossing conditions)
699 } else if
700 //PH (((AliRICHChamber*) (*fChambers)[idvol])
701 (((AliRICHChamber*)fChambers->At(idvol))
702 ->SigGenCond(localPos[0], localPos[2], localPos[1]))
703 {
704 //PH ((AliRICHChamber*) (*fChambers)[idvol])
705 ((AliRICHChamber*)fChambers->At(idvol))
706 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
707 if (eloss > 0)
708 {
709 if(gMC->TrackPid() == kNeutron)
710 printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
711 //nPads = Hits2SDigits(xhit,yhit,eloss,idvol,kMip);//for N kir
712 hits[17] = nPads;
713 //printf("Npads:%d",NPads);
714 }
715 xhit = localPos[0];
716 yhit = localPos[2];
717 eloss = destep;
718 tlength += step ;
719 //
720 // nothing special happened, add up energy loss
721 } else {
722 eloss += destep;
723 tlength += step ;
724 }
725 }
726 }
727 /*************************************************End of MIP treatment**************************************/
728}//void AliRICHv1::StepManager()
2463d46d 729
730
731
732
733
734
735