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