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