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