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