Default arguments declared once
[u/mrichter/AliRoot.git] / RICH / AliRICH.cxx
CommitLineData
4c039060 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/*
2e5f0f7b 17 $Log$
cc683707 18 Revision 1.50 2001/05/10 14:44:16 jbarbosa
19 Corrected some overlaps (thanks I. Hrivnacovna).
20
4591f067 21 Revision 1.49 2001/05/10 12:23:49 jbarbosa
22 Repositioned the RICH modules.
23 Eliminated magic numbers.
24 Incorporated diagnostics (from macros).
25
a3d71079 26 Revision 1.48 2001/03/15 10:35:00 jbarbosa
27 Corrected bug in MakeBranch (was using a different version of STEER)
28
c3ecfac9 29 Revision 1.47 2001/03/14 18:13:56 jbarbosa
30 Several changes to adapt to new IO.
31 Removed digitising function, using AliRICHMerger::Digitise from now on.
32
34ead2dd 33 Revision 1.46 2001/03/12 17:46:33 hristov
34 Changes needed on Sun with CC 5.0
35
5cf7bbad 36 Revision 1.45 2001/02/27 22:11:46 jbarbosa
37 Testing TreeS, removing of output.
38
633e9a38 39 Revision 1.44 2001/02/27 15:19:12 jbarbosa
40 Transition to SDigits.
41
b251a2b5 42 Revision 1.43 2001/02/23 17:19:06 jbarbosa
43 Corrected photocathode definition in BuildGeometry().
44
ca96c9ea 45 Revision 1.42 2001/02/13 20:07:23 jbarbosa
46 Parametrised definition of photcathode dimensions. New spacers. New data members in AliRICHHit to store particle momentum
47 when entering the freon. Corrected calls to particle stack.
48
c24372d0 49 Revision 1.41 2001/01/26 20:00:20 hristov
50 Major upgrade of AliRoot code
51
2ab0c725 52 Revision 1.40 2001/01/24 20:58:03 jbarbosa
53 Enhanced BuildGeometry. Now the photocathodes are drawn.
54
91975aa2 55 Revision 1.39 2001/01/22 21:40:24 jbarbosa
56 Removing magic numbers
57
ee2105a4 58 Revision 1.37 2000/12/20 14:07:25 jbarbosa
59 Removed dependencies on TGeant3 (thanks to F. Carminati and I. Hrivnacova)
60
5cfcdf54 61 Revision 1.36 2000/12/18 17:45:54 jbarbosa
62 Cleaned up PadHits object.
63
176a9917 64 Revision 1.35 2000/12/15 16:49:40 jbarbosa
65 Geometry and materials updates (wire supports, pcbs, backplane supports, frame).
66
6197ac87 67 Revision 1.34 2000/11/10 18:12:12 jbarbosa
68 Bug fix for AliRICHCerenkov (thanks to P. Hristov)
69
3cc8edee 70 Revision 1.33 2000/11/02 10:09:01 jbarbosa
71 Minor bug correction (some pointers were not initialised in the default constructor)
72
9b674b76 73 Revision 1.32 2000/11/01 15:32:55 jbarbosa
74 Updated to handle both reconstruction algorithms.
75
a4622d0f 76 Revision 1.31 2000/10/26 20:18:33 jbarbosa
77 Supports for methane and freon vessels
78
683b273b 79 Revision 1.30 2000/10/24 13:19:12 jbarbosa
80 Geometry updates.
81
e293afae 82 Revision 1.29 2000/10/19 19:39:25 jbarbosa
83 Some more changes to geometry. Further correction of digitisation "per part. type"
84
53a60579 85 Revision 1.28 2000/10/17 20:50:57 jbarbosa
86 Inversed digtise by particle type (now, only the selected particle type is not digitsed).
87 Corrected several geometry minor bugs.
88 Added new parameter (opaque quartz thickness).
89
3587e146 90 Revision 1.27 2000/10/11 10:33:55 jbarbosa
91 Corrected bug introduced by earlier revisions (CerenkovData array cannot be reset to zero on wach call of StepManager)
92
8cda28a5 93 Revision 1.26 2000/10/03 21:44:08 morsch
94 Use AliSegmentation and AliHit abstract base classes.
95
a2f7eaf6 96 Revision 1.25 2000/10/02 21:28:12 fca
97 Removal of useless dependecies via forward declarations
98
94de3818 99 Revision 1.24 2000/10/02 15:43:17 jbarbosa
100 Fixed forward declarations.
101 Fixed honeycomb density.
102 Fixed cerenkov storing.
103 New electronics.
104
9dda3582 105 Revision 1.23 2000/09/13 10:42:14 hristov
106 Minor corrections for HP, DEC and Sun; strings.h included
107
e813258a 108 Revision 1.22 2000/09/12 18:11:13 fca
109 zero hits area before using
110
e51a3fa9 111 Revision 1.21 2000/07/21 10:21:07 morsch
112 fNrawch = 0; and fNrechits = 0; in the default constructor.
113
db481a6e 114 Revision 1.20 2000/07/10 15:28:39 fca
115 Correction of the inheritance scheme
116
452a64c6 117 Revision 1.19 2000/06/30 16:29:51 dibari
118 Added kDebugLevel variable to control output size on demand
119
33984590 120 Revision 1.18 2000/06/12 15:15:46 jbarbosa
121 Cleaned up version.
122
237c933d 123 Revision 1.17 2000/06/09 14:58:37 jbarbosa
124 New digitisation per particle type
125
d28b5fc2 126 Revision 1.16 2000/04/19 12:55:43 morsch
127 Newly structured and updated version (JB, AM)
128
4c039060 129*/
130
2e5f0f7b 131
ddae0931 132////////////////////////////////////////////////
133// Manager and hits classes for set:RICH //
134////////////////////////////////////////////////
fe4da5cc 135
ddae0931 136#include <TBRIK.h>
137#include <TTUBE.h>
fe4da5cc 138#include <TNode.h>
139#include <TRandom.h>
ddae0931 140#include <TObject.h>
141#include <TVector.h>
142#include <TObjArray.h>
2e5f0f7b 143#include <TArrayF.h>
237c933d 144#include <TFile.h>
145#include <TParticle.h>
94de3818 146#include <TGeometry.h>
147#include <TTree.h>
a3d71079 148#include <TH1.h>
149#include <TH2.h>
150#include <TCanvas.h>
151//#include <TPad.h>
152#include <TF1.h>
94de3818 153
237c933d 154#include <iostream.h>
e813258a 155#include <strings.h>
fe4da5cc 156
fe4da5cc 157#include "AliRICH.h"
a2f7eaf6 158#include "AliSegmentation.h"
91975aa2 159#include "AliRICHSegmentationV0.h"
237c933d 160#include "AliRICHHit.h"
161#include "AliRICHCerenkov.h"
b251a2b5 162#include "AliRICHSDigit.h"
237c933d 163#include "AliRICHDigit.h"
164#include "AliRICHTransientDigit.h"
165#include "AliRICHRawCluster.h"
a4622d0f 166#include "AliRICHRecHit1D.h"
167#include "AliRICHRecHit3D.h"
237c933d 168#include "AliRICHHitMapA1.h"
2e5f0f7b 169#include "AliRICHClusterFinder.h"
34ead2dd 170#include "AliRICHMerger.h"
fe4da5cc 171#include "AliRun.h"
ddae0931 172#include "AliMC.h"
94de3818 173#include "AliMagF.h"
452a64c6 174#include "AliConst.h"
175#include "AliPDG.h"
2e5f0f7b 176#include "AliPoints.h"
ddae0931 177#include "AliCallf77.h"
237c933d 178
ddae0931 179
180// Static variables for the pad-hit iterator routines
181static Int_t sMaxIterPad=0;
182static Int_t sCurIterPad=0;
ddae0931 183
fe4da5cc 184ClassImp(AliRICH)
ddae0931 185
186//___________________________________________
fe4da5cc 187AliRICH::AliRICH()
188{
237c933d 189// Default constructor for RICH manager class
190
ddae0931 191 fIshunt = 0;
192 fHits = 0;
b251a2b5 193 fSDigits = 0;
194 fNSDigits = 0;
2e5f0f7b 195 fNcerenkovs = 0;
ddae0931 196 fDchambers = 0;
9b674b76 197 fRecHits1D = 0;
198 fRecHits3D = 0;
199 fRawClusters = 0;
200 fChambers = 0;
ddae0931 201 fCerenkovs = 0;
9dda3582 202 for (Int_t i=0; i<7; i++)
203 {
204 fNdch[i] = 0;
205 fNrawch[i] = 0;
a4622d0f 206 fNrechits1D[i] = 0;
207 fNrechits3D[i] = 0;
9dda3582 208 }
9b674b76 209
210 fFileName = 0;
fe4da5cc 211}
ddae0931 212
213//___________________________________________
fe4da5cc 214AliRICH::AliRICH(const char *name, const char *title)
ddae0931 215 : AliDetector(name,title)
fe4da5cc 216{
ddae0931 217//Begin_Html
218/*
219 <img src="gif/alirich.gif">
220*/
221//End_Html
222
2e5f0f7b 223 fHits = new TClonesArray("AliRICHHit",1000 );
1cedd08a 224 gAlice->AddHitList(fHits);
b251a2b5 225 fSDigits = new TClonesArray("AliRICHSDigit",100000);
2e5f0f7b 226 fCerenkovs = new TClonesArray("AliRICHCerenkov",1000);
227 gAlice->AddHitList(fCerenkovs);
228 //gAlice->AddHitList(fHits);
b251a2b5 229 fNSDigits = 0;
2e5f0f7b 230 fNcerenkovs = 0;
231 fIshunt = 0;
ddae0931 232
9dda3582 233 //fNdch = new Int_t[kNCH];
ddae0931 234
237c933d 235 fDchambers = new TObjArray(kNCH);
2e5f0f7b 236
a4622d0f 237 fRecHits1D = new TObjArray(kNCH);
238 fRecHits3D = new TObjArray(kNCH);
ddae0931 239
240 Int_t i;
241
237c933d 242 for (i=0; i<kNCH ;i++) {
2e5f0f7b 243 (*fDchambers)[i] = new TClonesArray("AliRICHDigit",10000);
ddae0931 244 fNdch[i]=0;
245 }
2e5f0f7b 246
9dda3582 247 //fNrawch = new Int_t[kNCH];
ddae0931 248
237c933d 249 fRawClusters = new TObjArray(kNCH);
d28b5fc2 250 //printf("Created fRwClusters with adress:%p",fRawClusters);
2e5f0f7b 251
237c933d 252 for (i=0; i<kNCH ;i++) {
2e5f0f7b 253 (*fRawClusters)[i] = new TClonesArray("AliRICHRawCluster",10000);
254 fNrawch[i]=0;
255 }
256
9dda3582 257 //fNrechits = new Int_t[kNCH];
ddae0931 258
237c933d 259 for (i=0; i<kNCH ;i++) {
a4622d0f 260 (*fRecHits1D)[i] = new TClonesArray("AliRICHRecHit1D",1000);
261 }
262 for (i=0; i<kNCH ;i++) {
263 (*fRecHits3D)[i] = new TClonesArray("AliRICHRecHit3D",1000);
2e5f0f7b 264 }
d28b5fc2 265 //printf("Created fRecHits with adress:%p",fRecHits);
2e5f0f7b 266
267
ddae0931 268 SetMarkerColor(kRed);
9b674b76 269
3cc8edee 270 /*fChambers = new TObjArray(kNCH);
9b674b76 271 for (i=0; i<kNCH; i++)
3cc8edee 272 (*fChambers)[i] = new AliRICHChamber();*/
9b674b76 273
274 fFileName = 0;
fe4da5cc 275}
276
237c933d 277AliRICH::AliRICH(const AliRICH& RICH)
278{
279// Copy Constructor
280}
281
282
ddae0931 283//___________________________________________
fe4da5cc 284AliRICH::~AliRICH()
285{
237c933d 286
287// Destructor of RICH manager class
288
ddae0931 289 fIshunt = 0;
290 delete fHits;
b251a2b5 291 delete fSDigits;
ddae0931 292 delete fCerenkovs;
ee2105a4 293
294 //PH Delete TObjArrays
295 if (fChambers) {
296 fChambers->Delete();
297 delete fChambers;
298 }
299 if (fDchambers) {
300 fDchambers->Delete();
301 delete fDchambers;
302 }
303 if (fRawClusters) {
304 fRawClusters->Delete();
305 delete fRawClusters;
306 }
307 if (fRecHits1D) {
308 fRecHits1D->Delete();
309 delete fRecHits1D;
310 }
311 if (fRecHits3D) {
312 fRecHits3D->Delete();
313 delete fRecHits3D;
314 }
315
fe4da5cc 316}
317
34ead2dd 318
319//_____________________________________________________________________________
320Int_t AliRICH::Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
321{
322//
323// Calls the charge disintegration method of the current chamber and adds
324// the simulated cluster to the root treee
325//
326 Int_t clhits[5];
327 Float_t newclust[4][500];
328 Int_t nnew;
329
330//
331// Integrated pulse height on chamber
332
333 clhits[0]=fNhits+1;
334
335 ((AliRICHChamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust, res);
336 Int_t ic=0;
337
338//
339// Add new clusters
340 for (Int_t i=0; i<nnew; i++) {
341 if (Int_t(newclust[0][i]) > 0) {
342 ic++;
343// Cluster Charge
344 clhits[1] = Int_t(newclust[0][i]);
345// Pad: ix
346 clhits[2] = Int_t(newclust[1][i]);
347// Pad: iy
348 clhits[3] = Int_t(newclust[2][i]);
349// Pad: chamber sector
350 clhits[4] = Int_t(newclust[3][i]);
351
352 //printf(" %d %d %d %d %d\n", clhits[0], clhits[1], clhits[2], clhits[3], clhits[4]);
353
354 AddSDigit(clhits);
355 }
356 }
357
358 if (gAlice->TreeS())
359 {
360 gAlice->TreeS()->Fill();
361 gAlice->TreeS()->Write(0,TObject::kOverwrite);
362 //printf("Filled SDigits...\n");
363 }
364
365return nnew;
366}
367//___________________________________________
368void AliRICH::Hits2SDigits()
369{
370
371// Dummy: sdigits are created during transport.
372// Called from alirun.
373
374 int nparticles = gAlice->GetNtrack();
375 cout << "Particles (RICH):" <<nparticles<<endl;
376 if (nparticles > 0) printf("SDigits were already generated.\n");
377
378}
379
380//___________________________________________
381void AliRICH::SDigits2Digits(Int_t nev, Int_t flag)
382{
383
384//
385// Generate digits.
386// Called from macro. Multiple events, more functionality.
387
388 AliRICHChamber* iChamber;
389
390 printf("Generating tresholds...\n");
391
392 for(Int_t i=0;i<7;i++) {
393 iChamber = &(Chamber(i));
394 iChamber->GenerateTresholds();
395 }
396
397 int nparticles = gAlice->GetNtrack();
398 if (nparticles > 0)
399 {
400 if (fMerger) {
401 fMerger->Init();
402 fMerger->Digitise(nev,flag);
403 }
404 }
405 //Digitise(nev,flag);
406}
407//___________________________________________
408void AliRICH::SDigits2Digits()
409{
410
411//
412// Generate digits
413// Called from alirun, single event only.
414
415 AliRICHChamber* iChamber;
416
417 printf("Generating tresholds...\n");
418
419 for(Int_t i=0;i<7;i++) {
420 iChamber = &(Chamber(i));
421 iChamber->GenerateTresholds();
422 }
423
424 int nparticles = gAlice->GetNtrack();
425 cout << "Particles (RICH):" <<nparticles<<endl;
426 if (nparticles > 0)
427 {
428 if (fMerger) {
429 fMerger->Init();
430 fMerger->Digitise(0,0);
431 }
432 }
433}
434//___________________________________________
435void AliRICH::Digits2Reco()
436{
437
438// Generate clusters
439// Called from alirun, single event only.
440
441 int nparticles = gAlice->GetNtrack();
442 cout << "Particles (RICH):" <<nparticles<<endl;
443 if (nparticles > 0) FindClusters(0,0);
444
445}
446
ddae0931 447//___________________________________________
fe4da5cc 448void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits)
449{
237c933d 450
451//
452// Adds a hit to the Hits list
237c933d 453
ddae0931 454 TClonesArray &lhits = *fHits;
2e5f0f7b 455 new(lhits[fNhits++]) AliRICHHit(fIshunt,track,vol,hits);
fe4da5cc 456}
fe4da5cc 457//_____________________________________________________________________________
ddae0931 458void AliRICH::AddCerenkov(Int_t track, Int_t *vol, Float_t *cerenkovs)
459{
237c933d 460
461//
462// Adds a RICH cerenkov hit to the Cerenkov Hits list
463//
464
ddae0931 465 TClonesArray &lcerenkovs = *fCerenkovs;
466 new(lcerenkovs[fNcerenkovs++]) AliRICHCerenkov(fIshunt,track,vol,cerenkovs);
2e5f0f7b 467 //printf ("Done for Cerenkov %d\n\n\n\n",fNcerenkovs);
fe4da5cc 468}
ddae0931 469//___________________________________________
b251a2b5 470void AliRICH::AddSDigit(Int_t *clhits)
ddae0931 471{
237c933d 472
473//
474// Add a RICH pad hit to the list
475//
476
b251a2b5 477 //printf("fsdigits:%p, data: %d\n",fSDigits,clhits[2]);
633e9a38 478 TClonesArray &lSDigits = *fSDigits;
479 new(lSDigits[fNSDigits++]) AliRICHSDigit(clhits);
33984590 480}
fe4da5cc 481//_____________________________________________________________________________
ddae0931 482void AliRICH::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
483{
237c933d 484
485 //
486 // Add a RICH digit to the list
487 //
ddae0931 488
633e9a38 489 //printf("fdigits:%p, data: %d\n",((TClonesArray*)(*fDchambers)[id]),digits[0]);
490 TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]);
491 new(ldigits[fNdch[id]++]) AliRICHDigit(tracks,charges,digits);
fe4da5cc 492}
493
494//_____________________________________________________________________________
2e5f0f7b 495void AliRICH::AddRawCluster(Int_t id, const AliRICHRawCluster& c)
fe4da5cc 496{
ddae0931 497 //
2e5f0f7b 498 // Add a RICH digit to the list
ddae0931 499 //
2e5f0f7b 500
501 TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]);
502 new(lrawcl[fNrawch[id]++]) AliRICHRawCluster(c);
fe4da5cc 503}
504
2e5f0f7b 505//_____________________________________________________________________________
a4622d0f 506void AliRICH::AddRecHit1D(Int_t id, Float_t *rechit, Float_t *photons, Int_t *padsx, Int_t* padsy)
507{
508
509 //
510 // Add a RICH reconstructed hit to the list
511 //
512
513 TClonesArray &lrec1D = *((TClonesArray*)(*fRecHits1D)[id]);
514 new(lrec1D[fNrechits1D[id]++]) AliRICHRecHit1D(id,rechit,photons,padsx,padsy);
515}
516
517//_____________________________________________________________________________
518void AliRICH::AddRecHit3D(Int_t id, Float_t *rechit)
2e5f0f7b 519{
237c933d 520
521 //
522 // Add a RICH reconstructed hit to the list
523 //
fe4da5cc 524
a4622d0f 525 TClonesArray &lrec3D = *((TClonesArray*)(*fRecHits3D)[id]);
526 new(lrec3D[fNrechits3D[id]++]) AliRICHRecHit3D(id,rechit);
2e5f0f7b 527}
fe4da5cc 528
ddae0931 529//___________________________________________
530void AliRICH::BuildGeometry()
531
fe4da5cc 532{
237c933d 533
534 //
535 // Builds a TNode geometry for event display
536 //
53d98323 537 TNode *node, *subnode, *top;
ddae0931 538
53d98323 539 const int kColorRICH = kRed;
ddae0931 540 //
237c933d 541 top=gAlice->GetGeometry()->GetNode("alice");
53d98323 542
543 AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
91975aa2 544 AliRICHSegmentationV0* segmentation;
53d98323 545 AliRICHChamber* iChamber;
a3d71079 546 AliRICHGeometry* geometry;
547
53d98323 548 iChamber = &(pRICH->Chamber(0));
91975aa2 549 segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel(0);
a3d71079 550 geometry=iChamber->GetGeometryModel();
ddae0931 551
552 new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
53d98323 553
c24372d0 554 Float_t padplane_width = segmentation->GetPadPlaneWidth();
555 Float_t padplane_length = segmentation->GetPadPlaneLength();
556
557 //printf("\n\n\n\n\n In BuildGeometry() npx: %d, npy: %d, dpx: %f, dpy:%f\n\n\n\n\n\n",segmentation->Npx(),segmentation->Npy(),segmentation->Dpx(),segmentation->Dpy());
91975aa2 558
c24372d0 559 new TBRIK("PHOTO","PHOTO","void", padplane_width/2,.1,padplane_length/2);
91975aa2 560
c24372d0 561 //printf("\n\n\n\n\n Padplane w: %f l: %f \n\n\n\n\n", padplane_width/2,padplane_length/2);
562 //printf("\n\n\n\n\n Padplane w: %f l: %f \n\n\n\n\n", segmentation->GetPadPlaneWidth(), segmentation->GetPadPlaneLength());
563
a3d71079 564 Float_t offset = 490 + 1.276 - geometry->GetGapThickness()/2; //distance from center of mother volume to methane
565 Float_t deltaphi = 19.5; //phi angle between center of chambers - z direction
566 Float_t deltatheta = 20; //theta angle between center of chambers - x direction
567 Float_t cosphi = TMath::Cos(deltaphi*TMath::Pi()/180);
568 Float_t sinphi = TMath::Sin(deltaphi*TMath::Pi()/180);
569 Float_t costheta = TMath::Cos(deltatheta*TMath::Pi()/180);
570 Float_t sintheta = TMath::Sin(deltatheta*TMath::Pi()/180);
571
572 //printf("\n\n%f %f %f %f %f %f %f\n\n",offset,deltatheta,deltaphi,cosphi,costheta,sinphi,sintheta);
573
574 new TRotMatrix("rot993","rot993",90., 0. , 90. - deltaphi, 90. , deltaphi, -90. );
575 new TRotMatrix("rot994","rot994",90., -deltatheta , 90. , 90.- deltatheta , 0. , 0. );
576 new TRotMatrix("rot995","rot995",90., 0. , 90. , 90. , 0. , 0. );
577 new TRotMatrix("rot996","rot996",90., deltatheta , 90. , 90 + deltatheta , 0. , 0. );
578 new TRotMatrix("rot997","rot997",90., 360. - deltatheta, 108.2 , 90.- deltatheta ,18.2 , 90 - deltatheta);
579 new TRotMatrix("rot998","rot998",90., 0. , 90 + deltaphi , 90. , deltaphi, 90. );
580 new TRotMatrix("rot999","rot999",90., deltatheta , 108.2 , 90.+ deltatheta ,18.2 , 90 + deltatheta);
581
582 Float_t pos1[3]={0. , offset*cosphi , offset*sinphi};
583 Float_t pos2[3]={offset*sintheta , offset*costheta , 0. };
584 Float_t pos3[3]={0. , offset , 0.};
585 Float_t pos4[3]={-offset*sintheta , offset*costheta , 0.};
586 Float_t pos5[3]={offset*sinphi , offset*costheta*cosphi, -offset*sinphi};
587 Float_t pos6[3]={0. , offset*cosphi , -offset*sinphi};
588 Float_t pos7[3]={ -offset*sinphi , offset*costheta*cosphi, -offset*sinphi};
589
91975aa2 590
237c933d 591 top->cd();
a3d71079 592 //Float_t pos1[3]={0,471.8999,165.2599};
2e5f0f7b 593 //Chamber(0).SetChamberTransform(pos1[0],pos1[1],pos1[2],
a3d71079 594 //new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90);
237c933d 595 node = new TNode("RICH1","RICH1","S_RICH",pos1[0],pos1[1],pos1[2],"rot993");
237c933d 596 node->SetLineColor(kColorRICH);
91975aa2 597 node->cd();
a3d71079 598 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 599 subnode->SetLineColor(kGreen);
600 fNodes->Add(subnode);
a3d71079 601 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 602 subnode->SetLineColor(kGreen);
603 fNodes->Add(subnode);
a3d71079 604 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 605 subnode->SetLineColor(kGreen);
606 fNodes->Add(subnode);
a3d71079 607 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 608 subnode->SetLineColor(kGreen);
609 fNodes->Add(subnode);
a3d71079 610 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 611 subnode->SetLineColor(kGreen);
612 fNodes->Add(subnode);
a3d71079 613 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 614 subnode->SetLineColor(kGreen);
615 fNodes->Add(subnode);
237c933d 616 fNodes->Add(node);
53d98323 617
618
91975aa2 619 top->cd();
a3d71079 620 //Float_t pos2[3]={171,470,0};
2e5f0f7b 621 //Chamber(1).SetChamberTransform(pos2[0],pos2[1],pos2[2],
a3d71079 622 //new TRotMatrix("rot994","rot994",90,-20,90,70,0,0);
237c933d 623 node = new TNode("RICH2","RICH2","S_RICH",pos2[0],pos2[1],pos2[2],"rot994");
237c933d 624 node->SetLineColor(kColorRICH);
91975aa2 625 node->cd();
a3d71079 626 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 627 subnode->SetLineColor(kGreen);
628 fNodes->Add(subnode);
a3d71079 629 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 630 subnode->SetLineColor(kGreen);
631 fNodes->Add(subnode);
a3d71079 632 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 633 subnode->SetLineColor(kGreen);
634 fNodes->Add(subnode);
a3d71079 635 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 636 subnode->SetLineColor(kGreen);
637 fNodes->Add(subnode);
a3d71079 638 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 639 subnode->SetLineColor(kGreen);
640 fNodes->Add(subnode);
a3d71079 641 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 642 subnode->SetLineColor(kGreen);
643 fNodes->Add(subnode);
237c933d 644 fNodes->Add(node);
53d98323 645
646
237c933d 647 top->cd();
a3d71079 648 //Float_t pos3[3]={0,500,0};
2e5f0f7b 649 //Chamber(2).SetChamberTransform(pos3[0],pos3[1],pos3[2],
a3d71079 650 //new TRotMatrix("rot995","rot995",90,0,90,90,0,0);
237c933d 651 node = new TNode("RICH3","RICH3","S_RICH",pos3[0],pos3[1],pos3[2],"rot995");
237c933d 652 node->SetLineColor(kColorRICH);
91975aa2 653 node->cd();
a3d71079 654 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 655 subnode->SetLineColor(kGreen);
656 fNodes->Add(subnode);
a3d71079 657 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 658 subnode->SetLineColor(kGreen);
659 fNodes->Add(subnode);
a3d71079 660 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 661 subnode->SetLineColor(kGreen);
662 fNodes->Add(subnode);
a3d71079 663 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 664 subnode->SetLineColor(kGreen);
665 fNodes->Add(subnode);
a3d71079 666 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 667 subnode->SetLineColor(kGreen);
668 fNodes->Add(subnode);
a3d71079 669 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 670 subnode->SetLineColor(kGreen);
671 fNodes->Add(subnode);
237c933d 672 fNodes->Add(node);
53d98323 673
237c933d 674 top->cd();
a3d71079 675 //Float_t pos4[3]={-171,470,0};
2e5f0f7b 676 //Chamber(3).SetChamberTransform(pos4[0],pos4[1],pos4[2],
a3d71079 677 //new TRotMatrix("rot996","rot996",90,20,90,110,0,0);
237c933d 678 node = new TNode("RICH4","RICH4","S_RICH",pos4[0],pos4[1],pos4[2],"rot996");
237c933d 679 node->SetLineColor(kColorRICH);
91975aa2 680 node->cd();
a3d71079 681 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 682 subnode->SetLineColor(kGreen);
683 fNodes->Add(subnode);
a3d71079 684 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 685 subnode->SetLineColor(kGreen);
686 fNodes->Add(subnode);
a3d71079 687 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 688 subnode->SetLineColor(kGreen);
689 fNodes->Add(subnode);
a3d71079 690 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 691 subnode->SetLineColor(kGreen);
692 fNodes->Add(subnode);
a3d71079 693 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 694 subnode->SetLineColor(kGreen);
695 fNodes->Add(subnode);
a3d71079 696 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 697 subnode->SetLineColor(kGreen);
698 fNodes->Add(subnode);
237c933d 699 fNodes->Add(node);
53d98323 700
701
237c933d 702 top->cd();
a3d71079 703 //Float_t pos5[3]={161.3999,443.3999,-165.3};
2e5f0f7b 704 //Chamber(4).SetChamberTransform(pos5[0],pos5[1],pos5[2],
a3d71079 705 //new TRotMatrix("rot997","rot997",90,340,108.1999,70,18.2,70);
237c933d 706 node = new TNode("RICH5","RICH5","S_RICH",pos5[0],pos5[1],pos5[2],"rot997");
237c933d 707 node->SetLineColor(kColorRICH);
91975aa2 708 node->cd();
a3d71079 709 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 710 subnode->SetLineColor(kGreen);
711 fNodes->Add(subnode);
a3d71079 712 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 713 subnode->SetLineColor(kGreen);
714 fNodes->Add(subnode);
a3d71079 715 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 716 subnode->SetLineColor(kGreen);
717 fNodes->Add(subnode);
a3d71079 718 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 719 subnode->SetLineColor(kGreen);
720 fNodes->Add(subnode);
a3d71079 721 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 722 subnode->SetLineColor(kGreen);
723 fNodes->Add(subnode);
a3d71079 724 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 725 subnode->SetLineColor(kGreen);
726 fNodes->Add(subnode);
237c933d 727 fNodes->Add(node);
53d98323 728
729
237c933d 730 top->cd();
a3d71079 731 //Float_t pos6[3]={0., 471.9, -165.3,};
2e5f0f7b 732 //Chamber(5).SetChamberTransform(pos6[0],pos6[1],pos6[2],
a3d71079 733 //new TRotMatrix("rot998","rot998",90,0,109.3099,90,19.30999,90);
237c933d 734 node = new TNode("RICH6","RICH6","S_RICH",pos6[0],pos6[1],pos6[2],"rot998");
237c933d 735 node->SetLineColor(kColorRICH);
91975aa2 736 fNodes->Add(node);node->cd();
a3d71079 737 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 738 subnode->SetLineColor(kGreen);
739 fNodes->Add(subnode);
a3d71079 740 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 741 subnode->SetLineColor(kGreen);
742 fNodes->Add(subnode);
a3d71079 743 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 744 subnode->SetLineColor(kGreen);
745 fNodes->Add(subnode);
a3d71079 746 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 747 subnode->SetLineColor(kGreen);
748 fNodes->Add(subnode);
a3d71079 749 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 750 subnode->SetLineColor(kGreen);
751 fNodes->Add(subnode);
a3d71079 752 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 753 subnode->SetLineColor(kGreen);
754 fNodes->Add(subnode);
53d98323 755
756
237c933d 757 top->cd();
a3d71079 758 //Float_t pos7[3]={-161.399,443.3999,-165.3};
2e5f0f7b 759 //Chamber(6).SetChamberTransform(pos7[0],pos7[1],pos7[2],
a3d71079 760 //new TRotMatrix("rot999","rot999",90,20,108.1999,110,18.2,110);
237c933d 761 node = new TNode("RICH7","RICH7","S_RICH",pos7[0],pos7[1],pos7[2],"rot999");
762 node->SetLineColor(kColorRICH);
91975aa2 763 node->cd();
a3d71079 764 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 765 subnode->SetLineColor(kGreen);
766 fNodes->Add(subnode);
a3d71079 767 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 768 subnode->SetLineColor(kGreen);
769 fNodes->Add(subnode);
a3d71079 770 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
91975aa2 771 subnode->SetLineColor(kGreen);
772 fNodes->Add(subnode);
a3d71079 773 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 774 subnode->SetLineColor(kGreen);
775 fNodes->Add(subnode);
a3d71079 776 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 777 subnode->SetLineColor(kGreen);
778 fNodes->Add(subnode);
a3d71079 779 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
91975aa2 780 subnode->SetLineColor(kGreen);
781 fNodes->Add(subnode);
237c933d 782 fNodes->Add(node);
ddae0931 783
fe4da5cc 784}
785
452a64c6 786//___________________________________________
787void AliRICH::CreateGeometry()
788{
789 //
790 // Create the geometry for RICH version 1
791 //
792 // Modified by: N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
793 // R.A. Fini (INFN - BARI, Rosanna.Fini@ba.infn.it)
794 // R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
795 //
796 //Begin_Html
797 /*
798 <img src="picts/AliRICHv1.gif">
799 */
800 //End_Html
801 //Begin_Html
802 /*
803 <img src="picts/AliRICHv1Tree.gif">
804 */
805 //End_Html
806
807 AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
c24372d0 808 AliRICHSegmentationV0* segmentation;
452a64c6 809 AliRICHGeometry* geometry;
810 AliRICHChamber* iChamber;
811
812 iChamber = &(pRICH->Chamber(0));
c24372d0 813 segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel(0);
452a64c6 814 geometry=iChamber->GetGeometryModel();
815
816 Float_t distance;
817 distance = geometry->GetFreonThickness()/2 + geometry->GetQuartzThickness() + geometry->GetGapThickness();
818 geometry->SetRadiatorToPads(distance);
819
3587e146 820 //Opaque quartz thickness
e293afae 821 Float_t oqua_thickness = .5;
683b273b 822 //CsI dimensions
823
53d98323 824 //Float_t csi_length = 160*.8 + 2.6;
825 //Float_t csi_width = 144*.84 + 2*2.6;
826
ca96c9ea 827 Float_t csi_width = segmentation->Npx()*segmentation->Dpx() + segmentation->DeadZone();
828 Float_t csi_length = segmentation->Npy()*segmentation->Dpy() + 2*segmentation->DeadZone();
c24372d0 829
830 //printf("\n\n\n\n\n In CreateGeometry() npx: %d, npy: %d, dpx: %f, dpy:%f deadzone: %f \n\n\n\n\n\n",segmentation->Npx(),segmentation->Npy(),segmentation->Dpx(),segmentation->Dpy(),segmentation->DeadZone());
53d98323 831
832 Int_t *idtmed = fIdtmed->GetArray()-999;
452a64c6 833
834 Int_t i;
835 Float_t zs;
836 Int_t idrotm[1099];
837 Float_t par[3];
838
839 // --- Define the RICH detector
840 // External aluminium box
e293afae 841 par[0] = 68.8;
6197ac87 842 par[1] = 13; //Original Settings
e293afae 843 par[2] = 70.86;
452a64c6 844 /*par[0] = 73.15;
845 par[1] = 11.5;
846 par[2] = 71.1;*/
847 gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
848
683b273b 849 // Air
850 par[0] = 66.3;
6197ac87 851 par[1] = 13; //Original Settings
683b273b 852 par[2] = 68.35;
452a64c6 853 /*par[0] = 66.55;
854 par[1] = 11.5;
855 par[2] = 64.8;*/
856 gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3);
857
6197ac87 858 // Air 2 (cutting the lower part of the box)
859
860 par[0] = 1.25;
861 par[1] = 3; //Original Settings
862 par[2] = 70.86;
863 gMC->Gsvolu("AIR2", "BOX ", idtmed[1000], par, 3);
864
865 // Air 3 (cutting the lower part of the box)
866
867 par[0] = 66.3;
868 par[1] = 3; //Original Settings
869 par[2] = 1.2505;
870 gMC->Gsvolu("AIR3", "BOX ", idtmed[1000], par, 3);
871
452a64c6 872 // Honeycomb
e293afae 873 par[0] = 66.3;
452a64c6 874 par[1] = .188; //Original Settings
e293afae 875 par[2] = 68.35;
452a64c6 876 /*par[0] = 66.55;
877 par[1] = .188;
878 par[2] = 63.1;*/
879 gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
880
881 // Aluminium sheet
e293afae 882 par[0] = 66.3;
452a64c6 883 par[1] = .025; //Original Settings
e293afae 884 par[2] = 68.35;
452a64c6 885 /*par[0] = 66.5;
886 par[1] = .025;
887 par[2] = 63.1;*/
888 gMC->Gsvolu("ALUM", "BOX ", idtmed[1009], par, 3);
889
890 // Quartz
891 par[0] = geometry->GetQuartzWidth()/2;
892 par[1] = geometry->GetQuartzThickness()/2;
893 par[2] = geometry->GetQuartzLength()/2;
894 /*par[0] = 63.1;
895 par[1] = .25; //Original Settings
896 par[2] = 65.5;*/
897 /*par[0] = geometry->GetQuartzWidth()/2;
898 par[1] = geometry->GetQuartzThickness()/2;
899 par[2] = geometry->GetQuartzLength()/2;*/
900 //printf("\n\n\n\n\n\n\n\\n\n\n\n Gap Thickness: %f %f %f\n\n\n\n\n\n\n\n\n\n\n\n\n\n",par[0],par[1],par[2]);
901 gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
902
903 // Spacers (cylinders)
904 par[0] = 0.;
905 par[1] = .5;
906 par[2] = geometry->GetFreonThickness()/2;
907 gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
908
683b273b 909 // Feet (freon slabs supports)
910
911 par[0] = .7;
912 par[1] = .3;
913 par[2] = 1.9;
914 gMC->Gsvolu("FOOT", "BOX", idtmed[1009], par, 3);
915
452a64c6 916 // Opaque quartz
e293afae 917 par[0] = geometry->GetQuartzWidth()/2;
918 par[1] = .2;
919 par[2] = geometry->GetQuartzLength()/2;
920 /*par[0] = 61.95;
452a64c6 921 par[1] = .2; //Original Settings
e293afae 922 par[2] = 66.5;*/
452a64c6 923 /*par[0] = 66.5;
924 par[1] = .2;
925 par[2] = 61.95;*/
926 gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
927
928 // Frame of opaque quartz
53a60579 929 par[0] = geometry->GetOuterFreonWidth()/2;
930 //+ oqua_thickness;
452a64c6 931 par[1] = geometry->GetFreonThickness()/2;
53a60579 932 par[2] = geometry->GetOuterFreonLength()/2;
933 //+ oqua_thickness;
452a64c6 934 /*par[0] = 20.65;
935 par[1] = .5; //Original Settings
936 par[2] = 66.5;*/
937 /*par[0] = 66.5;
938 par[1] = .5;
939 par[2] = 20.65;*/
940 gMC->Gsvolu("OQF1", "BOX ", idtmed[1007], par, 3);
941
e293afae 942 par[0] = geometry->GetInnerFreonWidth()/2;
452a64c6 943 par[1] = geometry->GetFreonThickness()/2;
e293afae 944 par[2] = geometry->GetInnerFreonLength()/2;
452a64c6 945 gMC->Gsvolu("OQF2", "BOX ", idtmed[1007], par, 3);
946
947 // Little bar of opaque quartz
e293afae 948 //par[0] = .275;
949 //par[1] = geometry->GetQuartzThickness()/2;
53a60579 950 //par[2] = geometry->GetInnerFreonLength()/2 - 2.4;
e293afae 951 //par[2] = geometry->GetInnerFreonLength()/2;
53a60579 952 //+ oqua_thickness;
452a64c6 953 /*par[0] = .275;
954 par[1] = .25; //Original Settings
955 par[2] = 63.1;*/
956 /*par[0] = 63.1;
957 par[1] = .25;
958 par[2] = .275;*/
e293afae 959 //gMC->Gsvolu("BARR", "BOX ", idtmed[1007], par, 3);
452a64c6 960
961 // Freon
e293afae 962 par[0] = geometry->GetOuterFreonWidth()/2 - oqua_thickness;
452a64c6 963 par[1] = geometry->GetFreonThickness()/2;
e293afae 964 par[2] = geometry->GetOuterFreonLength()/2 - 2*oqua_thickness;
452a64c6 965 /*par[0] = 20.15;
966 par[1] = .5; //Original Settings
967 par[2] = 65.5;*/
968 /*par[0] = 65.5;
969 par[1] = .5;
970 par[2] = 20.15;*/
971 gMC->Gsvolu("FRE1", "BOX ", idtmed[1003], par, 3);
972
e293afae 973 par[0] = geometry->GetInnerFreonWidth()/2 - oqua_thickness;
452a64c6 974 par[1] = geometry->GetFreonThickness()/2;
e293afae 975 par[2] = geometry->GetInnerFreonLength()/2 - 2*oqua_thickness;
452a64c6 976 gMC->Gsvolu("FRE2", "BOX ", idtmed[1003], par, 3);
977
978 // Methane
e293afae 979 //par[0] = 64.8;
683b273b 980 par[0] = csi_width/2;
452a64c6 981 par[1] = geometry->GetGapThickness()/2;
982 //printf("\n\n\n\n\n\n\n\\n\n\n\n Gap Thickness: %f\n\n\n\n\n\n\n\n\n\n\n\n\n\n",par[1]);
e293afae 983 //par[2] = 64.8;
683b273b 984 par[2] = csi_length/2;
452a64c6 985 gMC->Gsvolu("META", "BOX ", idtmed[1004], par, 3);
986
987 // Methane gap
e293afae 988 //par[0] = 64.8;
683b273b 989 par[0] = csi_width/2;
452a64c6 990 par[1] = geometry->GetProximityGapThickness()/2;
991 //printf("\n\n\n\n\n\n\n\\n\n\n\n Gap Thickness: %f\n\n\n\n\n\n\n\n\n\n\n\n\n\n",par[1]);
e293afae 992 //par[2] = 64.8;
683b273b 993 par[2] = csi_length/2;
452a64c6 994 gMC->Gsvolu("GAP ", "BOX ", idtmed[1008], par, 3);
995
996 // CsI photocathode
e293afae 997 //par[0] = 64.8;
683b273b 998 par[0] = csi_width/2;
452a64c6 999 par[1] = .25;
e293afae 1000 //par[2] = 64.8;
683b273b 1001 par[2] = csi_length/2;
452a64c6 1002 gMC->Gsvolu("CSI ", "BOX ", idtmed[1005], par, 3);
1003
1004 // Anode grid
1005 par[0] = 0.;
1006 par[1] = .001;
1007 par[2] = 20.;
1008 gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
683b273b 1009
6197ac87 1010 // Wire supports
1011 // Bar of metal
1012
1013 par[0] = csi_width/2;
1014 par[1] = 1.05;
1015 par[2] = 1.05;
1016 gMC->Gsvolu("WSMe", "BOX ", idtmed[1009], par, 3);
1017
1018 // Ceramic pick up (base)
1019
1020 par[0] = csi_width/2;
1021 par[1] = .25;
1022 par[2] = 1.05;
1023 gMC->Gsvolu("WSG1", "BOX ", idtmed[1010], par, 3);
1024
1025 // Ceramic pick up (head)
1026
1027 par[0] = csi_width/2;
1028 par[1] = .1;
1029 par[2] = .1;
1030 gMC->Gsvolu("WSG2", "BOX ", idtmed[1010], par, 3);
1031
683b273b 1032 // Aluminium supports for methane and CsI
1033 // Short bar
1034
1035 par[0] = csi_width/2;
1036 par[1] = geometry->GetGapThickness()/2 + .25;
1037 par[2] = (68.35 - csi_length/2)/2;
1038 gMC->Gsvolu("SMSH", "BOX", idtmed[1009], par, 3);
452a64c6 1039
683b273b 1040 // Long bar
1041
1042 par[0] = (66.3 - csi_width/2)/2;
1043 par[1] = geometry->GetGapThickness()/2 + .25;
1044 par[2] = csi_length/2 + 68.35 - csi_length/2;
1045 gMC->Gsvolu("SMLG", "BOX", idtmed[1009], par, 3);
1046
1047 // Aluminium supports for freon
1048 // Short bar
1049
1050 par[0] = geometry->GetQuartzWidth()/2;
1051 par[1] = .3;
1052 par[2] = (68.35 - geometry->GetQuartzLength()/2)/2;
1053 gMC->Gsvolu("SFSH", "BOX", idtmed[1009], par, 3);
1054
1055 // Long bar
1056
1057 par[0] = (66.3 - geometry->GetQuartzWidth()/2)/2;
1058 par[1] = .3;
1059 par[2] = geometry->GetQuartzLength()/2 + 68.35 - geometry->GetQuartzLength()/2;
1060 gMC->Gsvolu("SFLG", "BOX", idtmed[1009], par, 3);
1061
6197ac87 1062 // PCB backplane
683b273b 1063
6197ac87 1064 par[0] = csi_width/2;
1065 par[1] = .25;
1066 par[2] = csi_length/4 -.5025;
1067 gMC->Gsvolu("PCB ", "BOX", idtmed[1011], par, 3);
683b273b 1068
6197ac87 1069
1070 // Backplane supports
1071
1072 // Aluminium slab
1073
1074 par[0] = 33.15;
1075 par[1] = 2;
1076 par[2] = 21.65;
1077 gMC->Gsvolu("BACK", "BOX", idtmed[1009], par, 3);
1078
1079 // Big hole
1080
1081 par[0] = 9.05;
1082 par[1] = 2;
1083 par[2] = 4.4625;
1084 gMC->Gsvolu("BKHL", "BOX", idtmed[1000], par, 3);
1085
1086 // Small hole
1087
1088 par[0] = 5.7;
1089 par[1] = 2;
1090 par[2] = 4.4625;
1091 gMC->Gsvolu("BKHS", "BOX", idtmed[1000], par, 3);
1092
1093 // Place holes inside backplane support
1094
1095 gMC->Gspos("BKHS", 1, "BACK", .8 + 5.7,0., .6 + 4.4625, 0, "ONLY");
1096 gMC->Gspos("BKHS", 2, "BACK", -.8 - 5.7,0., .6 + 4.4625, 0, "ONLY");
1097 gMC->Gspos("BKHS", 3, "BACK", .8 + 5.7,0., -.6 - 4.4625, 0, "ONLY");
1098 gMC->Gspos("BKHS", 4, "BACK", -.8 - 5.7,0., -.6 - 4.4625, 0, "ONLY");
1099 gMC->Gspos("BKHS", 5, "BACK", .8 + 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1100 gMC->Gspos("BKHS", 6, "BACK", -.8 - 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1101 gMC->Gspos("BKHS", 7, "BACK", .8 + 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1102 gMC->Gspos("BKHS", 8, "BACK", -.8 - 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1103 gMC->Gspos("BKHL", 1, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., .6 + 4.4625, 0, "ONLY");
1104 gMC->Gspos("BKHL", 2, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 4.4625, 0, "ONLY");
1105 gMC->Gspos("BKHL", 3, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 4.4625, 0, "ONLY");
1106 gMC->Gspos("BKHL", 4, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 4.4625, 0, "ONLY");
1107 gMC->Gspos("BKHL", 5, "BACK", .8 + 11.4+ 1.6 + 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1108 gMC->Gspos("BKHL", 6, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1109 gMC->Gspos("BKHL", 7, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1110 gMC->Gspos("BKHL", 8, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1111
1112
1113
452a64c6 1114 // --- Places the detectors defined with GSVOLU
1115 // Place material inside RICH
6197ac87 1116 gMC->Gspos("SRIC", 1, "RICH", 0.,0., 0., 0, "ONLY");
4591f067 1117 gMC->Gspos("AIR2", 1, "RICH", 66.3 + 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 0., 0, "ONLY");
1118 gMC->Gspos("AIR2", 2, "RICH", -66.3 - 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 0., 0, "ONLY");
1119 gMC->Gspos("AIR3", 1, "RICH", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, -68.35 - 1.25, 0, "ONLY");
1120 gMC->Gspos("AIR3", 2, "RICH", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 68.35 + 1.25, 0, "ONLY");
6197ac87 1121
683b273b 1122
1123 gMC->Gspos("ALUM", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.025, 0., 0, "ONLY");
1124 gMC->Gspos("HONE", 1, "SRIC", 0., 1.276- geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .188, 0., 0, "ONLY");
1125 gMC->Gspos("ALUM", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .025, 0., 0, "ONLY");
1126 gMC->Gspos("FOOT", 1, "SRIC", 64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
1127 gMC->Gspos("FOOT", 2, "SRIC", 21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3 , 36.9, 0, "ONLY");
1128 gMC->Gspos("FOOT", 3, "SRIC", -21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
1129 gMC->Gspos("FOOT", 4, "SRIC", -64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
1130 gMC->Gspos("FOOT", 5, "SRIC", 64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1131 gMC->Gspos("FOOT", 6, "SRIC", 21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1132 gMC->Gspos("FOOT", 7, "SRIC", -21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1133 gMC->Gspos("FOOT", 8, "SRIC", -64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
452a64c6 1134 gMC->Gspos("OQUA", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .2, 0., 0, "ONLY");
1135
683b273b 1136 // Supports placing
1137
1138 // Methane supports
1139 gMC->Gspos("SMLG", 1, "SRIC", csi_width/2 + (66.3 - csi_width/2)/2, 1.276 + .25, 0., 0, "ONLY");
1140 gMC->Gspos("SMLG", 2, "SRIC", - csi_width/2 - (66.3 - csi_width/2)/2, 1.276 + .25, 0., 0, "ONLY");
1141 gMC->Gspos("SMSH", 1, "SRIC", 0., 1.276 + .25, csi_length/2 + (68.35 - csi_length/2)/2, 0, "ONLY");
1142 gMC->Gspos("SMSH", 2, "SRIC", 0., 1.276 + .25, - csi_length/2 - (68.35 - csi_length/2)/2, 0, "ONLY");
1143
1144 //Freon supports
1145
1146 Float_t supp_y = 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness() - .2 + .3; //y position of freon supports
1147
1148 gMC->Gspos("SFLG", 1, "SRIC", geometry->GetQuartzWidth()/2 + (66.3 - geometry->GetQuartzWidth()/2)/2, supp_y, 0., 0, "ONLY");
1149 gMC->Gspos("SFLG", 2, "SRIC", - geometry->GetQuartzWidth()/2 - (66.3 - geometry->GetQuartzWidth()/2)/2, supp_y, 0., 0, "ONLY");
1150 gMC->Gspos("SFSH", 1, "SRIC", 0., supp_y, geometry->GetQuartzLength()/2 + (68.35 - geometry->GetQuartzLength()/2)/2, 0, "ONLY");
1151 gMC->Gspos("SFSH", 2, "SRIC", 0., supp_y, - geometry->GetQuartzLength()/2 - (68.35 - geometry->GetQuartzLength()/2)/2, 0, "ONLY");
1152
452a64c6 1153 AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
1154
c24372d0 1155 //Placing of the spacers inside the freon slabs
1156
1157 Int_t nspacers = 30;
452a64c6 1158 //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n Spacers:%d\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n",nspacers);
1159
1160 //printf("Nspacers: %d", nspacers);
1161
c24372d0 1162 for (i = 0; i < nspacers/3; i++) {
1163 zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2;
1164 gMC->Gspos("SPAC", i, "FRE1", 10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
452a64c6 1165 }
c24372d0 1166
1167 for (i = nspacers/3; i < (nspacers*2)/3; i++) {
1168 zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1169 gMC->Gspos("SPAC", i, "FRE1", 0, 0., zs, idrotm[1019], "ONLY"); //Original settings
1170 }
1171
1172 for (i = (nspacers*2)/3; i < nspacers; ++i) {
1173 zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1174 gMC->Gspos("SPAC", i, "FRE1", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
452a64c6 1175 }
1176
c24372d0 1177 for (i = 0; i < nspacers/3; i++) {
1178 zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2;
1179 gMC->Gspos("SPAC", i, "FRE2", 10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
452a64c6 1180 }
c24372d0 1181
1182 for (i = nspacers/3; i < (nspacers*2)/3; i++) {
1183 zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1184 gMC->Gspos("SPAC", i, "FRE2", 0, 0., zs, idrotm[1019], "ONLY"); //Original settings
452a64c6 1185 }
1186
c24372d0 1187 for (i = (nspacers*2)/3; i < nspacers; ++i) {
1188 zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1189 gMC->Gspos("SPAC", i, "FRE2", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
1190 }
452a64c6 1191
c24372d0 1192
452a64c6 1193 gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY");
1194 gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY");
683b273b 1195 gMC->Gspos("OQF1", 1, "SRIC", geometry->GetOuterFreonWidth()/2 + geometry->GetInnerFreonWidth()/2 + 2, 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness()/2, 0., 0, "ONLY"); //Original settings (31.3)
633e9a38 1196// printf("Opaque quartz in SRIC %f\n", 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness()/2);
452a64c6 1197 gMC->Gspos("OQF2", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY"); //Original settings
683b273b 1198 gMC->Gspos("OQF1", 3, "SRIC", - (geometry->GetOuterFreonWidth()/2 + geometry->GetInnerFreonWidth()/2) - 2, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY"); //Original settings (-31.3)
e293afae 1199 //gMC->Gspos("BARR", 1, "QUAR", - geometry->GetInnerFreonWidth()/2 - oqua_thickness, 0., 0., 0, "ONLY"); //Original settings (-21.65)
1200 //gMC->Gspos("BARR", 2, "QUAR", geometry->GetInnerFreonWidth()/2 + oqua_thickness, 0., 0., 0, "ONLY"); //Original settings (21.65)
452a64c6 1201 gMC->Gspos("QUAR", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness()/2, 0., 0, "ONLY");
1202 gMC->Gspos("GAP ", 1, "META", 0., geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 - 0.0001, 0., 0, "ONLY");
1203 gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
1204 gMC->Gspos("CSI ", 1, "SRIC", 0., 1.276 + geometry->GetGapThickness()/2 + .25, 0., 0, "ONLY");
a3d71079 1205 printf("CSI pos: %f\n",1.276 + geometry->GetGapThickness()/2 + .25);
6197ac87 1206
1207 // Wire support placing
1208
1209 gMC->Gspos("WSG2", 1, "GAP ", 0., geometry->GetProximityGapThickness()/2 - .1, 0., 0, "ONLY");
1210 gMC->Gspos("WSG1", 1, "CSI ", 0., 0., 0., 0, "ONLY");
1211 gMC->Gspos("WSMe", 1, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, 0., 0, "ONLY");
1212
1213 // Backplane placing
1214
1215 gMC->Gspos("BACK", 1, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 43.3, 0, "ONLY");
1216 gMC->Gspos("BACK", 2, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2 , 43.3, 0, "ONLY");
1217 gMC->Gspos("BACK", 3, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
1218 gMC->Gspos("BACK", 4, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
1219 gMC->Gspos("BACK", 5, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
1220 gMC->Gspos("BACK", 6, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
1221
1222 // PCB placing
1223
4591f067 1224 gMC->Gspos("PCB ", 1, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, csi_width/4 + .5025 + 2.5, 0, "ONLY");
1225 gMC->Gspos("PCB ", 2, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, -csi_width/4 - .5025 - 2.5, 0, "ONLY");
6197ac87 1226
1227
452a64c6 1228
1229 //printf("Position of the gap: %f to %f\n", 1.276 + geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 - .2, 1.276 + geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 + .2);
1230
1231 // Place RICH inside ALICE apparatus
c24372d0 1232
a3d71079 1233 /* old values
1234
1235 AliMatrix(idrotm[1000], 90., 0., 70.69, 90., 19.31, -90.);
1236 AliMatrix(idrotm[1001], 90., -20., 90., 70., 0., 0.);
1237 AliMatrix(idrotm[1002], 90., 0., 90., 90., 0., 0.);
1238 AliMatrix(idrotm[1003], 90., 20., 90., 110., 0., 0.);
1239 AliMatrix(idrotm[1004], 90., 340., 108.2, 70., 18.2, 70.);
1240 AliMatrix(idrotm[1005], 90., 0., 109.31, 90., 19.31, 90.);
1241 AliMatrix(idrotm[1006], 90., 20., 108.2, 110., 18.2, 110.);
1242
1243 gMC->Gspos("RICH", 1, "ALIC", 0., 471.9, 165.26, idrotm[1000], "ONLY");
1244 gMC->Gspos("RICH", 2, "ALIC", 171., 470., 0., idrotm[1001], "ONLY");
1245 gMC->Gspos("RICH", 3, "ALIC", 0., 500., 0., idrotm[1002], "ONLY");
1246 gMC->Gspos("RICH", 4, "ALIC", -171., 470., 0., idrotm[1003], "ONLY");
1247 gMC->Gspos("RICH", 5, "ALIC", 161.4, 443.4, -165.3, idrotm[1004], "ONLY");
1248 gMC->Gspos("RICH", 6, "ALIC", 0., 471.9, -165.3, idrotm[1005], "ONLY");
1249 gMC->Gspos("RICH", 7, "ALIC", -161.4, 443.4, -165.3, idrotm[1006], "ONLY");*/
1250
1251 // The placing of the chambers is measured from the vertex to the base of the methane vessel (490 cm)
1252
1253 Float_t offset = 490 + 1.276 - geometry->GetGapThickness()/2; //distance from center of mother volume to methane
1254 Float_t deltaphi = 19.5; //phi angle between center of chambers - z direction
1255 Float_t deltatheta = 20; //theta angle between center of chambers - x direction
1256 Float_t cosphi = TMath::Cos(deltaphi*TMath::Pi()/180);
1257 Float_t sinphi = TMath::Sin(deltaphi*TMath::Pi()/180);
1258 Float_t costheta = TMath::Cos(deltatheta*TMath::Pi()/180);
1259 Float_t sintheta = TMath::Sin(deltatheta*TMath::Pi()/180);
1260
1261 //printf("\n\n%f %f %f %f %f %f %f\n\n",offset,deltatheta,deltaphi,cosphi,costheta,sinphi,sintheta);
1262
1263 AliMatrix(idrotm[1000], 90., 0. , 90. - deltaphi, 90. , deltaphi, -90. );
1264 AliMatrix(idrotm[1001], 90., -deltatheta , 90. , 90.- deltatheta , 0. , 0. );
1265 AliMatrix(idrotm[1002], 90., 0. , 90. , 90. , 0. , 0. );
1266 AliMatrix(idrotm[1003], 90., deltatheta , 90. , 90 + deltatheta , 0. , 0. );
1267 AliMatrix(idrotm[1004], 90., 360. - deltatheta, 108.2 , 90.- deltatheta ,18.2 , 90 - deltatheta);
1268 AliMatrix(idrotm[1005], 90., 0. , 90 + deltaphi , 90. , deltaphi, 90. );
1269 AliMatrix(idrotm[1006], 90., deltatheta , 108.2 , 90.+ deltatheta ,18.2 , 90 + deltatheta);
1270
1271 gMC->Gspos("RICH", 1, "ALIC", 0. , offset*cosphi , offset*sinphi ,idrotm[1000], "ONLY");
1272 gMC->Gspos("RICH", 2, "ALIC", (offset)*sintheta , offset*costheta , 0. ,idrotm[1001], "ONLY");
1273 gMC->Gspos("RICH", 3, "ALIC", 0. , offset , 0. ,idrotm[1002], "ONLY");
1274 gMC->Gspos("RICH", 4, "ALIC", -(offset)*sintheta, offset*costheta , 0. ,idrotm[1003], "ONLY");
1275 gMC->Gspos("RICH", 5, "ALIC", (offset)*sinphi , offset*costheta*cosphi, -offset*sinphi,idrotm[1004], "ONLY");
1276 gMC->Gspos("RICH", 6, "ALIC", 0. , offset*cosphi , -offset*sinphi,idrotm[1005], "ONLY");
1277 gMC->Gspos("RICH", 7, "ALIC", -(offset)*sinphi , offset*costheta*cosphi, -offset*sinphi,idrotm[1006], "ONLY");
452a64c6 1278
1279}
1280
1281
1282//___________________________________________
1283void AliRICH::CreateMaterials()
1284{
1285 //
1286 // *** DEFINITION OF AVAILABLE RICH MATERIALS ***
1287 // ORIGIN : NICK VAN EIJNDHOVEN
1288 // Modified by: N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
1289 // R.A. Fini (INFN - BARI, Rosanna.Fini@ba.infn.it)
1290 // R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
1291 //
1292 Int_t isxfld = gAlice->Field()->Integ();
1293 Float_t sxmgmx = gAlice->Field()->Max();
1294 Int_t i;
1295
1296 /************************************Antonnelo's Values (14-vectors)*****************************************/
1297 /*
1298 Float_t ppckov[14] = { 5.63e-9,5.77e-9,5.9e-9,6.05e-9,6.2e-9,6.36e-9,6.52e-9,
1299 6.7e-9,6.88e-9,7.08e-9,7.3e-9,7.51e-9,7.74e-9,8e-9 };
1300 Float_t rIndexQuarz[14] = { 1.528309,1.533333,
1301 1.538243,1.544223,1.550568,1.55777,
1302 1.565463,1.574765,1.584831,1.597027,
1303 1.611858,1.6277,1.6472,1.6724 };
1304 Float_t rIndexOpaqueQuarz[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1305 Float_t rIndexMethane[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1306 Float_t rIndexGrid[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1307 Float_t abscoFreon[14] = { 179.0987,179.0987,
1308 179.0987,179.0987,179.0987,142.92,56.65,13.95,10.43,7.07,2.03,.5773,.33496,0. };
1309 //Float_t abscoFreon[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
1310 // 1e-5,1e-5,1e-5,1e-5,1e-5 };
1311 Float_t abscoQuarz[14] = { 64.035,39.98,35.665,31.262,27.527,22.815,21.04,17.52,
1312 14.177,9.282,4.0925,1.149,.3627,.10857 };
1313 Float_t abscoOpaqueQuarz[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
1314 1e-5,1e-5,1e-5,1e-5,1e-5 };
1315 Float_t abscoCsI[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
1316 1e-4,1e-4,1e-4,1e-4 };
1317 Float_t abscoMethane[14] = { 1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,
1318 1e6,1e6,1e6 };
1319 Float_t abscoGrid[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
1320 1e-4,1e-4,1e-4,1e-4 };
1321 Float_t efficAll[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1322 Float_t efficCsI[14] = { 6e-4,.005,.0075,.01125,.045,.117,.135,.16575,
1323 .17425,.1785,.1836,.1904,.1938,.221 };
1324 Float_t efficGrid[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1325 */
1326
1327
1328 /**********************************End of Antonnelo's Values**********************************/
1329
1330 /**********************************Values from rich_media.f (31-vectors)**********************************/
1331
1332
1333 //Photons energy intervals
1334 Float_t ppckov[26];
1335 for (i=0;i<26;i++)
1336 {
1337 ppckov[i] = (Float_t(i)*0.1+5.5)*1e-9;
1338 //printf ("Energy intervals: %e\n",ppckov[i]);
1339 }
1340
1341
1342 //Refraction index for quarz
1343 Float_t rIndexQuarz[26];
1344 Float_t e1= 10.666;
1345 Float_t e2= 18.125;
1346 Float_t f1= 46.411;
1347 Float_t f2= 228.71;
1348 for (i=0;i<26;i++)
1349 {
1350 Float_t ene=ppckov[i]*1e9;
1351 Float_t a=f1/(e1*e1 - ene*ene);
1352 Float_t b=f2/(e2*e2 - ene*ene);
1353 rIndexQuarz[i] = TMath::Sqrt(1. + a + b );
1354 //printf ("rIndexQuarz: %e\n",rIndexQuarz[i]);
1355 }
1356
1357 //Refraction index for opaque quarz, methane and grid
1358 Float_t rIndexOpaqueQuarz[26];
1359 Float_t rIndexMethane[26];
1360 Float_t rIndexGrid[26];
1361 for (i=0;i<26;i++)
1362 {
1363 rIndexOpaqueQuarz[i]=1;
1364 rIndexMethane[i]=1.000444;
1365 rIndexGrid[i]=1;
1366 //printf ("rIndexOpaqueQuarz , etc: %e, %e, %e\n",rIndexOpaqueQuarz[i], rIndexMethane[i], rIndexGrid[i]=1);
1367 }
1368
1369 //Absorption index for freon
1370 Float_t abscoFreon[26] = {179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987,
1371 179.0987, 142.9206, 56.64957, 25.58622, 13.95293, 12.03905, 10.42953, 8.804196,
1372 7.069031, 4.461292, 2.028366, 1.293013, .577267, .40746, .334964, 0., 0., 0.};
1373
1374 //Absorption index for quarz
1375 /*Float_t Qzt [21] = {.0,.0,.005,.04,.35,.647,.769,.808,.829,.844,.853,.858,.869,.887,.903,.902,.902,
1376 .906,.907,.907,.907};
1377 Float_t Wavl2[] = {150.,155.,160.0,165.0,170.0,175.0,180.0,185.0,190.0,195.0,200.0,205.0,210.0,
1378 215.0,220.0,225.0,230.0,235.0,240.0,245.0,250.0};
1379 Float_t abscoQuarz[31];
1380 for (Int_t i=0;i<31;i++)
1381 {
1382 Float_t Xlam = 1237.79 / (ppckov[i]*1e9);
1383 if (Xlam <= 160) abscoQuarz[i] = 0;
1384 if (Xlam > 250) abscoQuarz[i] = 1;
1385 else
1386 {
1387 for (Int_t j=0;j<21;j++)
1388 {
1389 //printf ("Passed\n");
1390 if (Xlam > Wavl2[j] && Xlam < Wavl2[j+1])
1391 {
1392 Float_t Dabs = (Qzt[j+1] - Qzt[j])/(Wavl2[j+1] - Wavl2[j]);
1393 Float_t Abso = Qzt[j] + Dabs*(Xlam - Wavl2[j]);
1394 abscoQuarz[i] = -5.0/(TMath::Log(Abso));
1395 }
1396 }
1397 }
1398 printf ("abscoQuarz: %e abscoFreon: %e for energy: %e\n",abscoQuarz[i],abscoFreon[i],ppckov[i]);
1399 }*/
1400
1401 /*Float_t abscoQuarz[31] = {49.64211, 48.41296, 47.46989, 46.50492, 45.13682, 44.47883, 43.1929 , 41.30922, 40.5943 ,
1402 39.82956, 38.98623, 38.6247 , 38.43448, 37.41084, 36.22575, 33.74852, 30.73901, 24.25086,
1403 17.94531, 11.88753, 5.99128, 3.83503, 2.36661, 1.53155, 1.30582, 1.08574, .8779708,
1404 .675275, 0., 0., 0.};
1405
1406 for (Int_t i=0;i<31;i++)
1407 {
1408 abscoQuarz[i] = abscoQuarz[i]/10;
1409 }*/
1410
1411 Float_t abscoQuarz [26] = {105.8, 65.52, 48.58, 42.85, 35.79, 31.262, 28.598, 27.527, 25.007, 22.815, 21.004,
1412 19.266, 17.525, 15.878, 14.177, 11.719, 9.282, 6.62, 4.0925, 2.601, 1.149, .667, .3627,
1413 .192, .1497, .10857};
1414
1415 //Absorption index for methane
1416 Float_t abscoMethane[26];
1417 for (i=0;i<26;i++)
1418 {
1419 abscoMethane[i]=AbsoCH4(ppckov[i]*1e9);
1420 //printf("abscoMethane: %e for energy: %e\n", abscoMethane[i],ppckov[i]*1e9);
1421 }
1422
1423 //Absorption index for opaque quarz, csi and grid, efficiency for all and grid
1424 Float_t abscoOpaqueQuarz[26];
1425 Float_t abscoCsI[26];
1426 Float_t abscoGrid[26];
1427 Float_t efficAll[26];
1428 Float_t efficGrid[26];
1429 for (i=0;i<26;i++)
1430 {
1431 abscoOpaqueQuarz[i]=1e-5;
1432 abscoCsI[i]=1e-4;
1433 abscoGrid[i]=1e-4;
1434 efficAll[i]=1;
1435 efficGrid[i]=1;
1436 //printf ("All must be 1: %e, %e, %e, %e, %e\n",abscoOpaqueQuarz[i],abscoCsI[i],abscoGrid[i],efficAll[i],efficGrid[i]);
1437 }
1438
1439 //Efficiency for csi
1440
1441 Float_t efficCsI[26] = {0.000199999995, 0.000600000028, 0.000699999975, 0.00499999989, 0.00749999983, 0.010125,
1442 0.0242999997, 0.0405000001, 0.0688500032, 0.105299994, 0.121500008, 0.141749993, 0.157949999,
1443 0.162, 0.166050002, 0.167669997, 0.174299985, 0.176789999, 0.179279998, 0.182599992, 0.18592,
1444 0.187579989, 0.189239994, 0.190899998, 0.207499996, 0.215799987};
1445
1446
1447
1448 //FRESNEL LOSS CORRECTION FOR PERPENDICULAR INCIDENCE AND
1449 //UNPOLARIZED PHOTONS
1450
1451 for (i=0;i<26;i++)
1452 {
1453 efficCsI[i] = efficCsI[i]/(1.-Fresnel(ppckov[i]*1e9,1.,0));
1454 //printf ("Fresnel result: %e for energy: %e\n",Fresnel(ppckov[i]*1e9,1.,0),ppckov[i]*1e9);
1455 }
1456
1457 /*******************************************End of rich_media.f***************************************/
1458
1459
1460
1461
1462
1463
1464 Float_t afre[2], agri, amet[2], aqua[2], ahon, zfre[2], zgri, zhon,
1465 zmet[2], zqua[2];
1466 Int_t nlmatfre;
1467 Float_t densquao;
1468 Int_t nlmatmet, nlmatqua;
1469 Float_t wmatquao[2], rIndexFreon[26];
1470 Float_t aquao[2], epsil, stmin, zquao[2];
1471 Int_t nlmatquao;
1472 Float_t radlal, densal, tmaxfd, deemax, stemax;
1473 Float_t aal, zal, radlgri, densfre, radlhon, densgri, denshon,densqua, densmet, wmatfre[2], wmatmet[2], wmatqua[2];
1474
1475 Int_t *idtmed = fIdtmed->GetArray()-999;
1476
452a64c6 1477 // --- Photon energy (GeV)
1478 // --- Refraction indexes
1479 for (i = 0; i < 26; ++i) {
1480 rIndexFreon[i] = ppckov[i] * .0172 * 1e9 + 1.177;
1481 //rIndexFreon[i] = 1;
1482 //printf ("rIndexFreon: %e \n efficCsI: %e for energy: %e\n",rIndexFreon[i], efficCsI[i], ppckov[i]);
1483 }
1484
1485 // --- Detection efficiencies (quantum efficiency for CsI)
1486 // --- Define parameters for honeycomb.
1487 // Used carbon of equivalent rad. lenght
1488
1489 ahon = 12.01;
1490 zhon = 6.;
9dda3582 1491 denshon = 0.1;
452a64c6 1492 radlhon = 18.8;
1493
1494 // --- Parameters to include in GSMIXT, relative to Quarz (SiO2)
1495
1496 aqua[0] = 28.09;
1497 aqua[1] = 16.;
1498 zqua[0] = 14.;
1499 zqua[1] = 8.;
1500 densqua = 2.64;
1501 nlmatqua = -2;
1502 wmatqua[0] = 1.;
1503 wmatqua[1] = 2.;
1504
1505 // --- Parameters to include in GSMIXT, relative to opaque Quarz (SiO2)
1506
1507 aquao[0] = 28.09;
1508 aquao[1] = 16.;
1509 zquao[0] = 14.;
1510 zquao[1] = 8.;
1511 densquao = 2.64;
1512 nlmatquao = -2;
1513 wmatquao[0] = 1.;
1514 wmatquao[1] = 2.;
1515
1516 // --- Parameters to include in GSMIXT, relative to Freon (C6F14)
1517
1518 afre[0] = 12.;
1519 afre[1] = 19.;
1520 zfre[0] = 6.;
1521 zfre[1] = 9.;
1522 densfre = 1.7;
1523 nlmatfre = -2;
1524 wmatfre[0] = 6.;
1525 wmatfre[1] = 14.;
1526
1527 // --- Parameters to include in GSMIXT, relative to methane (CH4)
1528
1529 amet[0] = 12.01;
1530 amet[1] = 1.;
1531 zmet[0] = 6.;
1532 zmet[1] = 1.;
1533 densmet = 7.17e-4;
1534 nlmatmet = -2;
1535 wmatmet[0] = 1.;
1536 wmatmet[1] = 4.;
1537
1538 // --- Parameters to include in GSMIXT, relative to anode grid (Cu)
1539
1540 agri = 63.54;
1541 zgri = 29.;
1542 densgri = 8.96;
1543 radlgri = 1.43;
1544
1545 // --- Parameters to include in GSMATE related to aluminium sheet
1546
1547 aal = 26.98;
1548 zal = 13.;
1549 densal = 2.7;
1550 radlal = 8.9;
6197ac87 1551
1552 // --- Glass parameters
1553
1554 Float_t aglass[5]={12.01, 28.09, 16., 10.8, 23.};
1555 Float_t zglass[5]={ 6., 14., 8., 5., 11.};
1556 Float_t wglass[5]={ 0.5, 0.105, 0.355, 0.03, 0.01};
1557 Float_t dglass=1.74;
1558
452a64c6 1559
1560 AliMaterial(1, "Air $", 14.61, 7.3, .001205, 30420., 67500);
1561 AliMaterial(6, "HON", ahon, zhon, denshon, radlhon, 0);
1562 AliMaterial(16, "CSI", ahon, zhon, denshon, radlhon, 0);
1563 AliMixture(20, "QUA", aqua, zqua, densqua, nlmatqua, wmatqua);
1564 AliMixture(21, "QUAO", aquao, zquao, densquao, nlmatquao, wmatquao);
1565 AliMixture(30, "FRE", afre, zfre, densfre, nlmatfre, wmatfre);
1566 AliMixture(40, "MET", amet, zmet, densmet, nlmatmet, wmatmet);
1567 AliMixture(41, "METG", amet, zmet, densmet, nlmatmet, wmatmet);
1568 AliMaterial(11, "GRI", agri, zgri, densgri, radlgri, 0);
1569 AliMaterial(50, "ALUM", aal, zal, densal, radlal, 0);
6197ac87 1570 AliMixture(32, "GLASS",aglass, zglass, dglass, 5, wglass);
1571 AliMaterial(31, "COPPER$", 63.54, 29., 8.96, 1.4, 0.);
452a64c6 1572
1573 tmaxfd = -10.;
1574 stemax = -.1;
1575 deemax = -.2;
1576 epsil = .001;
1577 stmin = -.001;
1578
1579 AliMedium(1, "DEFAULT MEDIUM AIR$", 1, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1580 AliMedium(2, "HONEYCOMB$", 6, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1581 AliMedium(3, "QUARZO$", 20, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1582 AliMedium(4, "FREON$", 30, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1583 AliMedium(5, "METANO$", 40, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1584 AliMedium(6, "CSI$", 16, 1, isxfld, sxmgmx,tmaxfd, stemax, deemax, epsil, stmin);
1585 AliMedium(7, "GRIGLIA$", 11, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1586 AliMedium(8, "QUARZOO$", 21, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1587 AliMedium(9, "GAP$", 41, 1, isxfld, sxmgmx,tmaxfd, .1, -deemax, epsil, -stmin);
1588 AliMedium(10, "ALUMINUM$", 50, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
6197ac87 1589 AliMedium(11, "GLASS", 32, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1590 AliMedium(12, "PCB_COPPER", 31, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
452a64c6 1591
1592
5cfcdf54 1593 gMC->SetCerenkov(idtmed[1000], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1594 gMC->SetCerenkov(idtmed[1001], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1595 gMC->SetCerenkov(idtmed[1002], 26, ppckov, abscoQuarz, efficAll,rIndexQuarz);
1596 gMC->SetCerenkov(idtmed[1003], 26, ppckov, abscoFreon, efficAll,rIndexFreon);
1597 gMC->SetCerenkov(idtmed[1004], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1598 gMC->SetCerenkov(idtmed[1005], 26, ppckov, abscoCsI, efficCsI, rIndexMethane);
1599 gMC->SetCerenkov(idtmed[1006], 26, ppckov, abscoGrid, efficGrid, rIndexGrid);
1600 gMC->SetCerenkov(idtmed[1007], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz);
1601 gMC->SetCerenkov(idtmed[1008], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1602 gMC->SetCerenkov(idtmed[1009], 26, ppckov, abscoGrid, efficGrid, rIndexGrid);
1603 gMC->SetCerenkov(idtmed[1010], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz);
452a64c6 1604}
1605
1606//___________________________________________
1607
1608Float_t AliRICH::Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
1609{
1610
1611 //ENE(EV), PDOTI=COS(INC.ANG.), PDOTR=COS(POL.PLANE ROT.ANG.)
1612
1613 Float_t en[36] = {5.0,5.1,5.2,5.3,5.4,5.5,5.6,5.7,5.8,5.9,6.0,6.1,6.2,
1614 6.3,6.4,6.5,6.6,6.7,6.8,6.9,7.0,7.1,7.2,7.3,7.4,7.5,7.6,7.7,
1615 7.8,7.9,8.0,8.1,8.2,8.3,8.4,8.5};
1616
1617
1618 Float_t csin[36] = {2.14,2.21,2.33,2.48,2.76,2.97,2.99,2.59,2.81,3.05,
1619 2.86,2.53,2.55,2.66,2.79,2.96,3.18,3.05,2.84,2.81,2.38,2.11,
1620 2.01,2.13,2.39,2.73,3.08,3.15,2.95,2.73,2.56,2.41,2.12,1.95,
1621 1.72,1.53};
1622
1623 Float_t csik[36] = {0.,0.,0.,0.,0.,0.196,0.408,0.208,0.118,0.49,0.784,0.543,
1624 0.424,0.404,0.371,0.514,0.922,1.102,1.139,1.376,1.461,1.253,0.878,
1625 0.69,0.612,0.649,0.824,1.347,1.571,1.678,1.763,1.857,1.824,1.824,
1626 1.714,1.498};
1627 Float_t xe=ene;
1628 Int_t j=Int_t(xe*10)-49;
1629 Float_t cn=csin[j]+((csin[j+1]-csin[j])/0.1)*(xe-en[j]);
1630 Float_t ck=csik[j]+((csik[j+1]-csik[j])/0.1)*(xe-en[j]);
1631
1632 //FORMULAE FROM HANDBOOK OF OPTICS, 33.23 OR
1633 //W.R. HUNTER, J.O.S.A. 54 (1964),15 , J.O.S.A. 55(1965),1197
1634
1635 Float_t sinin=TMath::Sqrt(1-pdoti*pdoti);
1636 Float_t tanin=sinin/pdoti;
1637
1638 Float_t c1=cn*cn-ck*ck-sinin*sinin;
1639 Float_t c2=4*cn*cn*ck*ck;
1640 Float_t aO=TMath::Sqrt(0.5*(TMath::Sqrt(c1*c1+c2)+c1));
1641 Float_t b2=0.5*(TMath::Sqrt(c1*c1+c2)-c1);
1642
1643 Float_t rs=((aO-pdoti)*(aO-pdoti)+b2)/((aO+pdoti)*(aO+pdoti)+b2);
1644 Float_t rp=rs*((aO-sinin*tanin)*(aO-sinin*tanin)+b2)/((aO+sinin*tanin)*(aO+sinin*tanin)+b2);
1645
1646
1647 //CORRECTION FACTOR FOR SURFACE ROUGHNESS
1648 //B.J. STAGG APPLIED OPTICS, 30(1991),4113
1649
1650 Float_t sigraf=18.;
1651 Float_t lamb=1240/ene;
1652 Float_t fresn;
1653
1654 Float_t rO=TMath::Exp(-(4*TMath::Pi()*pdoti*sigraf/lamb)*(4*TMath::Pi()*pdoti*sigraf/lamb));
1655
1656 if(pola)
1657 {
1658 Float_t pdotr=0.8; //DEGREE OF POLARIZATION : 1->P , -1->S
1659 fresn=0.5*(rp*(1+pdotr)+rs*(1-pdotr));
1660 }
1661 else
1662 fresn=0.5*(rp+rs);
1663
1664 fresn = fresn*rO;
1665 return(fresn);
1666}
1667
1668//__________________________________________
1669Float_t AliRICH::AbsoCH4(Float_t x)
1670{
1671
1672 //KLOSCH,SCH4(9),WL(9),EM(9),ALENGTH(31)
1673 Float_t sch4[9] = {.12,.16,.23,.38,.86,2.8,7.9,28.,80.}; //MB X 10^22
1674 //Float_t wl[9] = {153.,152.,151.,150.,149.,148.,147.,146.,145};
1675 Float_t em[9] = {8.1,8.158,8.212,8.267,8.322,8.378,8.435,8.493,8.55};
1676 const Float_t kLosch=2.686763E19; // LOSCHMIDT NUMBER IN CM-3
1677 const Float_t kIgas1=100, kIgas2=0, kOxy=10., kWater=5., kPressure=750.,kTemperature=283.;
1678 Float_t pn=kPressure/760.;
1679 Float_t tn=kTemperature/273.16;
1680
1681
1682// ------- METHANE CROSS SECTION -----------------
1683// ASTROPH. J. 214, L47 (1978)
1684
1685 Float_t sm=0;
1686 if (x<7.75)
1687 sm=.06e-22;
1688
1689 if(x>=7.75 && x<=8.1)
1690 {
1691 Float_t c0=-1.655279e-1;
1692 Float_t c1=6.307392e-2;
1693 Float_t c2=-8.011441e-3;
1694 Float_t c3=3.392126e-4;
1695 sm=(c0+c1*x+c2*x*x+c3*x*x*x)*1.e-18;
1696 }
1697
1698 if (x> 8.1)
1699 {
1700 Int_t j=0;
1701 while (x<=em[j] && x>=em[j+1])
1702 {
1703 j++;
1704 Float_t a=(sch4[j+1]-sch4[j])/(em[j+1]-em[j]);
1705 sm=(sch4[j]+a*(x-em[j]))*1e-22;
1706 }
1707 }
1708
1709 Float_t dm=(kIgas1/100.)*(1.-((kOxy+kWater)/1.e6))*kLosch*pn/tn;
1710 Float_t abslm=1./sm/dm;
1711
1712// ------- ISOBUTHANE CROSS SECTION --------------
1713// i-C4H10 (ai) abs. length from curves in
1714// Lu-McDonald paper for BARI RICH workshop .
1715// -----------------------------------------------------------
1716
1717 Float_t ai;
1718 Float_t absli;
1719 if (kIgas2 != 0)
1720 {
1721 if (x<7.25)
1722 ai=100000000.;
1723
1724 if(x>=7.25 && x<7.375)
1725 ai=24.3;
1726
1727 if(x>=7.375)
1728 ai=.0000000001;
1729
1730 Float_t si = 1./(ai*kLosch*273.16/293.); // ISOB. CRO.SEC.IN CM2
1731 Float_t di=(kIgas2/100.)*(1.-((kOxy+kWater)/1.e6))*kLosch*pn/tn;
1732 absli =1./si/di;
1733 }
1734 else
1735 absli=1.e18;
1736// ---------------------------------------------------------
1737//
1738// transmission of O2
1739//
1740// y= path in cm, x=energy in eV
1741// so= cross section for UV absorption in cm2
1742// do= O2 molecular density in cm-3
1743// ---------------------------------------------------------
1744
1745 Float_t abslo;
1746 Float_t so=0;
1747 if(x>=6.0)
1748 {
1749 if(x>=6.0 && x<6.5)
1750 {
1751 so=3.392709e-13 * TMath::Exp(2.864104 *x);
1752 so=so*1e-18;
1753 }
1754
1755 if(x>=6.5 && x<7.0)
1756 {
1757 so=2.910039e-34 * TMath::Exp(10.3337*x);
1758 so=so*1e-18;
1759 }
1760
1761
1762 if (x>=7.0)
1763 {
1764 Float_t a0=-73770.76;
1765 Float_t a1=46190.69;
1766 Float_t a2=-11475.44;
1767 Float_t a3=1412.611;
1768 Float_t a4=-86.07027;
1769 Float_t a5=2.074234;
1770 so= a0+(a1*x)+(a2*x*x)+(a3*x*x*x)+(a4*x*x*x*x)+(a5*x*x*x*x*x);
1771 so=so*1e-18;
1772 }
1773
1774 Float_t dox=(kOxy/1e6)*kLosch*pn/tn;
1775 abslo=1./so/dox;
1776 }
1777 else
1778 abslo=1.e18;
1779// ---------------------------------------------------------
1780//
1781// transmission of H2O
1782//
1783// y= path in cm, x=energy in eV
1784// sw= cross section for UV absorption in cm2
1785// dw= H2O molecular density in cm-3
1786// ---------------------------------------------------------
1787
1788 Float_t abslw;
1789
1790 Float_t b0=29231.65;
1791 Float_t b1=-15807.74;
1792 Float_t b2=3192.926;
1793 Float_t b3=-285.4809;
1794 Float_t b4=9.533944;
1795
1796 if(x>6.75)
1797 {
1798 Float_t sw= b0+(b1*x)+(b2*x*x)+(b3*x*x*x)+(b4*x*x*x*x);
1799 sw=sw*1e-18;
1800 Float_t dw=(kWater/1e6)*kLosch*pn/tn;
1801 abslw=1./sw/dw;
1802 }
1803 else
1804 abslw=1.e18;
1805
1806// ---------------------------------------------------------
1807
1808 Float_t alength=1./(1./abslm+1./absli+1./abslo+1./abslw);
1809 return (alength);
1810}
1811
1812
1813
ddae0931 1814//___________________________________________
1815Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
1816{
237c933d 1817
1818// Default value
1819
ddae0931 1820 return 9999;
1821}
fe4da5cc 1822
ddae0931 1823//___________________________________________
2ab0c725 1824void AliRICH::MakeBranch(Option_t* option, char *file)
fe4da5cc 1825{
237c933d 1826 // Create Tree branches for the RICH.
fe4da5cc 1827
237c933d 1828 const Int_t kBufferSize = 4000;
ddae0931 1829 char branchname[20];
2ab0c725 1830
1831 AliDetector::MakeBranch(option,file);
1832
5cf7bbad 1833 const char *cH = strstr(option,"H");
1834 const char *cD = strstr(option,"D");
1835 const char *cR = strstr(option,"R");
34ead2dd 1836 const char *cS = strstr(option,"S");
1837
2ab0c725 1838
1839 if (cH) {
1840 sprintf(branchname,"%sCerenkov",GetName());
1841 if (fCerenkovs && gAlice->TreeH()) {
c3ecfac9 1842 //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeH(),branchname, &fCerenkovs, kBufferSize, file) ;
1843 gAlice->MakeBranchInTree(gAlice->TreeH(),branchname, &fCerenkovs, kBufferSize, file) ;
1844 //branch->SetAutoDelete(kFALSE);
b251a2b5 1845 }
1846 sprintf(branchname,"%sSDigits",GetName());
1847 if (fSDigits && gAlice->TreeH()) {
c3ecfac9 1848 //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeH(),branchname, &fSDigits, kBufferSize, file) ;
1849 gAlice->MakeBranchInTree(gAlice->TreeH(),branchname, &fSDigits, kBufferSize, file) ;
1850 //branch->SetAutoDelete(kFALSE);
34ead2dd 1851 //printf("Making branch %sSDigits in TreeH\n",GetName());
2ab0c725 1852 }
b251a2b5 1853 }
1854
34ead2dd 1855 if (cS) {
b251a2b5 1856 sprintf(branchname,"%sSDigits",GetName());
1857 if (fSDigits && gAlice->TreeS()) {
c3ecfac9 1858 //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeS(),branchname, &fSDigits, kBufferSize, file) ;
1859 gAlice->MakeBranchInTree(gAlice->TreeS(),branchname, &fSDigits, kBufferSize, file) ;
1860 //branch->SetAutoDelete(kFALSE);
34ead2dd 1861 //printf("Making branch %sSDigits in TreeS\n",GetName());
b251a2b5 1862 }
34ead2dd 1863 }
fe4da5cc 1864
2ab0c725 1865 if (cD) {
1866 //
1867 // one branch for digits per chamber
1868 //
1869 Int_t i;
1870
1871 for (i=0; i<kNCH ;i++) {
633e9a38 1872 sprintf(branchname,"%sDigits%d",GetName(),i+1);
1873 if (fDchambers && gAlice->TreeD()) {
c3ecfac9 1874 //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeD(),branchname, &((*fDchambers)[i]), kBufferSize, file) ;
1875 gAlice->MakeBranchInTree(gAlice->TreeD(),branchname, &((*fDchambers)[i]), kBufferSize, file) ;
1876 //branch->SetAutoDelete(kFALSE);
633e9a38 1877 //printf("Making Branch %sDigits%d\n",GetName(),i+1);
1878 }
2ab0c725 1879 }
fe4da5cc 1880 }
2e5f0f7b 1881
2ab0c725 1882 if (cR) {
1883 //
1884 // one branch for raw clusters per chamber
1885 //
34ead2dd 1886
1887 //printf("Called MakeBranch for TreeR\n");
1888
2ab0c725 1889 Int_t i;
1890
1891 for (i=0; i<kNCH ;i++) {
1892 sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
1893 if (fRawClusters && gAlice->TreeR()) {
c3ecfac9 1894 //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRawClusters)[i]), kBufferSize, file) ;
1895 gAlice->MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRawClusters)[i]), kBufferSize, file) ;
1896 //branch->SetAutoDelete(kFALSE);
2ab0c725 1897 }
1898 }
1899 //
1900 // one branch for rec hits per chamber
1901 //
1902 for (i=0; i<kNCH ;i++) {
1903 sprintf(branchname,"%sRecHits1D%d",GetName(),i+1);
1904 if (fRecHits1D && gAlice->TreeR()) {
c3ecfac9 1905 //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits1D)[i]), kBufferSize, file) ;
1906 gAlice->MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits1D)[i]), kBufferSize, file) ;
1907 //branch->SetAutoDelete(kFALSE);
2ab0c725 1908 }
1909 }
1910 for (i=0; i<kNCH ;i++) {
1911 sprintf(branchname,"%sRecHits3D%d",GetName(),i+1);
1912 if (fRecHits3D && gAlice->TreeR()) {
c3ecfac9 1913 //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits3D)[i]), kBufferSize, file) ;
a3d71079 1914 gAlice->MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits3D)[i]), kBufferSize, file) ;
c3ecfac9 1915 //branch->SetAutoDelete(kFALSE);
2e5f0f7b 1916 }
2ab0c725 1917 }
1918 }
fe4da5cc 1919}
1920
ddae0931 1921//___________________________________________
1922void AliRICH::SetTreeAddress()
fe4da5cc 1923{
237c933d 1924 // Set branch address for the Hits and Digits Tree.
2e5f0f7b 1925 char branchname[20];
1926 Int_t i;
1927
ddae0931 1928 AliDetector::SetTreeAddress();
1929
1930 TBranch *branch;
1931 TTree *treeH = gAlice->TreeH();
1932 TTree *treeD = gAlice->TreeD();
2e5f0f7b 1933 TTree *treeR = gAlice->TreeR();
34ead2dd 1934 TTree *treeS = gAlice->TreeS();
ddae0931 1935
1936 if (treeH) {
b251a2b5 1937 if (fCerenkovs) {
ddae0931 1938 branch = treeH->GetBranch("RICHCerenkov");
1939 if (branch) branch->SetAddress(&fCerenkovs);
fe4da5cc 1940 }
633e9a38 1941 if (fSDigits) {
b251a2b5 1942 branch = treeH->GetBranch("RICHSDigits");
34ead2dd 1943 if (branch)
1944 {
1945 branch->SetAddress(&fSDigits);
1946 //printf("Setting sdigits branch address at %p in TreeH\n",&fSDigits);
1947 }
b251a2b5 1948 }
ddae0931 1949 }
1950
34ead2dd 1951 if (treeS) {
b251a2b5 1952 if (fSDigits) {
1953 branch = treeS->GetBranch("RICHSDigits");
34ead2dd 1954 if (branch)
1955 {
1956 branch->SetAddress(&fSDigits);
1957 //printf("Setting sdigits branch address at %p in TreeS\n",&fSDigits);
1958 }
b251a2b5 1959 }
34ead2dd 1960 }
b251a2b5 1961
1962
ddae0931 1963 if (treeD) {
237c933d 1964 for (int i=0; i<kNCH; i++) {
ddae0931 1965 sprintf(branchname,"%sDigits%d",GetName(),i+1);
1966 if (fDchambers) {
34ead2dd 1967 branch = treeD->GetBranch(branchname);
1968 if (branch) branch->SetAddress(&((*fDchambers)[i]));
ddae0931 1969 }
fe4da5cc 1970 }
fe4da5cc 1971 }
2e5f0f7b 1972 if (treeR) {
237c933d 1973 for (i=0; i<kNCH; i++) {
2e5f0f7b 1974 sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
1975 if (fRawClusters) {
1976 branch = treeR->GetBranch(branchname);
1977 if (branch) branch->SetAddress(&((*fRawClusters)[i]));
1978 }
1979 }
1980
237c933d 1981 for (i=0; i<kNCH; i++) {
a4622d0f 1982 sprintf(branchname,"%sRecHits1D%d",GetName(),i+1);
1983 if (fRecHits1D) {
2e5f0f7b 1984 branch = treeR->GetBranch(branchname);
a4622d0f 1985 if (branch) branch->SetAddress(&((*fRecHits1D)[i]));
2e5f0f7b 1986 }
1987 }
1988
a4622d0f 1989 for (i=0; i<kNCH; i++) {
1990 sprintf(branchname,"%sRecHits3D%d",GetName(),i+1);
1991 if (fRecHits3D) {
1992 branch = treeR->GetBranch(branchname);
1993 if (branch) branch->SetAddress(&((*fRecHits3D)[i]));
1994 }
1995 }
1996
2e5f0f7b 1997 }
ddae0931 1998}
1999//___________________________________________
2000void AliRICH::ResetHits()
2001{
237c933d 2002 // Reset number of clusters and the cluster array for this detector
ddae0931 2003 AliDetector::ResetHits();
b251a2b5 2004 fNSDigits = 0;
2e5f0f7b 2005 fNcerenkovs = 0;
b251a2b5 2006 if (fSDigits) fSDigits->Clear();
ddae0931 2007 if (fCerenkovs) fCerenkovs->Clear();
fe4da5cc 2008}
2009
2e5f0f7b 2010
ddae0931 2011//____________________________________________
2012void AliRICH::ResetDigits()
2013{
237c933d 2014 //
2015 // Reset number of digits and the digits array for this detector
2016 //
2017 for ( int i=0;i<kNCH;i++ ) {
2ab0c725 2018 if (fDchambers && (*fDchambers)[i]) (*fDchambers)[i]->Clear();
ddae0931 2019 if (fNdch) fNdch[i]=0;
2020 }
2021}
2e5f0f7b 2022
ddae0931 2023//____________________________________________
2e5f0f7b 2024void AliRICH::ResetRawClusters()
fe4da5cc 2025{
237c933d 2026 //
2027 // Reset number of raw clusters and the raw clust array for this detector
2028 //
2029 for ( int i=0;i<kNCH;i++ ) {
2e5f0f7b 2030 if ((*fRawClusters)[i]) ((TClonesArray*)(*fRawClusters)[i])->Clear();
2031 if (fNrawch) fNrawch[i]=0;
fe4da5cc 2032 }
ddae0931 2033}
fe4da5cc 2034
2e5f0f7b 2035//____________________________________________
a4622d0f 2036void AliRICH::ResetRecHits1D()
fe4da5cc 2037{
237c933d 2038 //
2039 // Reset number of raw clusters and the raw clust array for this detector
2040 //
2e5f0f7b 2041
237c933d 2042 for ( int i=0;i<kNCH;i++ ) {
a4622d0f 2043 if ((*fRecHits1D)[i]) ((TClonesArray*)(*fRecHits1D)[i])->Clear();
2044 if (fNrechits1D) fNrechits1D[i]=0;
2045 }
2046}
2047
2048//____________________________________________
2049void AliRICH::ResetRecHits3D()
2050{
2051 //
2052 // Reset number of raw clusters and the raw clust array for this detector
2053 //
2054
2055 for ( int i=0;i<kNCH;i++ ) {
2056 if ((*fRecHits3D)[i]) ((TClonesArray*)(*fRecHits3D)[i])->Clear();
2057 if (fNrechits3D) fNrechits3D[i]=0;
2e5f0f7b 2058 }
fe4da5cc 2059}
2060
ddae0931 2061//___________________________________________
2e5f0f7b 2062void AliRICH::SetGeometryModel(Int_t id, AliRICHGeometry *geometry)
fe4da5cc 2063{
237c933d 2064
2065//
2066// Setter for the RICH geometry model
2067//
2068
2069
2e5f0f7b 2070 ((AliRICHChamber*) (*fChambers)[id])->GeometryModel(geometry);
fe4da5cc 2071}
2072
ddae0931 2073//___________________________________________
a2f7eaf6 2074void AliRICH::SetSegmentationModel(Int_t id, AliSegmentation *segmentation)
ddae0931 2075{
237c933d 2076
2077//
2078// Setter for the RICH segmentation model
2079//
2080
a2f7eaf6 2081 ((AliRICHChamber*) (*fChambers)[id])->SetSegmentationModel(segmentation);
ddae0931 2082}
fe4da5cc 2083
ddae0931 2084//___________________________________________
2e5f0f7b 2085void AliRICH::SetResponseModel(Int_t id, AliRICHResponse *response)
fe4da5cc 2086{
237c933d 2087
2088//
2089// Setter for the RICH response model
2090//
2091
2e5f0f7b 2092 ((AliRICHChamber*) (*fChambers)[id])->ResponseModel(response);
fe4da5cc 2093}
2094
2e5f0f7b 2095void AliRICH::SetReconstructionModel(Int_t id, AliRICHClusterFinder *reconst)
ddae0931 2096{
237c933d 2097
2098//
2099// Setter for the RICH reconstruction model (clusters)
2100//
2101
a2f7eaf6 2102 ((AliRICHChamber*) (*fChambers)[id])->SetReconstructionModel(reconst);
fe4da5cc 2103}
2104
ddae0931 2105//___________________________________________
ddae0931 2106void AliRICH::StepManager()
fe4da5cc 2107{
237c933d 2108
452a64c6 2109// Full Step Manager
2110
2111 Int_t copy, id;
2112 static Int_t idvol;
2113 static Int_t vol[2];
2114 Int_t ipart;
c24372d0 2115 static Float_t hits[22];
452a64c6 2116 static Float_t ckovData[19];
2117 TLorentzVector position;
2118 TLorentzVector momentum;
2119 Float_t pos[3];
2120 Float_t mom[4];
2121 Float_t localPos[3];
2122 Float_t localMom[4];
2123 Float_t localTheta,localPhi;
2124 Float_t theta,phi;
2125 Float_t destep, step;
2126 Float_t ranf[2];
2127 Int_t nPads;
2128 Float_t coscerenkov;
2129 static Float_t eloss, xhit, yhit, tlength;
2130 const Float_t kBig=1.e10;
2131
2132 TClonesArray &lhits = *fHits;
452a64c6 2133 TParticle *current = (TParticle*)(*gAlice->Particles())[gAlice->CurrentTrack()];
237c933d 2134
452a64c6 2135 //if (current->Energy()>1)
2136 //{
2137
2138 // Only gas gap inside chamber
2139 // Tag chambers and record hits when track enters
2140
2141 idvol=-1;
2142 id=gMC->CurrentVolID(copy);
2143 Float_t cherenkovLoss=0;
2144 //gAlice->KeepTrack(gAlice->CurrentTrack());
2145
2146 gMC->TrackPosition(position);
2147 pos[0]=position(0);
2148 pos[1]=position(1);
2149 pos[2]=position(2);
8cda28a5 2150 //bzero((char *)ckovData,sizeof(ckovData)*19);
452a64c6 2151 ckovData[1] = pos[0]; // X-position for hit
2152 ckovData[2] = pos[1]; // Y-position for hit
2153 ckovData[3] = pos[2]; // Z-position for hit
3cc8edee 2154 ckovData[6] = 0; // dummy track length
452a64c6 2155 //ckovData[11] = gAlice->CurrentTrack();
8cda28a5 2156
2157 //printf("\n+++++++++++\nTrack: %d\n++++++++++++\n",gAlice->CurrentTrack());
452a64c6 2158
2159 //AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH");
2160
2161 /********************Store production parameters for Cerenkov photons************************/
2162//is it a Cerenkov photon?
8cda28a5 2163 if (gMC->TrackPid() == 50000050) {
452a64c6 2164
2165 //if (gMC->VolId("GAP ")==gMC->CurrentVolID(copy))
2166 //{
2167 Float_t ckovEnergy = current->Energy();
2168 //energy interval for tracking
2169 if (ckovEnergy > 5.6e-09 && ckovEnergy < 7.8e-09 )
2170 //if (ckovEnergy > 0)
2171 {
8cda28a5 2172 if (gMC->IsTrackEntering()){ //is track entering?
2173 //printf("Track entered (1)\n");
452a64c6 2174 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2175 { //is it in freo?
5cfcdf54 2176 if (gMC->IsNewTrack()){ //is it the first step?
8cda28a5 2177 //printf("I'm in!\n");
452a64c6 2178 Int_t mother = current->GetFirstMother();
2179
2180 //printf("Second Mother:%d\n",current->GetSecondMother());
2181
2182 ckovData[10] = mother;
2183 ckovData[11] = gAlice->CurrentTrack();
2184 ckovData[12] = 1; //Media where photon was produced 1->Freon, 2->Quarz
8cda28a5 2185 //printf("Produced in FREO\n");
452a64c6 2186 fCkovNumber++;
2187 fFreonProd=1;
2188 //printf("Index: %d\n",fCkovNumber);
2189 } //first step question
2190 } //freo question
2191
5cfcdf54 2192 if (gMC->IsNewTrack()){ //is it first step?
452a64c6 2193 if (gMC->VolId("QUAR")==gMC->CurrentVolID(copy)) //is it in quarz?
2194 {
2195 ckovData[12] = 2;
8cda28a5 2196 //printf("Produced in QUAR\n");
452a64c6 2197 } //quarz question
2198 } //first step question
2199
2200 //printf("Before %d\n",fFreonProd);
2201 } //track entering question
2202
2203 if (ckovData[12] == 1) //was it produced in Freon?
2204 //if (fFreonProd == 1)
2205 {
2206 if (gMC->IsTrackEntering()){ //is track entering?
8cda28a5 2207 //printf("Track entered (2)\n");
2208 //printf("Current volume (should be META): %s\n",gMC->CurrentVolName());
2209 //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("META"),gMC->CurrentVolID(copy));
452a64c6 2210 if (gMC->VolId("META")==gMC->CurrentVolID(copy)) //is it in gap?
2211 {
8cda28a5 2212 //printf("Got in META\n");
452a64c6 2213 gMC->TrackMomentum(momentum);
2214 mom[0]=momentum(0);
2215 mom[1]=momentum(1);
2216 mom[2]=momentum(2);
2217 mom[3]=momentum(3);
2218 // Z-position for hit
2219
2220
2221 /**************** Photons lost in second grid have to be calculated by hand************/
2222
2223 Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
2224 Float_t t = (1. - .025 / cophi) * (1. - .05 / cophi);
2225 gMC->Rndm(ranf, 1);
2226 //printf("grid calculation:%f\n",t);
2227 if (ranf[0] > t) {
5cfcdf54 2228 gMC->StopTrack();
452a64c6 2229 ckovData[13] = 5;
2230 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
8cda28a5 2231 //printf("Added One (1)!\n");
452a64c6 2232 //printf("Lost one in grid\n");
2233 }
2234 /**********************************************************************************/
2235 } //gap
2236
8cda28a5 2237 //printf("Current volume (should be CSI) (1): %s\n",gMC->CurrentVolName());
2238 //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
452a64c6 2239 if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy)) //is it in csi?
2240 {
8cda28a5 2241 //printf("Got in CSI\n");
452a64c6 2242 gMC->TrackMomentum(momentum);
2243 mom[0]=momentum(0);
2244 mom[1]=momentum(1);
2245 mom[2]=momentum(2);
2246 mom[3]=momentum(3);
2247
2248 /********* Photons lost by Fresnel reflection have to be calculated by hand********/
2249 /***********************Cerenkov phtons (always polarised)*************************/
2250
2251 Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
2252 Float_t t = Fresnel(ckovEnergy*1e9,cophi,1);
2253 gMC->Rndm(ranf, 1);
2254 if (ranf[0] < t) {
5cfcdf54 2255 gMC->StopTrack();
452a64c6 2256 ckovData[13] = 6;
2257 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
8cda28a5 2258 //printf("Added One (2)!\n");
452a64c6 2259 //printf("Lost by Fresnel\n");
2260 }
2261 /**********************************************************************************/
2262 }
2263 } //track entering?
2264
2265
2266 /********************Evaluation of losses************************/
2267 /******************still in the old fashion**********************/
2268
5cfcdf54 2269 TArrayI procs;
2270 Int_t i1 = gMC->StepProcesses(procs); //number of physics mechanisms acting on the particle
452a64c6 2271 for (Int_t i = 0; i < i1; ++i) {
2272 // Reflection loss
5cfcdf54 2273 if (procs[i] == kPLightReflection) { //was it reflected
452a64c6 2274 ckovData[13]=10;
2275 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2276 ckovData[13]=1;
2277 if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
2278 ckovData[13]=2;
5cfcdf54 2279 //gMC->StopTrack();
9dda3582 2280 //AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
452a64c6 2281 } //reflection question
8cda28a5 2282
452a64c6 2283 // Absorption loss
5cfcdf54 2284 else if (procs[i] == kPLightAbsorption) { //was it absorbed?
8cda28a5 2285 //printf("Got in absorption\n");
452a64c6 2286 ckovData[13]=20;
2287 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2288 ckovData[13]=11;
2289 if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
2290 ckovData[13]=12;
2291 if (gMC->CurrentVolID(copy) == gMC->VolId("META"))
2292 ckovData[13]=13;
2293 if (gMC->CurrentVolID(copy) == gMC->VolId("GAP "))
2294 ckovData[13]=13;
2295
2296 if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC"))
2297 ckovData[13]=15;
2298
2299 // CsI inefficiency
2300 if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
2301 ckovData[13]=16;
2302 }
5cfcdf54 2303 gMC->StopTrack();
452a64c6 2304 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
8cda28a5 2305 //printf("Added One (3)!\n");
452a64c6 2306 //printf("Added cerenkov %d\n",fCkovNumber);
2307 } //absorption question
2308
2309
2310 // Photon goes out of tracking scope
5cfcdf54 2311 else if (procs[i] == kPStop) { //is it below energy treshold?
452a64c6 2312 ckovData[13]=21;
5cfcdf54 2313 gMC->StopTrack();
452a64c6 2314 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
8cda28a5 2315 //printf("Added One (4)!\n");
452a64c6 2316 } // energy treshold question
2317 } //number of mechanisms cycle
2318 /**********************End of evaluation************************/
2319 } //freon production question
2320 } //energy interval question
2321 //}//inside the proximity gap question
2322 } //cerenkov photon question
2323
2324 /**************************************End of Production Parameters Storing*********************/
2325
2326
2327 /*******************************Treat photons that hit the CsI (Ckovs and Feedbacks)************/
2328
2329 if (gMC->TrackPid() == 50000050 || gMC->TrackPid() == 50000051) {
2330 //printf("Cerenkov\n");
a3d71079 2331
2332 //if (gMC->TrackPid() == 50000051)
2333 //printf("Tracking a feedback\n");
2334
2335 if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
452a64c6 2336 {
8cda28a5 2337 //printf("Current volume (should be CSI) (2): %s\n",gMC->CurrentVolName());
2338 //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
2339 //printf("Got in CSI\n");
a3d71079 2340 //printf("Tracking a %d\n",gMC->TrackPid());
452a64c6 2341 if (gMC->Edep() > 0.){
2342 gMC->TrackPosition(position);
2343 gMC->TrackMomentum(momentum);
2344 pos[0]=position(0);
2345 pos[1]=position(1);
2346 pos[2]=position(2);
2347 mom[0]=momentum(0);
2348 mom[1]=momentum(1);
2349 mom[2]=momentum(2);
2350 mom[3]=momentum(3);
2351 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
2352 Double_t rt = TMath::Sqrt(tc);
2353 theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
2354 phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
2355 gMC->Gmtod(pos,localPos,1);
2356 gMC->Gmtod(mom,localMom,2);
2357
2358 gMC->CurrentVolOffID(2,copy);
2359 vol[0]=copy;
2360 idvol=vol[0]-1;
2361
2362 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
2363 //->Sector(localPos[0], localPos[2]);
2364 //printf("Sector:%d\n",sector);
2365
2366 /*if (gMC->TrackPid() == 50000051){
2367 fFeedbacks++;
2368 printf("Feedbacks:%d\n",fFeedbacks);
2369 }*/
2370
2371 ((AliRICHChamber*) (*fChambers)[idvol])
2372 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2373 if(idvol<kNCH) {
2374 ckovData[0] = gMC->TrackPid(); // particle type
2375 ckovData[1] = pos[0]; // X-position for hit
2376 ckovData[2] = pos[1]; // Y-position for hit
2377 ckovData[3] = pos[2]; // Z-position for hit
2378 ckovData[4] = theta; // theta angle of incidence
2379 ckovData[5] = phi; // phi angle of incidence
b251a2b5 2380 ckovData[8] = (Float_t) fNSDigits; // first sdigit
452a64c6 2381 ckovData[9] = -1; // last pad hit
2382 ckovData[13] = 4; // photon was detected
2383 ckovData[14] = mom[0];
2384 ckovData[15] = mom[1];
2385 ckovData[16] = mom[2];
2386
2387 destep = gMC->Edep();
2388 gMC->SetMaxStep(kBig);
2389 cherenkovLoss += destep;
2390 ckovData[7]=cherenkovLoss;
2391
b251a2b5 2392 nPads = Hits2SDigits(localPos[0],localPos[2],cherenkovLoss,idvol,kCerenkov);
ca96c9ea 2393
b251a2b5 2394 if (fNSDigits > (Int_t)ckovData[8]) {
452a64c6 2395 ckovData[8]= ckovData[8]+1;
b251a2b5 2396 ckovData[9]= (Float_t) fNSDigits;
452a64c6 2397 }
2398
ca96c9ea 2399 //printf("Cerenkov loss: %f\n", cherenkovLoss);
2400
452a64c6 2401 ckovData[17] = nPads;
2402 //printf("nPads:%d",nPads);
2403
2404 //TClonesArray *Hits = RICH->Hits();
2405 AliRICHHit *mipHit = (AliRICHHit*) (fHits->UncheckedAt(0));
2406 if (mipHit)
2407 {
2408 mom[0] = current->Px();
2409 mom[1] = current->Py();
2410 mom[2] = current->Pz();
2411 Float_t mipPx = mipHit->fMomX;
2412 Float_t mipPy = mipHit->fMomY;
2413 Float_t mipPz = mipHit->fMomZ;
2414
2415 Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
2416 Float_t rt = TMath::Sqrt(r);
2417 Float_t mipR = mipPx*mipPx + mipPy*mipPy + mipPz*mipPz;
2418 Float_t mipRt = TMath::Sqrt(mipR);
2419 if ((rt*mipRt) > 0)
2420 {
2421 coscerenkov = (mom[0]*mipPx + mom[1]*mipPy + mom[2]*mipPz)/(rt*mipRt);
2422 }
2423 else
2424 {
2425 coscerenkov = 0;
2426 }
2427 Float_t cherenkov = TMath::ACos(coscerenkov);
2428 ckovData[18]=cherenkov;
2429 }
2430 //if (sector != -1)
2431 //{
2432 AddHit(gAlice->CurrentTrack(),vol,ckovData);
2433 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
8cda28a5 2434 //printf("Added One (5)!\n");
452a64c6 2435 //}
2436 }
2437 }
2438 }
2439 }
2440
2441 /***********************************************End of photon hits*********************************************/
2442
2443
2444 /**********************************************Charged particles treatment*************************************/
2445
2446 else if (gMC->TrackCharge())
2447 //else if (1 == 1)
2448 {
2449//If MIP
2450 /*if (gMC->IsTrackEntering())
2451 {
2452 hits[13]=20;//is track entering?
2453 }*/
2454 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2455 {
c24372d0 2456 gMC->TrackMomentum(momentum);
2457 mom[0]=momentum(0);
2458 mom[1]=momentum(1);
2459 mom[2]=momentum(2);
2460 mom[3]=momentum(3);
2461 hits [19] = mom[0];
2462 hits [20] = mom[1];
2463 hits [21] = mom[2];
452a64c6 2464 fFreonProd=1;
2465 }
2466
2467 if (gMC->VolId("GAP ")== gMC->CurrentVolID(copy)) {
2468// Get current particle id (ipart), track position (pos) and momentum (mom)
2469
2470 gMC->CurrentVolOffID(3,copy);
2471 vol[0]=copy;
2472 idvol=vol[0]-1;
2473
2474 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
2475 //->Sector(localPos[0], localPos[2]);
2476 //printf("Sector:%d\n",sector);
2477
2478 gMC->TrackPosition(position);
2479 gMC->TrackMomentum(momentum);
2480 pos[0]=position(0);
2481 pos[1]=position(1);
2482 pos[2]=position(2);
2483 mom[0]=momentum(0);
2484 mom[1]=momentum(1);
2485 mom[2]=momentum(2);
2486 mom[3]=momentum(3);
2487 gMC->Gmtod(pos,localPos,1);
2488 gMC->Gmtod(mom,localMom,2);
2489
2490 ipart = gMC->TrackPid();
2491 //
2492 // momentum loss and steplength in last step
2493 destep = gMC->Edep();
2494 step = gMC->TrackStep();
2495
2496 //
2497 // record hits when track enters ...
2498 if( gMC->IsTrackEntering()) {
2499// gMC->SetMaxStep(fMaxStepGas);
2500 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
2501 Double_t rt = TMath::Sqrt(tc);
2502 theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
2503 phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
2504
2505
2506 Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
2507 Double_t localRt = TMath::Sqrt(localTc);
2508 localTheta = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])))*kRaddeg;
2509 localPhi = Float_t(TMath::ATan2(Double_t(localMom[2]),Double_t(localMom[0])))*kRaddeg;
2510
2511 hits[0] = Float_t(ipart); // particle type
2512 hits[1] = localPos[0]; // X-position for hit
2513 hits[2] = localPos[1]; // Y-position for hit
2514 hits[3] = localPos[2]; // Z-position for hit
2515 hits[4] = localTheta; // theta angle of incidence
2516 hits[5] = localPhi; // phi angle of incidence
b251a2b5 2517 hits[8] = (Float_t) fNSDigits; // first sdigit
452a64c6 2518 hits[9] = -1; // last pad hit
2519 hits[13] = fFreonProd; // did id hit the freon?
2520 hits[14] = mom[0];
2521 hits[15] = mom[1];
2522 hits[16] = mom[2];
c24372d0 2523 hits[18] = 0; // dummy cerenkov angle
452a64c6 2524
2525 tlength = 0;
2526 eloss = 0;
2527 fFreonProd = 0;
2528
2529 Chamber(idvol).LocaltoGlobal(localPos,hits+1);
2530
2531
2532 //To make chamber coordinates x-y had to pass localPos[0], localPos[2]
2533 xhit = localPos[0];
2534 yhit = localPos[2];
2535 // Only if not trigger chamber
2536 if(idvol<kNCH) {
2537 //
2538 // Initialize hit position (cursor) in the segmentation model
2539 ((AliRICHChamber*) (*fChambers)[idvol])
2540 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2541 }
2542 }
2543
2544 //
2545 // Calculate the charge induced on a pad (disintegration) in case
2546 //
2547 // Mip left chamber ...
2548 if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
2549 gMC->SetMaxStep(kBig);
2550 eloss += destep;
2551 tlength += step;
2552
2553
2554 // Only if not trigger chamber
2555 if(idvol<kNCH) {
2556 if (eloss > 0)
2557 {
2558 if(gMC->TrackPid() == kNeutron)
2559 printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
b251a2b5 2560 nPads = Hits2SDigits(xhit,yhit,eloss,idvol,kMip);
452a64c6 2561 hits[17] = nPads;
2562 //printf("nPads:%d",nPads);
2563 }
2564 }
2565
2566 hits[6]=tlength;
2567 hits[7]=eloss;
b251a2b5 2568 if (fNSDigits > (Int_t)hits[8]) {
452a64c6 2569 hits[8]= hits[8]+1;
b251a2b5 2570 hits[9]= (Float_t) fNSDigits;
452a64c6 2571 }
2572
2573 //if(sector !=-1)
2574 new(lhits[fNhits++]) AliRICHHit(fIshunt,gAlice->CurrentTrack(),vol,hits);
2575 eloss = 0;
2576 //
2577 // Check additional signal generation conditions
2578 // defined by the segmentation
2579 // model (boundary crossing conditions)
2580 } else if
2581 (((AliRICHChamber*) (*fChambers)[idvol])
2582 ->SigGenCond(localPos[0], localPos[2], localPos[1]))
2583 {
2584 ((AliRICHChamber*) (*fChambers)[idvol])
2585 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2586 if (eloss > 0)
2587 {
2588 if(gMC->TrackPid() == kNeutron)
2589 printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
b251a2b5 2590 nPads = Hits2SDigits(xhit,yhit,eloss,idvol,kMip);
452a64c6 2591 hits[17] = nPads;
2592 //printf("Npads:%d",NPads);
2593 }
2594 xhit = localPos[0];
2595 yhit = localPos[2];
2596 eloss = destep;
2597 tlength += step ;
2598 //
2599 // nothing special happened, add up energy loss
2600 } else {
2601 eloss += destep;
2602 tlength += step ;
2603 }
2604 }
2605 }
2606 /*************************************************End of MIP treatment**************************************/
2607 //}
ddae0931 2608}
2e5f0f7b 2609
237c933d 2610void AliRICH::FindClusters(Int_t nev,Int_t lastEntry)
ddae0931 2611{
2e5f0f7b 2612
ddae0931 2613//
2614// Loop on chambers and on cathode planes
2615//
2e5f0f7b 2616 for (Int_t icat=1;icat<2;icat++) {
2617 gAlice->ResetDigits();
34ead2dd 2618 gAlice->TreeD()->GetEvent(0);
237c933d 2619 for (Int_t ich=0;ich<kNCH;ich++) {
2e5f0f7b 2620 AliRICHChamber* iChamber=(AliRICHChamber*) (*fChambers)[ich];
237c933d 2621 TClonesArray *pRICHdigits = this->DigitsAddress(ich);
2622 if (pRICHdigits == 0)
2e5f0f7b 2623 continue;
2624 //
2625 // Get ready the current chamber stuff
2626 //
2627 AliRICHResponse* response = iChamber->GetResponseModel();
a2f7eaf6 2628 AliSegmentation* seg = iChamber->GetSegmentationModel();
2e5f0f7b 2629 AliRICHClusterFinder* rec = iChamber->GetReconstructionModel();
2630 if (seg) {
2631 rec->SetSegmentation(seg);
2632 rec->SetResponse(response);
237c933d 2633 rec->SetDigits(pRICHdigits);
2e5f0f7b 2634 rec->SetChamber(ich);
2635 if (nev==0) rec->CalibrateCOG();
2636 rec->FindRawClusters();
2637 }
2638 TClonesArray *fRch;
2639 fRch=RawClustAddress(ich);
2640 fRch->Sort();
2641 } // for ich
2642
2643 gAlice->TreeR()->Fill();
2644 TClonesArray *fRch;
237c933d 2645 for (int i=0;i<kNCH;i++) {
2e5f0f7b 2646 fRch=RawClustAddress(i);
2647 int nraw=fRch->GetEntriesFast();
2648 printf ("Chamber %d, raw clusters %d\n",i,nraw);
2649 }
2650
2651 ResetRawClusters();
2652
2653 } // for icat
2654
2655 char hname[30];
2656 sprintf(hname,"TreeR%d",nev);
33984590 2657 gAlice->TreeR()->Write(hname,kOverwrite,0);
2e5f0f7b 2658 gAlice->TreeR()->Reset();
2659
2660 //gObjectTable->Print();
ddae0931 2661}
2662
b251a2b5 2663AliRICHSDigit* AliRICH::FirstPad(AliRICHHit* hit,TClonesArray *clusters )
ddae0931 2664{
2665//
2666 // Initialise the pad iterator
b251a2b5 2667 // Return the address of the first sdigit for hit
ddae0931 2668 TClonesArray *theClusters = clusters;
2669 Int_t nclust = theClusters->GetEntriesFast();
2670 if (nclust && hit->fPHlast > 0) {
2e5f0f7b 2671 sMaxIterPad=Int_t(hit->fPHlast);
2672 sCurIterPad=Int_t(hit->fPHfirst);
b251a2b5 2673 return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1);
ddae0931 2674 } else {
2675 return 0;
fe4da5cc 2676 }
ddae0931 2677
fe4da5cc 2678}
2679
b251a2b5 2680AliRICHSDigit* AliRICH::NextPad(TClonesArray *clusters)
fe4da5cc 2681{
237c933d 2682
2683 // Iterates over pads
2684
ddae0931 2685 sCurIterPad++;
2686 if (sCurIterPad <= sMaxIterPad) {
b251a2b5 2687 return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1);
ddae0931 2688 } else {
2689 return 0;
2690 }
fe4da5cc 2691}
2692
237c933d 2693AliRICH& AliRICH::operator=(const AliRICH& rhs)
2694{
2695// Assignment operator
2696 return *this;
2697
2698}
ddae0931 2699
cc683707 2700void AliRICH::DiagnosticsFE(Int_t evNumber1,Int_t evNumber2)
a3d71079 2701{
2702
2703 Int_t NpadX = 162; // number of pads on X
2704 Int_t NpadY = 162; // number of pads on Y
2705
2706 Int_t Pad[162][162];
2707 for (Int_t i=0;i<NpadX;i++) {
2708 for (Int_t j=0;j<NpadY;j++) {
2709 Pad[i][j]=0;
2710 }
2711 }
2712
2713 // Create some histograms
2714
2715 TH1F *pionspectra1 = new TH1F("pionspectra1","Pion Spectra",200,-4,2);
2716 TH1F *pionspectra2 = new TH1F("pionspectra2","Pion Spectra",200,-4,2);
2717 TH1F *pionspectra3 = new TH1F("pionspectra3","Pion Spectra",200,-4,2);
2718 TH1F *protonspectra1 = new TH1F("protonspectra1","Proton Spectra",200,-4,2);
2719 TH1F *protonspectra2 = new TH1F("protonspectra2","Proton Spectra",200,-4,2);
2720 TH1F *protonspectra3 = new TH1F("protonspectra3","Proton Spectra",200,-4,2);
2721 TH1F *kaonspectra1 = new TH1F("kaonspectra1","Kaon Spectra",100,-4,2);
2722 TH1F *kaonspectra2 = new TH1F("kaonspectra2","Kaon Spectra",100,-4,2);
2723 TH1F *kaonspectra3 = new TH1F("kaonspectra3","Kaon Spectra",100,-4,2);
2724 TH1F *electronspectra1 = new TH1F("electronspectra1","Electron Spectra",100,-4,2);
2725 TH1F *electronspectra2 = new TH1F("electronspectra2","Electron Spectra",100,-4,2);
2726 TH1F *electronspectra3 = new TH1F("electronspectra3","Electron Spectra",100,-4,2);
2727 TH1F *muonspectra1 = new TH1F("muonspectra1","Muon Spectra",100,-4,2);
2728 TH1F *muonspectra2 = new TH1F("muonspectra2","Muon Spectra",100,-4,2);
2729 TH1F *muonspectra3 = new TH1F("muonspectra3","Muon Spectra",100,-4,2);
2730 TH1F *neutronspectra1 = new TH1F("neutronspectra1","Neutron Spectra",100,-4,2);
2731 TH1F *neutronspectra2 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
2732 TH1F *neutronspectra3 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
2733 TH1F *chargedspectra1 = new TH1F("chargedspectra1","Charged particles above 1 GeV Spectra",100,-1,3);
2734 TH1F *chargedspectra2 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
2735 TH1F *chargedspectra3 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
2736 TH1F *pionptspectrafinal = new TH1F("pionptspectrafinal","Primary Pions Transverse Momenta at HMPID",20,0,5);
2737 TH1F *pionptspectravertex = new TH1F("pionptspectravertex","Primary Pions Transverse Momenta at vertex",20,0,5);
2738 TH1F *kaonptspectrafinal = new TH1F("kaonptspectrafinal","Primary Kaons Transverse Momenta at HMPID",20,0,5);
2739 TH1F *kaonptspectravertex = new TH1F("kaonptspectravertex","Primary Kaons Transverse Momenta at vertex",20,0,5);
2740 //TH1F *hitsPhi = new TH1F("hitsPhi","Distribution of phi angle of incidence",100,-180,180);
2741 TH1F *hitsTheta = new TH1F("hitsTheta","Distribution of Theta angle of incidence, all tracks",100,0,50);
2742 TH1F *hitsTheta500MeV = new TH1F("hitsTheta500MeV","Distribution of Theta angle of incidence, 0.5-1 GeV primary tracks",100,0,50);
2743 TH1F *hitsTheta1GeV = new TH1F("hitsTheta1GeV","Distribution of Theta angle of incidence, 1-2 GeV primary tracks",100,0,50);
2744 TH1F *hitsTheta2GeV = new TH1F("hitsTheta2GeV","Distribution of Theta angle of incidence, 2-3 GeV primary tracks",100,0,50);
2745 TH1F *hitsTheta3GeV = new TH1F("hitsTheta3GeV","Distribution of Theta angle of incidence, >3 GeV primary tracks",100,0,50);
2746 TH2F *production = new TH2F("production","Mother production vertices",100,-300,300,100,0,600);
2747
2748
2749
2750
2751// Start loop over events
2752
2753 Int_t pion=0, kaon=0, proton=0, electron=0, positron=0, neutron=0, highneutrons=0, muon=0;
2754 Int_t chargedpions=0,primarypions=0,highprimarypions=0,chargedkaons=0,primarykaons=0,highprimarykaons=0;
2755 Int_t photons=0, primaryphotons=0, highprimaryphotons=0;
2756 TRandom* random=0;
2757
2758 for (int nev=0; nev<= evNumber2; nev++) {
2759 Int_t nparticles = gAlice->GetEvent(nev);
2760
2761
2762 printf ("Event number : %d\n",nev);
2763 printf ("Number of particles: %d\n",nparticles);
2764 if (nev < evNumber1) continue;
2765 if (nparticles <= 0) return;
2766
2767// Get pointers to RICH detector and Hits containers
2768
2769 AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
2770
2771 TTree *treeH = gAlice->TreeH();
2772 Int_t ntracks =(Int_t) treeH->GetEntries();
2773
2774// Start loop on tracks in the hits containers
2775
2776 for (Int_t track=0; track<ntracks;track++) {
2777 printf ("Processing Track: %d\n",track);
2778 gAlice->ResetHits();
2779 treeH->GetEvent(track);
2780
2781 for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
2782 mHit;
2783 mHit=(AliRICHHit*)pRICH->NextHit())
2784 {
2785 //Int_t nch = mHit->fChamber; // chamber number
2786 //Float_t x = mHit->X(); // x-pos of hit
2787 //Float_t y = mHit->Z(); // y-pos
2788 //Float_t z = mHit->Y();
2789 //Float_t phi = mHit->fPhi; //Phi angle of incidence
2790 Float_t theta = mHit->fTheta; //Theta angle of incidence
2791 Float_t px = mHit->MomX();
2792 Float_t py = mHit->MomY();
2793 Int_t index = mHit->Track();
2794 Int_t particle = (Int_t)(mHit->fParticle);
2795 Float_t R;
2796 Float_t PTfinal;
2797 Float_t PTvertex;
2798
2799 TParticle *current = gAlice->Particle(index);
2800
2801 //Float_t energy=current->Energy();
2802
2803 R=TMath::Sqrt(current->Vx()*current->Vx() + current->Vy()*current->Vy());
2804 PTfinal=TMath::Sqrt(px*px + py*py);
2805 PTvertex=TMath::Sqrt(current->Px()*current->Px() + current->Py()*current->Py());
2806
2807
2808
2809 if (TMath::Abs(particle) < 10000000)
2810 {
2811 hitsTheta->Fill(theta,(float) 1);
2812 if (R<5)
2813 {
2814 if (PTvertex>.5 && PTvertex<=1)
2815 {
2816 hitsTheta500MeV->Fill(theta,(float) 1);
2817 }
2818 if (PTvertex>1 && PTvertex<=2)
2819 {
2820 hitsTheta1GeV->Fill(theta,(float) 1);
2821 }
2822 if (PTvertex>2 && PTvertex<=3)
2823 {
2824 hitsTheta2GeV->Fill(theta,(float) 1);
2825 }
2826 if (PTvertex>3)
2827 {
2828 hitsTheta3GeV->Fill(theta,(float) 1);
2829 }
2830 }
2831
2832 }
2833
2834 //if (nch == 3)
2835 //{
2836
2837 //printf("Particle type: %d\n",current->GetPdgCode());
2838 if (TMath::Abs(particle) < 50000051)
2839 {
2840 //if (TMath::Abs(particle) == 50000050 || TMath::Abs(particle) == 2112)
2841 if (TMath::Abs(particle) == 2112 || TMath::Abs(particle) == 50000050)
2842 {
2843 //gMC->Rndm(&random, 1);
2844 if (random->Rndm() < .1)
2845 production->Fill(current->Vz(),R,(float) 1);
2846 if (TMath::Abs(particle) == 50000050)
2847 //if (TMath::Abs(particle) > 50000000)
2848 {
2849 photons +=1;
2850 if (R<5)
2851 {
2852 primaryphotons +=1;
2853 if (current->Energy()>0.001)
2854 highprimaryphotons +=1;
2855 }
2856 }
2857 if (TMath::Abs(particle) == 2112)
2858 {
2859 neutron +=1;
2860 if (current->Energy()>0.0001)
2861 highneutrons +=1;
2862 }
2863 }
2864 if (TMath::Abs(particle) < 50000000)
2865 {
2866 production->Fill(current->Vz(),R,(float) 1);
2867 //printf("Adding %d at %f\n",particle,R);
2868 }
2869 //mip->Fill(x,y,(float) 1);
2870 }
2871
2872 if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
2873 {
2874 if (R<5)
2875 {
2876 pionptspectravertex->Fill(PTvertex,(float) 1);
2877 pionptspectrafinal->Fill(PTfinal,(float) 1);
2878 }
2879 }
2880
2881 if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
2882 || TMath::Abs(particle)==311)
2883 {
2884 if (R<5)
2885 {
2886 kaonptspectravertex->Fill(PTvertex,(float) 1);
2887 kaonptspectrafinal->Fill(PTfinal,(float) 1);
2888 }
2889 }
2890
2891
2892 if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
2893 {
2894 pionspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2895 //printf ("fParticle: %d, PDG code:%d\n",particle,current->GetPdgCode());
2896 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2897 pionspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2898 if (R>250 && R<450)
2899 {
2900 pionspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2901 //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\R:%f\n\n\n\n\n\n\n\n\n",R);
2902 }
2903 //printf("Pion mass: %e\n",current->GetCalcMass());
2904 pion +=1;
2905 if (TMath::Abs(particle)==211)
2906 {
2907 chargedpions +=1;
2908 if (R<5)
2909 {
2910 primarypions +=1;
2911 if (current->Energy()>1)
2912 highprimarypions +=1;
2913 }
2914 }
2915 }
2916 if (TMath::Abs(particle)==2212)
2917 {
2918 protonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2919 //ptspectra->Fill(Pt,(float) 1);
2920 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2921 protonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2922 if (R>250 && R<450)
2923 protonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2924 //printf("\n\n\n\n\n\n\nProton mass: %e\n\n\n\n\n\n\n\n\n",current->GetCalcMass());
2925 proton +=1;
2926 }
2927 if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
2928 || TMath::Abs(particle)==311)
2929 {
2930 kaonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2931 //ptspectra->Fill(Pt,(float) 1);
2932 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2933 kaonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2934 if (R>250 && R<450)
2935 kaonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2936 //printf("Kaon mass: %e\n",current->GetCalcMass());
2937 kaon +=1;
2938 if (TMath::Abs(particle)==321)
2939 {
2940 chargedkaons +=1;
2941 if (R<5)
2942 {
2943 primarykaons +=1;
2944 if (current->Energy()>1)
2945 highprimarykaons +=1;
2946 }
2947 }
2948 }
2949 if (TMath::Abs(particle)==11)
2950 {
2951 electronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2952 //ptspectra->Fill(Pt,(float) 1);
2953 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2954 electronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2955 if (R>250 && R<450)
2956 electronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2957 //printf("Electron mass: %e\n",current->GetCalcMass());
2958 if (particle == 11)
2959 electron +=1;
2960 if (particle == -11)
2961 positron +=1;
2962 }
2963 if (TMath::Abs(particle)==13)
2964 {
2965 muonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2966 //ptspectra->Fill(Pt,(float) 1);
2967 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2968 muonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2969 if (R>250 && R<450)
2970 muonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2971 //printf("Muon mass: %e\n",current->GetCalcMass());
2972 muon +=1;
2973 }
2974 if (TMath::Abs(particle)==2112)
2975 {
2976 neutronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2977 //ptspectra->Fill(Pt,(float) 1);
2978 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2979 neutronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2980 if (R>250 && R<450)
2981 {
2982 neutronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2983 //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\R:%f\n\n\n\n\n\n\n\n\n",R);
2984 }
2985 //printf("Neutron mass: %e\n",current->GetCalcMass());
2986 neutron +=1;
2987 }
2988 if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
2989 {
2990 if (current->Energy()-current->GetCalcMass()>1)
2991 {
2992 chargedspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2993 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2994 chargedspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2995 if (R>250 && R<450)
2996 chargedspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2997 }
2998 }
2999 //printf("Hits:%d\n",hit);
3000 //printf ("Chamber number:%d x:%f y:%f\n",nch,x,y);
3001 // Fill the histograms
3002 //Nh1+=nhits;
3003 //h->Fill(x,y,(float) 1);
3004 //}
3005 //}
3006 }
3007
3008 }
3009
3010 }
3011 // }
3012
3013 //Create canvases, set the view range, show histograms
3014
3015 TCanvas *c2 = new TCanvas("c2","Angles of incidence",150,150,100,150);
3016 c2->Divide(2,2);
3017 //c2->SetFillColor(42);
3018
3019 c2->cd(1);
3020 hitsTheta500MeV->SetFillColor(5);
3021 hitsTheta500MeV->Draw();
3022 c2->cd(2);
3023 hitsTheta1GeV->SetFillColor(5);
3024 hitsTheta1GeV->Draw();
3025 c2->cd(3);
3026 hitsTheta2GeV->SetFillColor(5);
3027 hitsTheta2GeV->Draw();
3028 c2->cd(4);
3029 hitsTheta3GeV->SetFillColor(5);
3030 hitsTheta3GeV->Draw();
3031
3032
3033
3034 TCanvas *c15 = new TCanvas("c15","Mothers Production Vertices",50,50,600,600);
3035 c15->cd();
3036 production->SetFillColor(42);
3037 production->SetXTitle("z (m)");
3038 production->SetYTitle("R (m)");
3039 production->Draw();
3040
3041 TCanvas *c10 = new TCanvas("c10","Pt Spectra",50,50,600,700);
3042 c10->Divide(2,2);
3043 c10->cd(1);
3044 pionptspectravertex->SetFillColor(5);
3045 pionptspectravertex->SetXTitle("Pt (GeV)");
3046 pionptspectravertex->Draw();
3047 c10->cd(2);
3048 pionptspectrafinal->SetFillColor(5);
3049 pionptspectrafinal->SetXTitle("Pt (GeV)");
3050 pionptspectrafinal->Draw();
3051 c10->cd(3);
3052 kaonptspectravertex->SetFillColor(5);
3053 kaonptspectravertex->SetXTitle("Pt (GeV)");
3054 kaonptspectravertex->Draw();
3055 c10->cd(4);
3056 kaonptspectrafinal->SetFillColor(5);
3057 kaonptspectrafinal->SetXTitle("Pt (GeV)");
3058 kaonptspectrafinal->Draw();
3059
3060
3061 TCanvas *c16 = new TCanvas("c16","Particles Spectra II",150,150,600,350);
3062 c16->Divide(2,1);
3063
3064 c16->cd(1);
3065 //TCanvas *c13 = new TCanvas("c13","Electron Spectra",400,10,600,700);
3066 electronspectra1->SetFillColor(5);
3067 electronspectra1->SetXTitle("log(GeV)");
3068 electronspectra2->SetFillColor(46);
3069 electronspectra2->SetXTitle("log(GeV)");
3070 electronspectra3->SetFillColor(10);
3071 electronspectra3->SetXTitle("log(GeV)");
3072 //c13->SetLogx();
3073 electronspectra1->Draw();
3074 electronspectra2->Draw("same");
3075 electronspectra3->Draw("same");
3076
3077 c16->cd(2);
3078 //TCanvas *c14 = new TCanvas("c14","Muon Spectra",400,10,600,700);
3079 muonspectra1->SetFillColor(5);
3080 muonspectra1->SetXTitle("log(GeV)");
3081 muonspectra2->SetFillColor(46);
3082 muonspectra2->SetXTitle("log(GeV)");
3083 muonspectra3->SetFillColor(10);
3084 muonspectra3->SetXTitle("log(GeV)");
3085 //c14->SetLogx();
3086 muonspectra1->Draw();
3087 muonspectra2->Draw("same");
3088 muonspectra3->Draw("same");
3089
3090 //c16->cd(3);
3091 //TCanvas *c16 = new TCanvas("c16","Neutron Spectra",400,10,600,700);
3092 //neutronspectra1->SetFillColor(42);
3093 //neutronspectra1->SetXTitle("log(GeV)");
3094 //neutronspectra2->SetFillColor(46);
3095 //neutronspectra2->SetXTitle("log(GeV)");
3096 //neutronspectra3->SetFillColor(10);
3097 //neutronspectra3->SetXTitle("log(GeV)");
3098 //c16->SetLogx();
3099 //neutronspectra1->Draw();
3100 //neutronspectra2->Draw("same");
3101 //neutronspectra3->Draw("same");
3102
3103 TCanvas *c9 = new TCanvas("c9","Particles Spectra",150,150,600,700);
3104 //TCanvas *c9 = new TCanvas("c9","Pion Spectra",400,10,600,700);
3105 c9->Divide(2,2);
3106
3107 c9->cd(1);
3108 pionspectra1->SetFillColor(5);
3109 pionspectra1->SetXTitle("log(GeV)");
3110 pionspectra2->SetFillColor(46);
3111 pionspectra2->SetXTitle("log(GeV)");
3112 pionspectra3->SetFillColor(10);
3113 pionspectra3->SetXTitle("log(GeV)");
3114 //c9->SetLogx();
3115 pionspectra1->Draw();
3116 pionspectra2->Draw("same");
3117 pionspectra3->Draw("same");
3118
3119 c9->cd(2);
3120 //TCanvas *c10 = new TCanvas("c10","Proton Spectra",400,10,600,700);
3121 protonspectra1->SetFillColor(5);
3122 protonspectra1->SetXTitle("log(GeV)");
3123 protonspectra2->SetFillColor(46);
3124 protonspectra2->SetXTitle("log(GeV)");
3125 protonspectra3->SetFillColor(10);
3126 protonspectra3->SetXTitle("log(GeV)");
3127 //c10->SetLogx();
3128 protonspectra1->Draw();
3129 protonspectra2->Draw("same");
3130 protonspectra3->Draw("same");
3131
3132 c9->cd(3);
3133 //TCanvas *c11 = new TCanvas("c11","Kaon Spectra",400,10,600,700);
3134 kaonspectra1->SetFillColor(5);
3135 kaonspectra1->SetXTitle("log(GeV)");
3136 kaonspectra2->SetFillColor(46);
3137 kaonspectra2->SetXTitle("log(GeV)");
3138 kaonspectra3->SetFillColor(10);
3139 kaonspectra3->SetXTitle("log(GeV)");
3140 //c11->SetLogx();
3141 kaonspectra1->Draw();
3142 kaonspectra2->Draw("same");
3143 kaonspectra3->Draw("same");
3144
3145 c9->cd(4);
3146 //TCanvas *c12 = new TCanvas("c12","Charged Particles Spectra",400,10,600,700);
3147 chargedspectra1->SetFillColor(5);
3148 chargedspectra1->SetXTitle("log(GeV)");
3149 chargedspectra2->SetFillColor(46);
3150 chargedspectra2->SetXTitle("log(GeV)");
3151 chargedspectra3->SetFillColor(10);
3152 chargedspectra3->SetXTitle("log(GeV)");
3153 //c12->SetLogx();
3154 chargedspectra1->Draw();
3155 chargedspectra2->Draw("same");
3156 chargedspectra3->Draw("same");
3157
3158
3159
3160 printf("*****************************************\n");
3161 printf("* Particle * Counts *\n");
3162 printf("*****************************************\n");
3163
3164 printf("* Pions: * %4d *\n",pion);
3165 printf("* Charged Pions: * %4d *\n",chargedpions);
3166 printf("* Primary Pions: * %4d *\n",primarypions);
3167 printf("* Primary Pions (p>1GeV/c): * %4d *\n",highprimarypions);
3168 printf("* Kaons: * %4d *\n",kaon);
3169 printf("* Charged Kaons: * %4d *\n",chargedkaons);
3170 printf("* Primary Kaons: * %4d *\n",primarykaons);
3171 printf("* Primary Kaons (p>1GeV/c): * %4d *\n",highprimarykaons);
3172 printf("* Muons: * %4d *\n",muon);
3173 printf("* Electrons: * %4d *\n",electron);
3174 printf("* Positrons: * %4d *\n",positron);
3175 printf("* Protons: * %4d *\n",proton);
3176 printf("* All Charged: * %4d *\n",(chargedpions+chargedkaons+muon+electron+positron+proton));
3177 printf("*****************************************\n");
3178 //printf("* Photons: * %3.1f *\n",photons);
3179 //printf("* Primary Photons: * %3.1f *\n",primaryphotons);
3180 //printf("* Primary Photons (p>1MeV/c):* %3.1f *\n",highprimaryphotons);
3181 //printf("*****************************************\n");
3182 //printf("* Neutrons: * %3.1f *\n",neutron);
3183 //printf("* Neutrons (p>100keV/c): * %3.1f *\n",highneutrons);
3184 //printf("*****************************************\n");
3185
3186 if (gAlice->TreeD())
3187 {
3188 gAlice->TreeD()->GetEvent(0);
3189
3190 Float_t occ[7];
3191 Float_t sum=0;
3192 Float_t mean=0;
3193 printf("\n*****************************************\n");
3194 printf("* Chamber * Digits * Occupancy *\n");
3195 printf("*****************************************\n");
3196
3197 for (Int_t ich=0;ich<7;ich++)
3198 {
3199 TClonesArray *Digits = DigitsAddress(ich); // Raw clusters branch
3200 Int_t ndigits = Digits->GetEntriesFast();
3201 occ[ich] = Float_t(ndigits)/(160*144);
3202 sum += Float_t(ndigits)/(160*144);
3203 printf("* %d * %d * %3.1f%% *\n",ich,ndigits,occ[ich]*100);
3204 }
3205 mean = sum/7;
3206 printf("*****************************************\n");
3207 printf("* Mean occupancy * %3.1f%% *\n",mean*100);
3208 printf("*****************************************\n");
3209 }
3210
3211 printf("\nEnd of analysis\n");
3212
3213}
3214
3215//_________________________________________________________________________________________________
3216
3217
cc683707 3218void AliRICH::DiagnosticsSE(Int_t diaglevel,Int_t evNumber1,Int_t evNumber2)
a3d71079 3219{
3220
3221AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH");
3222 AliRICHSegmentationV0* segmentation;
3223 AliRICHChamber* chamber;
3224
3225 chamber = &(pRICH->Chamber(0));
3226 segmentation=(AliRICHSegmentationV0*) chamber->GetSegmentationModel(0);
3227
3228 Int_t NpadX = segmentation->Npx(); // number of pads on X
3229 Int_t NpadY = segmentation->Npy(); // number of pads on Y
3230
3231 //Int_t Pad[144][160];
3232 /*for (Int_t i=0;i<NpadX;i++) {
3233 for (Int_t j=0;j<NpadY;j++) {
3234 Pad[i][j]=0;
3235 }
3236 } */
3237
3238
3239 Int_t xmin= -NpadX/2;
3240 Int_t xmax= NpadX/2;
3241 Int_t ymin= -NpadY/2;
3242 Int_t ymax= NpadY/2;
3243
3244 TH2F *feedback = 0;
3245 TH2F *mip = 0;
3246 TH2F *cerenkov = 0;
3247 TH2F *h = 0;
3248 TH1F *hitsX = 0;
3249 TH1F *hitsY = 0;
3250
3251 TH2F *hc0 = new TH2F("hc0","Zoom on center of central chamber",150,-30,30,150,-50,10);
3252
3253 if (diaglevel == 1)
3254 {
3255 printf("Single Ring Hits\n");
3256 feedback = new TH2F("feedback","Feedback hit distribution",150,-30,30,150,-50,10);
3257 mip = new TH2F("mip","Mip hit distribution",150,-30,30,150,-50,10);
3258 cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-30,30,150,-50,10);
3259 h = new TH2F("h","Detector hit distribution",150,-30,30,150,-50,10);
3260 hitsX = new TH1F("hitsX","Distribution of hits along x-axis",150,-30,30);
3261 hitsY = new TH1F("hitsY","Distribution of hits along z-axis",150,-50,10);
3262 }
3263 else
3264 {
3265 printf("Full Event Hits\n");
3266
3267 feedback = new TH2F("feedback","Feedback hit distribution",150,-300,300,150,-300,300);
3268 mip = new TH2F("mip","Mip hit distribution",150,-300,300,150,-300,300);
3269 cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-300,300,150,-300,300);
3270 h = new TH2F("h","Detector hit distribution",150,-300,300,150,-300,300);
3271 hitsX = new TH1F("digitsX","Distribution of hits along x-axis",200,-300,300);
3272 hitsY = new TH1F("digitsY","Distribution of hits along z-axis",200,-300,300);
3273 }
3274
ca96c9ea 3275
ddae0931 3276
a3d71079 3277 TH2F *hc1 = new TH2F("hc1","Chamber 1 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3278 TH2F *hc2 = new TH2F("hc2","Chamber 2 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3279 TH2F *hc3 = new TH2F("hc3","Chamber 3 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3280 TH2F *hc4 = new TH2F("hc4","Chamber 4 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3281 TH2F *hc5 = new TH2F("hc5","Chamber 5 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3282 TH2F *hc6 = new TH2F("hc6","Chamber 6 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3283 TH2F *hc7 = new TH2F("hc7","Chamber 7 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3284
3285 TH1F *Clcharge = new TH1F("Clcharge","Cluster Charge Distribution",500,0.,500.);
3286 TH1F *ckovangle = new TH1F("ckovangle","Cerenkov angle per photon",200,.3,1);
3287 TH1F *hckphi = new TH1F("hckphi","Cerenkov phi angle per photon",620,-3.1,3.1);
3288 TH1F *mother = new TH1F("mother","Cerenkovs per Mip",75,0.,75.);
3289 TH1F *radius = new TH1F("radius","Mean distance to Mip",100,0.,20.);
3290 TH1F *phspectra1 = new TH1F("phspectra1","Detected Photon Spectra",200,5.,10.);
3291 TH1F *phspectra2 = new TH1F("phspectra2","Produced Photon Spectra",200,5.,10.);
3292 TH1F *totalphotonstrack = new TH1F("totalphotonstrack","Produced Photons per Mip",100,200,700.);
3293 TH1F *totalphotonsevent = new TH1F("totalphotonsevent","Produced Photons per Mip",100,200,700.);
3294 //TH1F *feedbacks = new TH1F("feedbacks","Produced Feedbacks per Mip",50,0.5,50.);
3295 TH1F *padnumber = new TH1F("padnumber","Number of pads per cluster",50,-0.5,50.);
3296 TH1F *padsev = new TH1F("padsev","Number of pads hit per MIP",50,0.5,100.);
3297 TH1F *clusev = new TH1F("clusev","Number of clusters per MIP",50,0.5,50.);
3298 TH1F *photev = new TH1F("photev","Number of detected photons per MIP",50,0.5,50.);
3299 TH1F *feedev = new TH1F("feedev","Number of feedbacks per MIP",50,0.5,50.);
3300 TH1F *padsmip = new TH1F("padsmip","Number of pads per event inside MIP region",50,0.5,50.);
3301 TH1F *padscl = new TH1F("padscl","Number of pads per event from cluster count",50,0.5,100.);
3302 TH1F *pionspectra = new TH1F("pionspectra","Pion Spectra",200,.5,10.);
3303 TH1F *protonspectra = new TH1F("protonspectra","Proton Spectra",200,.5,10.);
3304 TH1F *kaonspectra = new TH1F("kaonspectra","Kaon Spectra",100,.5,10.);
3305 TH1F *chargedspectra = new TH1F("chargedspectra","Charged particles above 1 GeV Spectra",100,.5,10.);
3306 TH1F *hitsPhi = new TH1F("hitsPhi","Distribution of phi angle of incidence",100,0,360);
3307 TH1F *hitsTheta = new TH1F("hitsTheta","Distribution of Theta angle of incidence",100,0,15);
3308 TH1F *Omega1D = new TH1F("omega","Reconstructed Cerenkov angle per track",200,.5,1);
3309 TH1F *Theta = new TH1F("theta","Reconstructed theta incidence angle per track",200,0,15);
3310 TH1F *Phi = new TH1F("phi","Reconstructed phi incidence per track",200,0,360);
3311 TH1F *Omega3D = new TH1F("omega","Reconstructed Cerenkov angle per track",200,.3,1);
3312 TH1F *PhotonCer = new TH1F("photoncer","Reconstructed Cerenkov angle per photon",200,.3,1);
3313 TH2F *PadsUsed = new TH2F("padsused","Pads Used for Reconstruction",100,-30,30,100,-30,30);
3314 TH1F *MeanRadius = new TH1F("radius","Mean Radius for reconstructed track",100,0.,20.);
3315
3316// Start loop over events
3317
3318 Int_t Nh=0;
3319 Int_t pads=0;
3320 Int_t Nh1=0;
3321 Int_t mothers[80000];
3322 Int_t mothers2[80000];
3323 Float_t mom[3];
3324 Int_t nraw=0;
3325 Int_t phot=0;
3326 Int_t feed=0;
3327 Int_t padmip=0;
3328 Float_t x=0,y=0;
3329
3330 for (Int_t i=0;i<100;i++) mothers[i]=0;
3331
3332 for (int nev=0; nev<= evNumber2; nev++) {
3333 Int_t nparticles = gAlice->GetEvent(nev);
3334
3335
3336 //cout<<"nev "<<nev<<endl;
3337 printf ("\n**********************************\nProcessing Event: %d\n",nev);
3338 //cout<<"nparticles "<<nparticles<<endl;
3339 printf ("Particles : %d\n\n",nparticles);
3340 if (nev < evNumber1) continue;
3341 if (nparticles <= 0) return;
3342
3343// Get pointers to RICH detector and Hits containers
3344
3345
3346 TTree *TH = gAlice->TreeH();
3347 Stat_t ntracks = TH->GetEntries();
3348
3349 // Start loop on tracks in the hits containers
3350 //Int_t Nc=0;
3351 for (Int_t track=0; track<ntracks;track++) {
3352
3353 printf ("\nProcessing Track: %d\n",track);
3354 gAlice->ResetHits();
3355 TH->GetEvent(track);
3356 Int_t nhits = pRICH->Hits()->GetEntriesFast();
3357 if (nhits) Nh+=nhits;
3358 printf("Hits : %d\n",nhits);
3359 for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
3360 mHit;
3361 mHit=(AliRICHHit*)pRICH->NextHit())
3362 {
3363 //Int_t nch = mHit->fChamber; // chamber number
3364 x = mHit->X(); // x-pos of hit
3365 y = mHit->Z(); // y-pos
3366 Float_t phi = mHit->fPhi; //Phi angle of incidence
3367 Float_t theta = mHit->fTheta; //Theta angle of incidence
3368 Int_t index = mHit->Track();
3369 Int_t particle = (Int_t)(mHit->fParticle);
3370 //Int_t freon = (Int_t)(mHit->fLoss);
3371
3372 hitsX->Fill(x,(float) 1);
3373 hitsY->Fill(y,(float) 1);
3374
3375 //printf("Particle:%9d\n",particle);
3376
3377 TParticle *current = (TParticle*)gAlice->Particle(index);
3378 //printf("Particle type: %d\n",sizeoff(Particles));
3379
3380 hitsTheta->Fill(theta,(float) 1);
3381 //hitsPhi->Fill(phi,(float) 1);
3382 //if (pRICH->GetDebugLevel() == -1)
3383 //printf("Theta:%f, Phi:%f\n",theta,phi);
3384
3385 //printf("Debug Level:%d\n",pRICH->GetDebugLevel());
3386
3387 if (current->GetPdgCode() < 10000000)
3388 {
3389 mip->Fill(x,y,(float) 1);
3390 //printf("adding mip\n");
3391 //if (current->Energy() - current->GetCalcMass()>1 && freon==1)
3392 //{
3393 hitsPhi->Fill(TMath::Abs(phi),(float) 1);
3394 //hitsTheta->Fill(theta,(float) 1);
3395 //printf("Theta:%f, Phi:%f\n",theta,phi);
3396 //}
3397 }
3398
3399 if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
3400 {
3401 pionspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
3402 }
3403 if (TMath::Abs(particle)==2212)
3404 {
3405 protonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
3406 }
3407 if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
3408 || TMath::Abs(particle)==311)
3409 {
3410 kaonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
3411 }
3412 if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
3413 {
3414 if (current->Energy() - current->GetCalcMass()>1)
3415 chargedspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
3416 }