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