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