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