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