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