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