]> 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
Corrected invertes arguments in pad size setter
[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$
ee2105a4 18 Revision 1.37 2000/12/20 14:07:25 jbarbosa
19 Removed dependencies on TGeant3 (thanks to F. Carminati and I. Hrivnacova)
20
5cfcdf54 21 Revision 1.36 2000/12/18 17:45:54 jbarbosa
22 Cleaned up PadHits object.
23
176a9917 24 Revision 1.35 2000/12/15 16:49:40 jbarbosa
25 Geometry and materials updates (wire supports, pcbs, backplane supports, frame).
26
6197ac87 27 Revision 1.34 2000/11/10 18:12:12 jbarbosa
28 Bug fix for AliRICHCerenkov (thanks to P. Hristov)
29
3cc8edee 30 Revision 1.33 2000/11/02 10:09:01 jbarbosa
31 Minor bug correction (some pointers were not initialised in the default constructor)
32
9b674b76 33 Revision 1.32 2000/11/01 15:32:55 jbarbosa
34 Updated to handle both reconstruction algorithms.
35
a4622d0f 36 Revision 1.31 2000/10/26 20:18:33 jbarbosa
37 Supports for methane and freon vessels
38
683b273b 39 Revision 1.30 2000/10/24 13:19:12 jbarbosa
40 Geometry updates.
41
e293afae 42 Revision 1.29 2000/10/19 19:39:25 jbarbosa
43 Some more changes to geometry. Further correction of digitisation "per part. type"
44
53a60579 45 Revision 1.28 2000/10/17 20:50:57 jbarbosa
46 Inversed digtise by particle type (now, only the selected particle type is not digitsed).
47 Corrected several geometry minor bugs.
48 Added new parameter (opaque quartz thickness).
49
3587e146 50 Revision 1.27 2000/10/11 10:33:55 jbarbosa
51 Corrected bug introduced by earlier revisions (CerenkovData array cannot be reset to zero on wach call of StepManager)
52
8cda28a5 53 Revision 1.26 2000/10/03 21:44:08 morsch
54 Use AliSegmentation and AliHit abstract base classes.
55
a2f7eaf6 56 Revision 1.25 2000/10/02 21:28:12 fca
57 Removal of useless dependecies via forward declarations
58
94de3818 59 Revision 1.24 2000/10/02 15:43:17 jbarbosa
60 Fixed forward declarations.
61 Fixed honeycomb density.
62 Fixed cerenkov storing.
63 New electronics.
64
9dda3582 65 Revision 1.23 2000/09/13 10:42:14 hristov
66 Minor corrections for HP, DEC and Sun; strings.h included
67
e813258a 68 Revision 1.22 2000/09/12 18:11:13 fca
69 zero hits area before using
70
e51a3fa9 71 Revision 1.21 2000/07/21 10:21:07 morsch
72 fNrawch = 0; and fNrechits = 0; in the default constructor.
73
db481a6e 74 Revision 1.20 2000/07/10 15:28:39 fca
75 Correction of the inheritance scheme
76
452a64c6 77 Revision 1.19 2000/06/30 16:29:51 dibari
78 Added kDebugLevel variable to control output size on demand
79
33984590 80 Revision 1.18 2000/06/12 15:15:46 jbarbosa
81 Cleaned up version.
82
237c933d 83 Revision 1.17 2000/06/09 14:58:37 jbarbosa
84 New digitisation per particle type
85
d28b5fc2 86 Revision 1.16 2000/04/19 12:55:43 morsch
87 Newly structured and updated version (JB, AM)
88
4c039060 89*/
90
2e5f0f7b 91
ddae0931 92////////////////////////////////////////////////
93// Manager and hits classes for set:RICH //
94////////////////////////////////////////////////
fe4da5cc 95
ddae0931 96#include <TBRIK.h>
97#include <TTUBE.h>
fe4da5cc 98#include <TNode.h>
99#include <TRandom.h>
ddae0931 100#include <TObject.h>
101#include <TVector.h>
102#include <TObjArray.h>
2e5f0f7b 103#include <TArrayF.h>
237c933d 104#include <TFile.h>
105#include <TParticle.h>
94de3818 106#include <TGeometry.h>
107#include <TTree.h>
108
237c933d 109#include <iostream.h>
e813258a 110#include <strings.h>
fe4da5cc 111
fe4da5cc 112#include "AliRICH.h"
a2f7eaf6 113#include "AliSegmentation.h"
237c933d 114#include "AliRICHHit.h"
115#include "AliRICHCerenkov.h"
116#include "AliRICHPadHit.h"
117#include "AliRICHDigit.h"
118#include "AliRICHTransientDigit.h"
119#include "AliRICHRawCluster.h"
a4622d0f 120#include "AliRICHRecHit1D.h"
121#include "AliRICHRecHit3D.h"
237c933d 122#include "AliRICHHitMapA1.h"
2e5f0f7b 123#include "AliRICHClusterFinder.h"
fe4da5cc 124#include "AliRun.h"
ddae0931 125#include "AliMC.h"
94de3818 126#include "AliMagF.h"
452a64c6 127#include "AliConst.h"
128#include "AliPDG.h"
2e5f0f7b 129#include "AliPoints.h"
ddae0931 130#include "AliCallf77.h"
237c933d 131
ddae0931 132
133// Static variables for the pad-hit iterator routines
134static Int_t sMaxIterPad=0;
135static Int_t sCurIterPad=0;
136static TClonesArray *fClusters2;
137static TClonesArray *fHits2;
2e5f0f7b 138static TTree *TrH1;
ddae0931 139
fe4da5cc 140ClassImp(AliRICH)
ddae0931 141
142//___________________________________________
fe4da5cc 143AliRICH::AliRICH()
144{
237c933d 145// Default constructor for RICH manager class
146
ddae0931 147 fIshunt = 0;
148 fHits = 0;
2e5f0f7b 149 fPadHits = 0;
150 fNPadHits = 0;
151 fNcerenkovs = 0;
ddae0931 152 fDchambers = 0;
9b674b76 153 fRecHits1D = 0;
154 fRecHits3D = 0;
155 fRawClusters = 0;
156 fChambers = 0;
ddae0931 157 fCerenkovs = 0;
9dda3582 158 for (Int_t i=0; i<7; i++)
159 {
160 fNdch[i] = 0;
161 fNrawch[i] = 0;
a4622d0f 162 fNrechits1D[i] = 0;
163 fNrechits3D[i] = 0;
9dda3582 164 }
9b674b76 165
166 fFileName = 0;
fe4da5cc 167}
ddae0931 168
169//___________________________________________
fe4da5cc 170AliRICH::AliRICH(const char *name, const char *title)
ddae0931 171 : AliDetector(name,title)
fe4da5cc 172{
ddae0931 173//Begin_Html
174/*
175 <img src="gif/alirich.gif">
176*/
177//End_Html
178
2e5f0f7b 179 fHits = new TClonesArray("AliRICHHit",1000 );
1cedd08a 180 gAlice->AddHitList(fHits);
2e5f0f7b 181 fPadHits = new TClonesArray("AliRICHPadHit",100000);
182 fCerenkovs = new TClonesArray("AliRICHCerenkov",1000);
183 gAlice->AddHitList(fCerenkovs);
184 //gAlice->AddHitList(fHits);
185 fNPadHits = 0;
186 fNcerenkovs = 0;
187 fIshunt = 0;
ddae0931 188
9dda3582 189 //fNdch = new Int_t[kNCH];
ddae0931 190
237c933d 191 fDchambers = new TObjArray(kNCH);
2e5f0f7b 192
a4622d0f 193 fRecHits1D = new TObjArray(kNCH);
194 fRecHits3D = new TObjArray(kNCH);
ddae0931 195
196 Int_t i;
197
237c933d 198 for (i=0; i<kNCH ;i++) {
2e5f0f7b 199 (*fDchambers)[i] = new TClonesArray("AliRICHDigit",10000);
ddae0931 200 fNdch[i]=0;
201 }
2e5f0f7b 202
9dda3582 203 //fNrawch = new Int_t[kNCH];
ddae0931 204
237c933d 205 fRawClusters = new TObjArray(kNCH);
d28b5fc2 206 //printf("Created fRwClusters with adress:%p",fRawClusters);
2e5f0f7b 207
237c933d 208 for (i=0; i<kNCH ;i++) {
2e5f0f7b 209 (*fRawClusters)[i] = new TClonesArray("AliRICHRawCluster",10000);
210 fNrawch[i]=0;
211 }
212
9dda3582 213 //fNrechits = new Int_t[kNCH];
ddae0931 214
237c933d 215 for (i=0; i<kNCH ;i++) {
a4622d0f 216 (*fRecHits1D)[i] = new TClonesArray("AliRICHRecHit1D",1000);
217 }
218 for (i=0; i<kNCH ;i++) {
219 (*fRecHits3D)[i] = new TClonesArray("AliRICHRecHit3D",1000);
2e5f0f7b 220 }
d28b5fc2 221 //printf("Created fRecHits with adress:%p",fRecHits);
2e5f0f7b 222
223
ddae0931 224 SetMarkerColor(kRed);
9b674b76 225
3cc8edee 226 /*fChambers = new TObjArray(kNCH);
9b674b76 227 for (i=0; i<kNCH; i++)
3cc8edee 228 (*fChambers)[i] = new AliRICHChamber();*/
9b674b76 229
230 fFileName = 0;
231
fe4da5cc 232}
233
237c933d 234AliRICH::AliRICH(const AliRICH& RICH)
235{
236// Copy Constructor
237}
238
239
ddae0931 240//___________________________________________
fe4da5cc 241AliRICH::~AliRICH()
242{
237c933d 243
244// Destructor of RICH manager class
245
ddae0931 246 fIshunt = 0;
247 delete fHits;
2e5f0f7b 248 delete fPadHits;
ddae0931 249 delete fCerenkovs;
ee2105a4 250
251 //PH Delete TObjArrays
252 if (fChambers) {
253 fChambers->Delete();
254 delete fChambers;
255 }
256 if (fDchambers) {
257 fDchambers->Delete();
258 delete fDchambers;
259 }
260 if (fRawClusters) {
261 fRawClusters->Delete();
262 delete fRawClusters;
263 }
264 if (fRecHits1D) {
265 fRecHits1D->Delete();
266 delete fRecHits1D;
267 }
268 if (fRecHits3D) {
269 fRecHits3D->Delete();
270 delete fRecHits3D;
271 }
272
fe4da5cc 273}
274
ddae0931 275//___________________________________________
fe4da5cc 276void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits)
277{
237c933d 278
279//
280// Adds a hit to the Hits list
281//
282
ddae0931 283 TClonesArray &lhits = *fHits;
2e5f0f7b 284 new(lhits[fNhits++]) AliRICHHit(fIshunt,track,vol,hits);
fe4da5cc 285}
fe4da5cc 286//_____________________________________________________________________________
ddae0931 287void AliRICH::AddCerenkov(Int_t track, Int_t *vol, Float_t *cerenkovs)
288{
237c933d 289
290//
291// Adds a RICH cerenkov hit to the Cerenkov Hits list
292//
293
ddae0931 294 TClonesArray &lcerenkovs = *fCerenkovs;
295 new(lcerenkovs[fNcerenkovs++]) AliRICHCerenkov(fIshunt,track,vol,cerenkovs);
2e5f0f7b 296 //printf ("Done for Cerenkov %d\n\n\n\n",fNcerenkovs);
fe4da5cc 297}
ddae0931 298//___________________________________________
2e5f0f7b 299void AliRICH::AddPadHit(Int_t *clhits)
ddae0931 300{
237c933d 301
302//
303// Add a RICH pad hit to the list
304//
305
2e5f0f7b 306 TClonesArray &lPadHits = *fPadHits;
307 new(lPadHits[fNPadHits++]) AliRICHPadHit(clhits);
33984590 308}
fe4da5cc 309//_____________________________________________________________________________
ddae0931 310void AliRICH::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
311{
237c933d 312
313 //
314 // Add a RICH digit to the list
315 //
ddae0931 316
317 TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]);
2e5f0f7b 318 new(ldigits[fNdch[id]++]) AliRICHDigit(tracks,charges,digits);
fe4da5cc 319}
320
321//_____________________________________________________________________________
2e5f0f7b 322void AliRICH::AddRawCluster(Int_t id, const AliRICHRawCluster& c)
fe4da5cc 323{
ddae0931 324 //
2e5f0f7b 325 // Add a RICH digit to the list
ddae0931 326 //
2e5f0f7b 327
328 TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]);
329 new(lrawcl[fNrawch[id]++]) AliRICHRawCluster(c);
fe4da5cc 330}
331
2e5f0f7b 332//_____________________________________________________________________________
a4622d0f 333void AliRICH::AddRecHit1D(Int_t id, Float_t *rechit, Float_t *photons, Int_t *padsx, Int_t* padsy)
334{
335
336 //
337 // Add a RICH reconstructed hit to the list
338 //
339
340 TClonesArray &lrec1D = *((TClonesArray*)(*fRecHits1D)[id]);
341 new(lrec1D[fNrechits1D[id]++]) AliRICHRecHit1D(id,rechit,photons,padsx,padsy);
342}
343
344//_____________________________________________________________________________
345void AliRICH::AddRecHit3D(Int_t id, Float_t *rechit)
2e5f0f7b 346{
237c933d 347
348 //
349 // Add a RICH reconstructed hit to the list
350 //
fe4da5cc 351
a4622d0f 352 TClonesArray &lrec3D = *((TClonesArray*)(*fRecHits3D)[id]);
353 new(lrec3D[fNrechits3D[id]++]) AliRICHRecHit3D(id,rechit);
2e5f0f7b 354}
fe4da5cc 355
ddae0931 356//___________________________________________
357void AliRICH::BuildGeometry()
358
fe4da5cc 359{
237c933d 360
361 //
362 // Builds a TNode geometry for event display
363 //
364 TNode *node, *top;
ddae0931 365
366 const int kColorRICH = kGreen;
367 //
237c933d 368 top=gAlice->GetGeometry()->GetNode("alice");
ddae0931 369
370
371 new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
372
237c933d 373 top->cd();
ddae0931 374 Float_t pos1[3]={0,471.8999,165.2599};
2e5f0f7b 375 //Chamber(0).SetChamberTransform(pos1[0],pos1[1],pos1[2],
376 new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90);
237c933d 377 node = new TNode("RICH1","RICH1","S_RICH",pos1[0],pos1[1],pos1[2],"rot993");
ddae0931 378
379
237c933d 380 node->SetLineColor(kColorRICH);
381 fNodes->Add(node);
382 top->cd();
ddae0931 383
384 Float_t pos2[3]={171,470,0};
2e5f0f7b 385 //Chamber(1).SetChamberTransform(pos2[0],pos2[1],pos2[2],
386 new TRotMatrix("rot994","rot994",90,-20,90,70,0,0);
237c933d 387 node = new TNode("RICH2","RICH2","S_RICH",pos2[0],pos2[1],pos2[2],"rot994");
ddae0931 388
389
237c933d 390 node->SetLineColor(kColorRICH);
391 fNodes->Add(node);
392 top->cd();
ddae0931 393 Float_t pos3[3]={0,500,0};
2e5f0f7b 394 //Chamber(2).SetChamberTransform(pos3[0],pos3[1],pos3[2],
395 new TRotMatrix("rot995","rot995",90,0,90,90,0,0);
237c933d 396 node = new TNode("RICH3","RICH3","S_RICH",pos3[0],pos3[1],pos3[2],"rot995");
ddae0931 397
398
237c933d 399 node->SetLineColor(kColorRICH);
400 fNodes->Add(node);
401 top->cd();
ddae0931 402 Float_t pos4[3]={-171,470,0};
2e5f0f7b 403 //Chamber(3).SetChamberTransform(pos4[0],pos4[1],pos4[2],
404 new TRotMatrix("rot996","rot996",90,20,90,110,0,0);
237c933d 405 node = new TNode("RICH4","RICH4","S_RICH",pos4[0],pos4[1],pos4[2],"rot996");
ddae0931 406
407
237c933d 408 node->SetLineColor(kColorRICH);
409 fNodes->Add(node);
410 top->cd();
ddae0931 411 Float_t pos5[3]={161.3999,443.3999,-165.3};
2e5f0f7b 412 //Chamber(4).SetChamberTransform(pos5[0],pos5[1],pos5[2],
413 new TRotMatrix("rot997","rot997",90,340,108.1999,70,18.2,70);
237c933d 414 node = new TNode("RICH5","RICH5","S_RICH",pos5[0],pos5[1],pos5[2],"rot997");
ddae0931 415
237c933d 416 node->SetLineColor(kColorRICH);
417 fNodes->Add(node);
418 top->cd();
ddae0931 419 Float_t pos6[3]={0., 471.9, -165.3,};
2e5f0f7b 420 //Chamber(5).SetChamberTransform(pos6[0],pos6[1],pos6[2],
421 new TRotMatrix("rot998","rot998",90,0,109.3099,90,19.30999,90);
237c933d 422 node = new TNode("RICH6","RICH6","S_RICH",pos6[0],pos6[1],pos6[2],"rot998");
ddae0931 423
424
237c933d 425 node->SetLineColor(kColorRICH);
426 fNodes->Add(node);
427 top->cd();
ddae0931 428 Float_t pos7[3]={-161.399,443.3999,-165.3};
2e5f0f7b 429 //Chamber(6).SetChamberTransform(pos7[0],pos7[1],pos7[2],
430 new TRotMatrix("rot999","rot999",90,20,108.1999,110,18.2,110);
237c933d 431 node = new TNode("RICH7","RICH7","S_RICH",pos7[0],pos7[1],pos7[2],"rot999");
432 node->SetLineColor(kColorRICH);
433 fNodes->Add(node);
ddae0931 434
fe4da5cc 435}
436
452a64c6 437//___________________________________________
438void AliRICH::CreateGeometry()
439{
440 //
441 // Create the geometry for RICH version 1
442 //
443 // Modified by: N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
444 // R.A. Fini (INFN - BARI, Rosanna.Fini@ba.infn.it)
445 // R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
446 //
447 //Begin_Html
448 /*
449 <img src="picts/AliRICHv1.gif">
450 */
451 //End_Html
452 //Begin_Html
453 /*
454 <img src="picts/AliRICHv1Tree.gif">
455 */
456 //End_Html
457
458 AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
a2f7eaf6 459 AliSegmentation* segmentation;
452a64c6 460 AliRICHGeometry* geometry;
461 AliRICHChamber* iChamber;
462
463 iChamber = &(pRICH->Chamber(0));
464 segmentation=iChamber->GetSegmentationModel(0);
465 geometry=iChamber->GetGeometryModel();
466
467 Float_t distance;
468 distance = geometry->GetFreonThickness()/2 + geometry->GetQuartzThickness() + geometry->GetGapThickness();
469 geometry->SetRadiatorToPads(distance);
470
3587e146 471 //Opaque quartz thickness
e293afae 472 Float_t oqua_thickness = .5;
683b273b 473 //CsI dimensions
474
475 Float_t csi_length = 160*.8 + 2.6;
476 Float_t csi_width = 144*.84 + 2*2.6;
452a64c6 477
478 Int_t *idtmed = fIdtmed->GetArray()-999;
479
480 Int_t i;
481 Float_t zs;
482 Int_t idrotm[1099];
483 Float_t par[3];
484
485 // --- Define the RICH detector
486 // External aluminium box
e293afae 487 par[0] = 68.8;
6197ac87 488 par[1] = 13; //Original Settings
e293afae 489 par[2] = 70.86;
452a64c6 490 /*par[0] = 73.15;
491 par[1] = 11.5;
492 par[2] = 71.1;*/
493 gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
494
683b273b 495 // Air
496 par[0] = 66.3;
6197ac87 497 par[1] = 13; //Original Settings
683b273b 498 par[2] = 68.35;
452a64c6 499 /*par[0] = 66.55;
500 par[1] = 11.5;
501 par[2] = 64.8;*/
502 gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3);
503
6197ac87 504 // Air 2 (cutting the lower part of the box)
505
506 par[0] = 1.25;
507 par[1] = 3; //Original Settings
508 par[2] = 70.86;
509 gMC->Gsvolu("AIR2", "BOX ", idtmed[1000], par, 3);
510
511 // Air 3 (cutting the lower part of the box)
512
513 par[0] = 66.3;
514 par[1] = 3; //Original Settings
515 par[2] = 1.2505;
516 gMC->Gsvolu("AIR3", "BOX ", idtmed[1000], par, 3);
517
452a64c6 518 // Honeycomb
e293afae 519 par[0] = 66.3;
452a64c6 520 par[1] = .188; //Original Settings
e293afae 521 par[2] = 68.35;
452a64c6 522 /*par[0] = 66.55;
523 par[1] = .188;
524 par[2] = 63.1;*/
525 gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
526
527 // Aluminium sheet
e293afae 528 par[0] = 66.3;
452a64c6 529 par[1] = .025; //Original Settings
e293afae 530 par[2] = 68.35;
452a64c6 531 /*par[0] = 66.5;
532 par[1] = .025;
533 par[2] = 63.1;*/
534 gMC->Gsvolu("ALUM", "BOX ", idtmed[1009], par, 3);
535
536 // Quartz
537 par[0] = geometry->GetQuartzWidth()/2;
538 par[1] = geometry->GetQuartzThickness()/2;
539 par[2] = geometry->GetQuartzLength()/2;
540 /*par[0] = 63.1;
541 par[1] = .25; //Original Settings
542 par[2] = 65.5;*/
543 /*par[0] = geometry->GetQuartzWidth()/2;
544 par[1] = geometry->GetQuartzThickness()/2;
545 par[2] = geometry->GetQuartzLength()/2;*/
546 //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]);
547 gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
548
549 // Spacers (cylinders)
550 par[0] = 0.;
551 par[1] = .5;
552 par[2] = geometry->GetFreonThickness()/2;
553 gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
554
683b273b 555 // Feet (freon slabs supports)
556
557 par[0] = .7;
558 par[1] = .3;
559 par[2] = 1.9;
560 gMC->Gsvolu("FOOT", "BOX", idtmed[1009], par, 3);
561
452a64c6 562 // Opaque quartz
e293afae 563 par[0] = geometry->GetQuartzWidth()/2;
564 par[1] = .2;
565 par[2] = geometry->GetQuartzLength()/2;
566 /*par[0] = 61.95;
452a64c6 567 par[1] = .2; //Original Settings
e293afae 568 par[2] = 66.5;*/
452a64c6 569 /*par[0] = 66.5;
570 par[1] = .2;
571 par[2] = 61.95;*/
572 gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
573
574 // Frame of opaque quartz
53a60579 575 par[0] = geometry->GetOuterFreonWidth()/2;
576 //+ oqua_thickness;
452a64c6 577 par[1] = geometry->GetFreonThickness()/2;
53a60579 578 par[2] = geometry->GetOuterFreonLength()/2;
579 //+ oqua_thickness;
452a64c6 580 /*par[0] = 20.65;
581 par[1] = .5; //Original Settings
582 par[2] = 66.5;*/
583 /*par[0] = 66.5;
584 par[1] = .5;
585 par[2] = 20.65;*/
586 gMC->Gsvolu("OQF1", "BOX ", idtmed[1007], par, 3);
587
e293afae 588 par[0] = geometry->GetInnerFreonWidth()/2;
452a64c6 589 par[1] = geometry->GetFreonThickness()/2;
e293afae 590 par[2] = geometry->GetInnerFreonLength()/2;
452a64c6 591 gMC->Gsvolu("OQF2", "BOX ", idtmed[1007], par, 3);
592
593 // Little bar of opaque quartz
e293afae 594 //par[0] = .275;
595 //par[1] = geometry->GetQuartzThickness()/2;
53a60579 596 //par[2] = geometry->GetInnerFreonLength()/2 - 2.4;
e293afae 597 //par[2] = geometry->GetInnerFreonLength()/2;
53a60579 598 //+ oqua_thickness;
452a64c6 599 /*par[0] = .275;
600 par[1] = .25; //Original Settings
601 par[2] = 63.1;*/
602 /*par[0] = 63.1;
603 par[1] = .25;
604 par[2] = .275;*/
e293afae 605 //gMC->Gsvolu("BARR", "BOX ", idtmed[1007], par, 3);
452a64c6 606
607 // Freon
e293afae 608 par[0] = geometry->GetOuterFreonWidth()/2 - oqua_thickness;
452a64c6 609 par[1] = geometry->GetFreonThickness()/2;
e293afae 610 par[2] = geometry->GetOuterFreonLength()/2 - 2*oqua_thickness;
452a64c6 611 /*par[0] = 20.15;
612 par[1] = .5; //Original Settings
613 par[2] = 65.5;*/
614 /*par[0] = 65.5;
615 par[1] = .5;
616 par[2] = 20.15;*/
617 gMC->Gsvolu("FRE1", "BOX ", idtmed[1003], par, 3);
618
e293afae 619 par[0] = geometry->GetInnerFreonWidth()/2 - oqua_thickness;
452a64c6 620 par[1] = geometry->GetFreonThickness()/2;
e293afae 621 par[2] = geometry->GetInnerFreonLength()/2 - 2*oqua_thickness;
452a64c6 622 gMC->Gsvolu("FRE2", "BOX ", idtmed[1003], par, 3);
623
624 // Methane
e293afae 625 //par[0] = 64.8;
683b273b 626 par[0] = csi_width/2;
452a64c6 627 par[1] = geometry->GetGapThickness()/2;
628 //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 629 //par[2] = 64.8;
683b273b 630 par[2] = csi_length/2;
452a64c6 631 gMC->Gsvolu("META", "BOX ", idtmed[1004], par, 3);
632
633 // Methane gap
e293afae 634 //par[0] = 64.8;
683b273b 635 par[0] = csi_width/2;
452a64c6 636 par[1] = geometry->GetProximityGapThickness()/2;
637 //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 638 //par[2] = 64.8;
683b273b 639 par[2] = csi_length/2;
452a64c6 640 gMC->Gsvolu("GAP ", "BOX ", idtmed[1008], par, 3);
641
642 // CsI photocathode
e293afae 643 //par[0] = 64.8;
683b273b 644 par[0] = csi_width/2;
452a64c6 645 par[1] = .25;
e293afae 646 //par[2] = 64.8;
683b273b 647 par[2] = csi_length/2;
452a64c6 648 gMC->Gsvolu("CSI ", "BOX ", idtmed[1005], par, 3);
649
650 // Anode grid
651 par[0] = 0.;
652 par[1] = .001;
653 par[2] = 20.;
654 gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
683b273b 655
6197ac87 656 // Wire supports
657 // Bar of metal
658
659 par[0] = csi_width/2;
660 par[1] = 1.05;
661 par[2] = 1.05;
662 gMC->Gsvolu("WSMe", "BOX ", idtmed[1009], par, 3);
663
664 // Ceramic pick up (base)
665
666 par[0] = csi_width/2;
667 par[1] = .25;
668 par[2] = 1.05;
669 gMC->Gsvolu("WSG1", "BOX ", idtmed[1010], par, 3);
670
671 // Ceramic pick up (head)
672
673 par[0] = csi_width/2;
674 par[1] = .1;
675 par[2] = .1;
676 gMC->Gsvolu("WSG2", "BOX ", idtmed[1010], par, 3);
677
683b273b 678 // Aluminium supports for methane and CsI
679 // Short bar
680
681 par[0] = csi_width/2;
682 par[1] = geometry->GetGapThickness()/2 + .25;
683 par[2] = (68.35 - csi_length/2)/2;
684 gMC->Gsvolu("SMSH", "BOX", idtmed[1009], par, 3);
452a64c6 685
683b273b 686 // Long bar
687
688 par[0] = (66.3 - csi_width/2)/2;
689 par[1] = geometry->GetGapThickness()/2 + .25;
690 par[2] = csi_length/2 + 68.35 - csi_length/2;
691 gMC->Gsvolu("SMLG", "BOX", idtmed[1009], par, 3);
692
693 // Aluminium supports for freon
694 // Short bar
695
696 par[0] = geometry->GetQuartzWidth()/2;
697 par[1] = .3;
698 par[2] = (68.35 - geometry->GetQuartzLength()/2)/2;
699 gMC->Gsvolu("SFSH", "BOX", idtmed[1009], par, 3);
700
701 // Long bar
702
703 par[0] = (66.3 - geometry->GetQuartzWidth()/2)/2;
704 par[1] = .3;
705 par[2] = geometry->GetQuartzLength()/2 + 68.35 - geometry->GetQuartzLength()/2;
706 gMC->Gsvolu("SFLG", "BOX", idtmed[1009], par, 3);
707
6197ac87 708 // PCB backplane
683b273b 709
6197ac87 710 par[0] = csi_width/2;
711 par[1] = .25;
712 par[2] = csi_length/4 -.5025;
713 gMC->Gsvolu("PCB ", "BOX", idtmed[1011], par, 3);
683b273b 714
6197ac87 715
716 // Backplane supports
717
718 // Aluminium slab
719
720 par[0] = 33.15;
721 par[1] = 2;
722 par[2] = 21.65;
723 gMC->Gsvolu("BACK", "BOX", idtmed[1009], par, 3);
724
725 // Big hole
726
727 par[0] = 9.05;
728 par[1] = 2;
729 par[2] = 4.4625;
730 gMC->Gsvolu("BKHL", "BOX", idtmed[1000], par, 3);
731
732 // Small hole
733
734 par[0] = 5.7;
735 par[1] = 2;
736 par[2] = 4.4625;
737 gMC->Gsvolu("BKHS", "BOX", idtmed[1000], par, 3);
738
739 // Place holes inside backplane support
740
741 gMC->Gspos("BKHS", 1, "BACK", .8 + 5.7,0., .6 + 4.4625, 0, "ONLY");
742 gMC->Gspos("BKHS", 2, "BACK", -.8 - 5.7,0., .6 + 4.4625, 0, "ONLY");
743 gMC->Gspos("BKHS", 3, "BACK", .8 + 5.7,0., -.6 - 4.4625, 0, "ONLY");
744 gMC->Gspos("BKHS", 4, "BACK", -.8 - 5.7,0., -.6 - 4.4625, 0, "ONLY");
745 gMC->Gspos("BKHS", 5, "BACK", .8 + 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
746 gMC->Gspos("BKHS", 6, "BACK", -.8 - 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
747 gMC->Gspos("BKHS", 7, "BACK", .8 + 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
748 gMC->Gspos("BKHS", 8, "BACK", -.8 - 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
749 gMC->Gspos("BKHL", 1, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., .6 + 4.4625, 0, "ONLY");
750 gMC->Gspos("BKHL", 2, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 4.4625, 0, "ONLY");
751 gMC->Gspos("BKHL", 3, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 4.4625, 0, "ONLY");
752 gMC->Gspos("BKHL", 4, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 4.4625, 0, "ONLY");
753 gMC->Gspos("BKHL", 5, "BACK", .8 + 11.4+ 1.6 + 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
754 gMC->Gspos("BKHL", 6, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
755 gMC->Gspos("BKHL", 7, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
756 gMC->Gspos("BKHL", 8, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
757
758
759
452a64c6 760 // --- Places the detectors defined with GSVOLU
761 // Place material inside RICH
6197ac87 762 gMC->Gspos("SRIC", 1, "RICH", 0.,0., 0., 0, "ONLY");
763 gMC->Gspos("AIR2", 1, "RICH", 66.3 + 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.5, 0., 0, "ONLY");
764 gMC->Gspos("AIR2", 2, "RICH", -66.3 - 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.5, 0., 0, "ONLY");
765 gMC->Gspos("AIR3", 1, "RICH", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.5, -68.35 - 1.25, 0, "ONLY");
766 gMC->Gspos("AIR3", 2, "RICH", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.5, 68.35 + 1.25, 0, "ONLY");
767
683b273b 768
769 gMC->Gspos("ALUM", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.025, 0., 0, "ONLY");
770 gMC->Gspos("HONE", 1, "SRIC", 0., 1.276- geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .188, 0., 0, "ONLY");
771 gMC->Gspos("ALUM", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .025, 0., 0, "ONLY");
772 gMC->Gspos("FOOT", 1, "SRIC", 64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
773 gMC->Gspos("FOOT", 2, "SRIC", 21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3 , 36.9, 0, "ONLY");
774 gMC->Gspos("FOOT", 3, "SRIC", -21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
775 gMC->Gspos("FOOT", 4, "SRIC", -64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
776 gMC->Gspos("FOOT", 5, "SRIC", 64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
777 gMC->Gspos("FOOT", 6, "SRIC", 21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
778 gMC->Gspos("FOOT", 7, "SRIC", -21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
779 gMC->Gspos("FOOT", 8, "SRIC", -64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
452a64c6 780 gMC->Gspos("OQUA", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .2, 0., 0, "ONLY");
781
683b273b 782 // Supports placing
783
784 // Methane supports
785 gMC->Gspos("SMLG", 1, "SRIC", csi_width/2 + (66.3 - csi_width/2)/2, 1.276 + .25, 0., 0, "ONLY");
786 gMC->Gspos("SMLG", 2, "SRIC", - csi_width/2 - (66.3 - csi_width/2)/2, 1.276 + .25, 0., 0, "ONLY");
787 gMC->Gspos("SMSH", 1, "SRIC", 0., 1.276 + .25, csi_length/2 + (68.35 - csi_length/2)/2, 0, "ONLY");
788 gMC->Gspos("SMSH", 2, "SRIC", 0., 1.276 + .25, - csi_length/2 - (68.35 - csi_length/2)/2, 0, "ONLY");
789
790 //Freon supports
791
792 Float_t supp_y = 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness() - .2 + .3; //y position of freon supports
793
794 gMC->Gspos("SFLG", 1, "SRIC", geometry->GetQuartzWidth()/2 + (66.3 - geometry->GetQuartzWidth()/2)/2, supp_y, 0., 0, "ONLY");
795 gMC->Gspos("SFLG", 2, "SRIC", - geometry->GetQuartzWidth()/2 - (66.3 - geometry->GetQuartzWidth()/2)/2, supp_y, 0., 0, "ONLY");
796 gMC->Gspos("SFSH", 1, "SRIC", 0., supp_y, geometry->GetQuartzLength()/2 + (68.35 - geometry->GetQuartzLength()/2)/2, 0, "ONLY");
797 gMC->Gspos("SFSH", 2, "SRIC", 0., supp_y, - geometry->GetQuartzLength()/2 - (68.35 - geometry->GetQuartzLength()/2)/2, 0, "ONLY");
798
452a64c6 799 AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
800
e293afae 801 //Int_t nspacers = (Int_t)(TMath::Abs(geometry->GetInnerFreonLength()/14.4));
802 Int_t nspacers = 9;
452a64c6 803 //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);
804
805 //printf("Nspacers: %d", nspacers);
806
807 //for (i = 1; i <= 9; ++i) {
808 //zs = (5 - i) * 14.4; //Original settings
809 for (i = 0; i < nspacers; i++) {
810 zs = (TMath::Abs(nspacers/2) - i) * 14.4;
811 gMC->Gspos("SPAC", i, "FRE1", 6.7, 0., zs, idrotm[1019], "ONLY"); //Original settings
812 //gMC->Gspos("SPAC", i, "FRE1", zs, 0., 6.7, idrotm[1019], "ONLY");
813 }
814 //for (i = 10; i <= 18; ++i) {
815 //zs = (14 - i) * 14.4; //Original settings
816 for (i = nspacers; i < nspacers*2; ++i) {
817 zs = (nspacers + TMath::Abs(nspacers/2) - i) * 14.4;
818 gMC->Gspos("SPAC", i, "FRE1", -6.7, 0., zs, idrotm[1019], "ONLY"); //Original settings
819 //gMC->Gspos("SPAC", i, "FRE1", zs, 0., -6.7, idrotm[1019], "ONLY");
820 }
821
822 //for (i = 1; i <= 9; ++i) {
823 //zs = (5 - i) * 14.4; //Original settings
824 for (i = 0; i < nspacers; i++) {
825 zs = (TMath::Abs(nspacers/2) - i) * 14.4;
826 gMC->Gspos("SPAC", i, "FRE2", 6.7, 0., zs, idrotm[1019], "ONLY"); //Original settings
827 //gMC->Gspos("SPAC", i, "FRE2", zs, 0., 6.7, idrotm[1019], "ONLY");
828 }
829 //for (i = 10; i <= 18; ++i) {
830 //zs = (5 - i) * 14.4; //Original settings
831 for (i = nspacers; i < nspacers*2; ++i) {
832 zs = (nspacers + TMath::Abs(nspacers/2) - i) * 14.4;
833 gMC->Gspos("SPAC", i, "FRE2", -6.7, 0., zs, idrotm[1019], "ONLY"); //Original settings
834 //gMC->Gspos("SPAC", i, "FRE2", zs, 0., -6.7, idrotm[1019], "ONLY");
835 }
836
837 /*gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY");
838 gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY");
839 gMC->Gspos("OQF1", 1, "SRIC", 31.3, -4.724, 41.3, 0, "ONLY");
840 gMC->Gspos("OQF2", 2, "SRIC", 0., -4.724, 0., 0, "ONLY");
841 gMC->Gspos("OQF1", 3, "SRIC", -31.3, -4.724, -41.3, 0, "ONLY");
842 gMC->Gspos("BARR", 1, "QUAR", -21.65, 0., 0., 0, "ONLY"); //Original settings
843 gMC->Gspos("BARR", 2, "QUAR", 21.65, 0., 0., 0, "ONLY"); //Original settings
844 gMC->Gspos("QUAR", 1, "SRIC", 0., -3.974, 0., 0, "ONLY");
845 gMC->Gspos("GAP ", 1, "META", 0., 4.8, 0., 0, "ONLY");
846 gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
847 gMC->Gspos("CSI ", 1, "SRIC", 0., 6.526, 0., 0, "ONLY");*/
848
849
850 gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY");
851 gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY");
683b273b 852 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)
452a64c6 853 gMC->Gspos("OQF2", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY"); //Original settings
683b273b 854 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 855 //gMC->Gspos("BARR", 1, "QUAR", - geometry->GetInnerFreonWidth()/2 - oqua_thickness, 0., 0., 0, "ONLY"); //Original settings (-21.65)
856 //gMC->Gspos("BARR", 2, "QUAR", geometry->GetInnerFreonWidth()/2 + oqua_thickness, 0., 0., 0, "ONLY"); //Original settings (21.65)
452a64c6 857 gMC->Gspos("QUAR", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness()/2, 0., 0, "ONLY");
858 gMC->Gspos("GAP ", 1, "META", 0., geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 - 0.0001, 0., 0, "ONLY");
859 gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
860 gMC->Gspos("CSI ", 1, "SRIC", 0., 1.276 + geometry->GetGapThickness()/2 + .25, 0., 0, "ONLY");
6197ac87 861
862 // Wire support placing
863
864 gMC->Gspos("WSG2", 1, "GAP ", 0., geometry->GetProximityGapThickness()/2 - .1, 0., 0, "ONLY");
865 gMC->Gspos("WSG1", 1, "CSI ", 0., 0., 0., 0, "ONLY");
866 gMC->Gspos("WSMe", 1, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, 0., 0, "ONLY");
867
868 // Backplane placing
869
870 gMC->Gspos("BACK", 1, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 43.3, 0, "ONLY");
871 gMC->Gspos("BACK", 2, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2 , 43.3, 0, "ONLY");
872 gMC->Gspos("BACK", 3, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
873 gMC->Gspos("BACK", 4, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
874 gMC->Gspos("BACK", 5, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
875 gMC->Gspos("BACK", 6, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
876
877 // PCB placing
878
879 gMC->Gspos("PCB ", 1, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, csi_width/4 + .5025 + 1.2, 0, "ONLY");
880 gMC->Gspos("PCB ", 2, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, -csi_width/4 - .5025 - 1.2, 0, "ONLY");
881
882
452a64c6 883
884 //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);
885
886 // Place RICH inside ALICE apparatus
887
888 AliMatrix(idrotm[1000], 90., 0., 70.69, 90., 19.31, -90.);
889 AliMatrix(idrotm[1001], 90., -20., 90., 70., 0., 0.);
890 AliMatrix(idrotm[1002], 90., 0., 90., 90., 0., 0.);
891 AliMatrix(idrotm[1003], 90., 20., 90., 110., 0., 0.);
892 AliMatrix(idrotm[1004], 90., 340., 108.2, 70., 18.2, 70.);
893 AliMatrix(idrotm[1005], 90., 0., 109.31, 90., 19.31, 90.);
894 AliMatrix(idrotm[1006], 90., 20., 108.2, 110., 18.2, 110.);
895
896 gMC->Gspos("RICH", 1, "ALIC", 0., 471.9, 165.26, idrotm[1000], "ONLY");
897 gMC->Gspos("RICH", 2, "ALIC", 171., 470., 0., idrotm[1001], "ONLY");
898 gMC->Gspos("RICH", 3, "ALIC", 0., 500., 0., idrotm[1002], "ONLY");
899 gMC->Gspos("RICH", 4, "ALIC", -171., 470., 0., idrotm[1003], "ONLY");
900 gMC->Gspos("RICH", 5, "ALIC", 161.4, 443.4, -165.3, idrotm[1004], "ONLY");
901 gMC->Gspos("RICH", 6, "ALIC", 0., 471.9, -165.3, idrotm[1005], "ONLY");
902 gMC->Gspos("RICH", 7, "ALIC", -161.4, 443.4, -165.3, idrotm[1006], "ONLY");
903
904}
905
906
907//___________________________________________
908void AliRICH::CreateMaterials()
909{
910 //
911 // *** DEFINITION OF AVAILABLE RICH MATERIALS ***
912 // ORIGIN : NICK VAN EIJNDHOVEN
913 // Modified by: N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
914 // R.A. Fini (INFN - BARI, Rosanna.Fini@ba.infn.it)
915 // R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
916 //
917 Int_t isxfld = gAlice->Field()->Integ();
918 Float_t sxmgmx = gAlice->Field()->Max();
919 Int_t i;
920
921 /************************************Antonnelo's Values (14-vectors)*****************************************/
922 /*
923 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,
924 6.7e-9,6.88e-9,7.08e-9,7.3e-9,7.51e-9,7.74e-9,8e-9 };
925 Float_t rIndexQuarz[14] = { 1.528309,1.533333,
926 1.538243,1.544223,1.550568,1.55777,
927 1.565463,1.574765,1.584831,1.597027,
928 1.611858,1.6277,1.6472,1.6724 };
929 Float_t rIndexOpaqueQuarz[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
930 Float_t rIndexMethane[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
931 Float_t rIndexGrid[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
932 Float_t abscoFreon[14] = { 179.0987,179.0987,
933 179.0987,179.0987,179.0987,142.92,56.65,13.95,10.43,7.07,2.03,.5773,.33496,0. };
934 //Float_t abscoFreon[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
935 // 1e-5,1e-5,1e-5,1e-5,1e-5 };
936 Float_t abscoQuarz[14] = { 64.035,39.98,35.665,31.262,27.527,22.815,21.04,17.52,
937 14.177,9.282,4.0925,1.149,.3627,.10857 };
938 Float_t abscoOpaqueQuarz[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
939 1e-5,1e-5,1e-5,1e-5,1e-5 };
940 Float_t abscoCsI[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
941 1e-4,1e-4,1e-4,1e-4 };
942 Float_t abscoMethane[14] = { 1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,
943 1e6,1e6,1e6 };
944 Float_t abscoGrid[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
945 1e-4,1e-4,1e-4,1e-4 };
946 Float_t efficAll[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
947 Float_t efficCsI[14] = { 6e-4,.005,.0075,.01125,.045,.117,.135,.16575,
948 .17425,.1785,.1836,.1904,.1938,.221 };
949 Float_t efficGrid[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
950 */
951
952
953 /**********************************End of Antonnelo's Values**********************************/
954
955 /**********************************Values from rich_media.f (31-vectors)**********************************/
956
957
958 //Photons energy intervals
959 Float_t ppckov[26];
960 for (i=0;i<26;i++)
961 {
962 ppckov[i] = (Float_t(i)*0.1+5.5)*1e-9;
963 //printf ("Energy intervals: %e\n",ppckov[i]);
964 }
965
966
967 //Refraction index for quarz
968 Float_t rIndexQuarz[26];
969 Float_t e1= 10.666;
970 Float_t e2= 18.125;
971 Float_t f1= 46.411;
972 Float_t f2= 228.71;
973 for (i=0;i<26;i++)
974 {
975 Float_t ene=ppckov[i]*1e9;
976 Float_t a=f1/(e1*e1 - ene*ene);
977 Float_t b=f2/(e2*e2 - ene*ene);
978 rIndexQuarz[i] = TMath::Sqrt(1. + a + b );
979 //printf ("rIndexQuarz: %e\n",rIndexQuarz[i]);
980 }
981
982 //Refraction index for opaque quarz, methane and grid
983 Float_t rIndexOpaqueQuarz[26];
984 Float_t rIndexMethane[26];
985 Float_t rIndexGrid[26];
986 for (i=0;i<26;i++)
987 {
988 rIndexOpaqueQuarz[i]=1;
989 rIndexMethane[i]=1.000444;
990 rIndexGrid[i]=1;
991 //printf ("rIndexOpaqueQuarz , etc: %e, %e, %e\n",rIndexOpaqueQuarz[i], rIndexMethane[i], rIndexGrid[i]=1);
992 }
993
994 //Absorption index for freon
995 Float_t abscoFreon[26] = {179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987,
996 179.0987, 142.9206, 56.64957, 25.58622, 13.95293, 12.03905, 10.42953, 8.804196,
997 7.069031, 4.461292, 2.028366, 1.293013, .577267, .40746, .334964, 0., 0., 0.};
998
999 //Absorption index for quarz
1000 /*Float_t Qzt [21] = {.0,.0,.005,.04,.35,.647,.769,.808,.829,.844,.853,.858,.869,.887,.903,.902,.902,
1001 .906,.907,.907,.907};
1002 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,
1003 215.0,220.0,225.0,230.0,235.0,240.0,245.0,250.0};
1004 Float_t abscoQuarz[31];
1005 for (Int_t i=0;i<31;i++)
1006 {
1007 Float_t Xlam = 1237.79 / (ppckov[i]*1e9);
1008 if (Xlam <= 160) abscoQuarz[i] = 0;
1009 if (Xlam > 250) abscoQuarz[i] = 1;
1010 else
1011 {
1012 for (Int_t j=0;j<21;j++)
1013 {
1014 //printf ("Passed\n");
1015 if (Xlam > Wavl2[j] && Xlam < Wavl2[j+1])
1016 {
1017 Float_t Dabs = (Qzt[j+1] - Qzt[j])/(Wavl2[j+1] - Wavl2[j]);
1018 Float_t Abso = Qzt[j] + Dabs*(Xlam - Wavl2[j]);
1019 abscoQuarz[i] = -5.0/(TMath::Log(Abso));
1020 }
1021 }
1022 }
1023 printf ("abscoQuarz: %e abscoFreon: %e for energy: %e\n",abscoQuarz[i],abscoFreon[i],ppckov[i]);
1024 }*/
1025
1026 /*Float_t abscoQuarz[31] = {49.64211, 48.41296, 47.46989, 46.50492, 45.13682, 44.47883, 43.1929 , 41.30922, 40.5943 ,
1027 39.82956, 38.98623, 38.6247 , 38.43448, 37.41084, 36.22575, 33.74852, 30.73901, 24.25086,
1028 17.94531, 11.88753, 5.99128, 3.83503, 2.36661, 1.53155, 1.30582, 1.08574, .8779708,
1029 .675275, 0., 0., 0.};
1030
1031 for (Int_t i=0;i<31;i++)
1032 {
1033 abscoQuarz[i] = abscoQuarz[i]/10;
1034 }*/
1035
1036 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,
1037 19.266, 17.525, 15.878, 14.177, 11.719, 9.282, 6.62, 4.0925, 2.601, 1.149, .667, .3627,
1038 .192, .1497, .10857};
1039
1040 //Absorption index for methane
1041 Float_t abscoMethane[26];
1042 for (i=0;i<26;i++)
1043 {
1044 abscoMethane[i]=AbsoCH4(ppckov[i]*1e9);
1045 //printf("abscoMethane: %e for energy: %e\n", abscoMethane[i],ppckov[i]*1e9);
1046 }
1047
1048 //Absorption index for opaque quarz, csi and grid, efficiency for all and grid
1049 Float_t abscoOpaqueQuarz[26];
1050 Float_t abscoCsI[26];
1051 Float_t abscoGrid[26];
1052 Float_t efficAll[26];
1053 Float_t efficGrid[26];
1054 for (i=0;i<26;i++)
1055 {
1056 abscoOpaqueQuarz[i]=1e-5;
1057 abscoCsI[i]=1e-4;
1058 abscoGrid[i]=1e-4;
1059 efficAll[i]=1;
1060 efficGrid[i]=1;
1061 //printf ("All must be 1: %e, %e, %e, %e, %e\n",abscoOpaqueQuarz[i],abscoCsI[i],abscoGrid[i],efficAll[i],efficGrid[i]);
1062 }
1063
1064 //Efficiency for csi
1065
1066 Float_t efficCsI[26] = {0.000199999995, 0.000600000028, 0.000699999975, 0.00499999989, 0.00749999983, 0.010125,
1067 0.0242999997, 0.0405000001, 0.0688500032, 0.105299994, 0.121500008, 0.141749993, 0.157949999,
1068 0.162, 0.166050002, 0.167669997, 0.174299985, 0.176789999, 0.179279998, 0.182599992, 0.18592,
1069 0.187579989, 0.189239994, 0.190899998, 0.207499996, 0.215799987};
1070
1071
1072
1073 //FRESNEL LOSS CORRECTION FOR PERPENDICULAR INCIDENCE AND
1074 //UNPOLARIZED PHOTONS
1075
1076 for (i=0;i<26;i++)
1077 {
1078 efficCsI[i] = efficCsI[i]/(1.-Fresnel(ppckov[i]*1e9,1.,0));
1079 //printf ("Fresnel result: %e for energy: %e\n",Fresnel(ppckov[i]*1e9,1.,0),ppckov[i]*1e9);
1080 }
1081
1082 /*******************************************End of rich_media.f***************************************/
1083
1084
1085
1086
1087
1088
1089 Float_t afre[2], agri, amet[2], aqua[2], ahon, zfre[2], zgri, zhon,
1090 zmet[2], zqua[2];
1091 Int_t nlmatfre;
1092 Float_t densquao;
1093 Int_t nlmatmet, nlmatqua;
1094 Float_t wmatquao[2], rIndexFreon[26];
1095 Float_t aquao[2], epsil, stmin, zquao[2];
1096 Int_t nlmatquao;
1097 Float_t radlal, densal, tmaxfd, deemax, stemax;
1098 Float_t aal, zal, radlgri, densfre, radlhon, densgri, denshon,densqua, densmet, wmatfre[2], wmatmet[2], wmatqua[2];
1099
1100 Int_t *idtmed = fIdtmed->GetArray()-999;
1101
452a64c6 1102 // --- Photon energy (GeV)
1103 // --- Refraction indexes
1104 for (i = 0; i < 26; ++i) {
1105 rIndexFreon[i] = ppckov[i] * .0172 * 1e9 + 1.177;
1106 //rIndexFreon[i] = 1;
1107 //printf ("rIndexFreon: %e \n efficCsI: %e for energy: %e\n",rIndexFreon[i], efficCsI[i], ppckov[i]);
1108 }
1109
1110 // --- Detection efficiencies (quantum efficiency for CsI)
1111 // --- Define parameters for honeycomb.
1112 // Used carbon of equivalent rad. lenght
1113
1114 ahon = 12.01;
1115 zhon = 6.;
9dda3582 1116 denshon = 0.1;
452a64c6 1117 radlhon = 18.8;
1118
1119 // --- Parameters to include in GSMIXT, relative to Quarz (SiO2)
1120
1121 aqua[0] = 28.09;
1122 aqua[1] = 16.;
1123 zqua[0] = 14.;
1124 zqua[1] = 8.;
1125 densqua = 2.64;
1126 nlmatqua = -2;
1127 wmatqua[0] = 1.;
1128 wmatqua[1] = 2.;
1129
1130 // --- Parameters to include in GSMIXT, relative to opaque Quarz (SiO2)
1131
1132 aquao[0] = 28.09;
1133 aquao[1] = 16.;
1134 zquao[0] = 14.;
1135 zquao[1] = 8.;
1136 densquao = 2.64;
1137 nlmatquao = -2;
1138 wmatquao[0] = 1.;
1139 wmatquao[1] = 2.;
1140
1141 // --- Parameters to include in GSMIXT, relative to Freon (C6F14)
1142
1143 afre[0] = 12.;
1144 afre[1] = 19.;
1145 zfre[0] = 6.;
1146 zfre[1] = 9.;
1147 densfre = 1.7;
1148 nlmatfre = -2;
1149 wmatfre[0] = 6.;
1150 wmatfre[1] = 14.;
1151
1152 // --- Parameters to include in GSMIXT, relative to methane (CH4)
1153
1154 amet[0] = 12.01;
1155 amet[1] = 1.;
1156 zmet[0] = 6.;
1157 zmet[1] = 1.;
1158 densmet = 7.17e-4;
1159 nlmatmet = -2;
1160 wmatmet[0] = 1.;
1161 wmatmet[1] = 4.;
1162
1163 // --- Parameters to include in GSMIXT, relative to anode grid (Cu)
1164
1165 agri = 63.54;
1166 zgri = 29.;
1167 densgri = 8.96;
1168 radlgri = 1.43;
1169
1170 // --- Parameters to include in GSMATE related to aluminium sheet
1171
1172 aal = 26.98;
1173 zal = 13.;
1174 densal = 2.7;
1175 radlal = 8.9;
6197ac87 1176
1177 // --- Glass parameters
1178
1179 Float_t aglass[5]={12.01, 28.09, 16., 10.8, 23.};
1180 Float_t zglass[5]={ 6., 14., 8., 5., 11.};
1181 Float_t wglass[5]={ 0.5, 0.105, 0.355, 0.03, 0.01};
1182 Float_t dglass=1.74;
1183
452a64c6 1184
1185 AliMaterial(1, "Air $", 14.61, 7.3, .001205, 30420., 67500);
1186 AliMaterial(6, "HON", ahon, zhon, denshon, radlhon, 0);
1187 AliMaterial(16, "CSI", ahon, zhon, denshon, radlhon, 0);
1188 AliMixture(20, "QUA", aqua, zqua, densqua, nlmatqua, wmatqua);
1189 AliMixture(21, "QUAO", aquao, zquao, densquao, nlmatquao, wmatquao);
1190 AliMixture(30, "FRE", afre, zfre, densfre, nlmatfre, wmatfre);
1191 AliMixture(40, "MET", amet, zmet, densmet, nlmatmet, wmatmet);
1192 AliMixture(41, "METG", amet, zmet, densmet, nlmatmet, wmatmet);
1193 AliMaterial(11, "GRI", agri, zgri, densgri, radlgri, 0);
1194 AliMaterial(50, "ALUM", aal, zal, densal, radlal, 0);
6197ac87 1195 AliMixture(32, "GLASS",aglass, zglass, dglass, 5, wglass);
1196 AliMaterial(31, "COPPER$", 63.54, 29., 8.96, 1.4, 0.);
452a64c6 1197
1198 tmaxfd = -10.;
1199 stemax = -.1;
1200 deemax = -.2;
1201 epsil = .001;
1202 stmin = -.001;
1203
1204 AliMedium(1, "DEFAULT MEDIUM AIR$", 1, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1205 AliMedium(2, "HONEYCOMB$", 6, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1206 AliMedium(3, "QUARZO$", 20, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1207 AliMedium(4, "FREON$", 30, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1208 AliMedium(5, "METANO$", 40, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1209 AliMedium(6, "CSI$", 16, 1, isxfld, sxmgmx,tmaxfd, stemax, deemax, epsil, stmin);
1210 AliMedium(7, "GRIGLIA$", 11, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1211 AliMedium(8, "QUARZOO$", 21, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1212 AliMedium(9, "GAP$", 41, 1, isxfld, sxmgmx,tmaxfd, .1, -deemax, epsil, -stmin);
1213 AliMedium(10, "ALUMINUM$", 50, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
6197ac87 1214 AliMedium(11, "GLASS", 32, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1215 AliMedium(12, "PCB_COPPER", 31, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
452a64c6 1216
1217
5cfcdf54 1218 gMC->SetCerenkov(idtmed[1000], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1219 gMC->SetCerenkov(idtmed[1001], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1220 gMC->SetCerenkov(idtmed[1002], 26, ppckov, abscoQuarz, efficAll,rIndexQuarz);
1221 gMC->SetCerenkov(idtmed[1003], 26, ppckov, abscoFreon, efficAll,rIndexFreon);
1222 gMC->SetCerenkov(idtmed[1004], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1223 gMC->SetCerenkov(idtmed[1005], 26, ppckov, abscoCsI, efficCsI, rIndexMethane);
1224 gMC->SetCerenkov(idtmed[1006], 26, ppckov, abscoGrid, efficGrid, rIndexGrid);
1225 gMC->SetCerenkov(idtmed[1007], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz);
1226 gMC->SetCerenkov(idtmed[1008], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1227 gMC->SetCerenkov(idtmed[1009], 26, ppckov, abscoGrid, efficGrid, rIndexGrid);
1228 gMC->SetCerenkov(idtmed[1010], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz);
452a64c6 1229}
1230
1231//___________________________________________
1232
1233Float_t AliRICH::Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
1234{
1235
1236 //ENE(EV), PDOTI=COS(INC.ANG.), PDOTR=COS(POL.PLANE ROT.ANG.)
1237
1238 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,
1239 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,
1240 7.8,7.9,8.0,8.1,8.2,8.3,8.4,8.5};
1241
1242
1243 Float_t csin[36] = {2.14,2.21,2.33,2.48,2.76,2.97,2.99,2.59,2.81,3.05,
1244 2.86,2.53,2.55,2.66,2.79,2.96,3.18,3.05,2.84,2.81,2.38,2.11,
1245 2.01,2.13,2.39,2.73,3.08,3.15,2.95,2.73,2.56,2.41,2.12,1.95,
1246 1.72,1.53};
1247
1248 Float_t csik[36] = {0.,0.,0.,0.,0.,0.196,0.408,0.208,0.118,0.49,0.784,0.543,
1249 0.424,0.404,0.371,0.514,0.922,1.102,1.139,1.376,1.461,1.253,0.878,
1250 0.69,0.612,0.649,0.824,1.347,1.571,1.678,1.763,1.857,1.824,1.824,
1251 1.714,1.498};
1252 Float_t xe=ene;
1253 Int_t j=Int_t(xe*10)-49;
1254 Float_t cn=csin[j]+((csin[j+1]-csin[j])/0.1)*(xe-en[j]);
1255 Float_t ck=csik[j]+((csik[j+1]-csik[j])/0.1)*(xe-en[j]);
1256
1257 //FORMULAE FROM HANDBOOK OF OPTICS, 33.23 OR
1258 //W.R. HUNTER, J.O.S.A. 54 (1964),15 , J.O.S.A. 55(1965),1197
1259
1260 Float_t sinin=TMath::Sqrt(1-pdoti*pdoti);
1261 Float_t tanin=sinin/pdoti;
1262
1263 Float_t c1=cn*cn-ck*ck-sinin*sinin;
1264 Float_t c2=4*cn*cn*ck*ck;
1265 Float_t aO=TMath::Sqrt(0.5*(TMath::Sqrt(c1*c1+c2)+c1));
1266 Float_t b2=0.5*(TMath::Sqrt(c1*c1+c2)-c1);
1267
1268 Float_t rs=((aO-pdoti)*(aO-pdoti)+b2)/((aO+pdoti)*(aO+pdoti)+b2);
1269 Float_t rp=rs*((aO-sinin*tanin)*(aO-sinin*tanin)+b2)/((aO+sinin*tanin)*(aO+sinin*tanin)+b2);
1270
1271
1272 //CORRECTION FACTOR FOR SURFACE ROUGHNESS
1273 //B.J. STAGG APPLIED OPTICS, 30(1991),4113
1274
1275 Float_t sigraf=18.;
1276 Float_t lamb=1240/ene;
1277 Float_t fresn;
1278
1279 Float_t rO=TMath::Exp(-(4*TMath::Pi()*pdoti*sigraf/lamb)*(4*TMath::Pi()*pdoti*sigraf/lamb));
1280
1281 if(pola)
1282 {
1283 Float_t pdotr=0.8; //DEGREE OF POLARIZATION : 1->P , -1->S
1284 fresn=0.5*(rp*(1+pdotr)+rs*(1-pdotr));
1285 }
1286 else
1287 fresn=0.5*(rp+rs);
1288
1289 fresn = fresn*rO;
1290 return(fresn);
1291}
1292
1293//__________________________________________
1294Float_t AliRICH::AbsoCH4(Float_t x)
1295{
1296
1297 //KLOSCH,SCH4(9),WL(9),EM(9),ALENGTH(31)
1298 Float_t sch4[9] = {.12,.16,.23,.38,.86,2.8,7.9,28.,80.}; //MB X 10^22
1299 //Float_t wl[9] = {153.,152.,151.,150.,149.,148.,147.,146.,145};
1300 Float_t em[9] = {8.1,8.158,8.212,8.267,8.322,8.378,8.435,8.493,8.55};
1301 const Float_t kLosch=2.686763E19; // LOSCHMIDT NUMBER IN CM-3
1302 const Float_t kIgas1=100, kIgas2=0, kOxy=10., kWater=5., kPressure=750.,kTemperature=283.;
1303 Float_t pn=kPressure/760.;
1304 Float_t tn=kTemperature/273.16;
1305
1306
1307// ------- METHANE CROSS SECTION -----------------
1308// ASTROPH. J. 214, L47 (1978)
1309
1310 Float_t sm=0;
1311 if (x<7.75)
1312 sm=.06e-22;
1313
1314 if(x>=7.75 && x<=8.1)
1315 {
1316 Float_t c0=-1.655279e-1;
1317 Float_t c1=6.307392e-2;
1318 Float_t c2=-8.011441e-3;
1319 Float_t c3=3.392126e-4;
1320 sm=(c0+c1*x+c2*x*x+c3*x*x*x)*1.e-18;
1321 }
1322
1323 if (x> 8.1)
1324 {
1325 Int_t j=0;
1326 while (x<=em[j] && x>=em[j+1])
1327 {
1328 j++;
1329 Float_t a=(sch4[j+1]-sch4[j])/(em[j+1]-em[j]);
1330 sm=(sch4[j]+a*(x-em[j]))*1e-22;
1331 }
1332 }
1333
1334 Float_t dm=(kIgas1/100.)*(1.-((kOxy+kWater)/1.e6))*kLosch*pn/tn;
1335 Float_t abslm=1./sm/dm;
1336
1337// ------- ISOBUTHANE CROSS SECTION --------------
1338// i-C4H10 (ai) abs. length from curves in
1339// Lu-McDonald paper for BARI RICH workshop .
1340// -----------------------------------------------------------
1341
1342 Float_t ai;
1343 Float_t absli;
1344 if (kIgas2 != 0)
1345 {
1346 if (x<7.25)
1347 ai=100000000.;
1348
1349 if(x>=7.25 && x<7.375)
1350 ai=24.3;
1351
1352 if(x>=7.375)
1353 ai=.0000000001;
1354
1355 Float_t si = 1./(ai*kLosch*273.16/293.); // ISOB. CRO.SEC.IN CM2
1356 Float_t di=(kIgas2/100.)*(1.-((kOxy+kWater)/1.e6))*kLosch*pn/tn;
1357 absli =1./si/di;
1358 }
1359 else
1360 absli=1.e18;
1361// ---------------------------------------------------------
1362//
1363// transmission of O2
1364//
1365// y= path in cm, x=energy in eV
1366// so= cross section for UV absorption in cm2
1367// do= O2 molecular density in cm-3
1368// ---------------------------------------------------------
1369
1370 Float_t abslo;
1371 Float_t so=0;
1372 if(x>=6.0)
1373 {
1374 if(x>=6.0 && x<6.5)
1375 {
1376 so=3.392709e-13 * TMath::Exp(2.864104 *x);
1377 so=so*1e-18;
1378 }
1379
1380 if(x>=6.5 && x<7.0)
1381 {
1382 so=2.910039e-34 * TMath::Exp(10.3337*x);
1383 so=so*1e-18;
1384 }
1385
1386
1387 if (x>=7.0)
1388 {
1389 Float_t a0=-73770.76;
1390 Float_t a1=46190.69;
1391 Float_t a2=-11475.44;
1392 Float_t a3=1412.611;
1393 Float_t a4=-86.07027;
1394 Float_t a5=2.074234;
1395 so= a0+(a1*x)+(a2*x*x)+(a3*x*x*x)+(a4*x*x*x*x)+(a5*x*x*x*x*x);
1396 so=so*1e-18;
1397 }
1398
1399 Float_t dox=(kOxy/1e6)*kLosch*pn/tn;
1400 abslo=1./so/dox;
1401 }
1402 else
1403 abslo=1.e18;
1404// ---------------------------------------------------------
1405//
1406// transmission of H2O
1407//
1408// y= path in cm, x=energy in eV
1409// sw= cross section for UV absorption in cm2
1410// dw= H2O molecular density in cm-3
1411// ---------------------------------------------------------
1412
1413 Float_t abslw;
1414
1415 Float_t b0=29231.65;
1416 Float_t b1=-15807.74;
1417 Float_t b2=3192.926;
1418 Float_t b3=-285.4809;
1419 Float_t b4=9.533944;
1420
1421 if(x>6.75)
1422 {
1423 Float_t sw= b0+(b1*x)+(b2*x*x)+(b3*x*x*x)+(b4*x*x*x*x);
1424 sw=sw*1e-18;
1425 Float_t dw=(kWater/1e6)*kLosch*pn/tn;
1426 abslw=1./sw/dw;
1427 }
1428 else
1429 abslw=1.e18;
1430
1431// ---------------------------------------------------------
1432
1433 Float_t alength=1./(1./abslm+1./absli+1./abslo+1./abslw);
1434 return (alength);
1435}
1436
1437
1438
ddae0931 1439//___________________________________________
1440Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
1441{
237c933d 1442
1443// Default value
1444
ddae0931 1445 return 9999;
1446}
fe4da5cc 1447
ddae0931 1448//___________________________________________
1449void AliRICH::MakeBranch(Option_t* option)
fe4da5cc 1450{
237c933d 1451 // Create Tree branches for the RICH.
fe4da5cc 1452
237c933d 1453 const Int_t kBufferSize = 4000;
ddae0931 1454 char branchname[20];
fe4da5cc 1455
ddae0931 1456
1457 AliDetector::MakeBranch(option);
1458 sprintf(branchname,"%sCerenkov",GetName());
1459 if (fCerenkovs && gAlice->TreeH()) {
237c933d 1460 gAlice->TreeH()->Branch(branchname,&fCerenkovs, kBufferSize);
ddae0931 1461 printf("Making Branch %s for Cerenkov Hits\n",branchname);
fe4da5cc 1462 }
ddae0931 1463
2e5f0f7b 1464 sprintf(branchname,"%sPadHits",GetName());
1465 if (fPadHits && gAlice->TreeH()) {
237c933d 1466 gAlice->TreeH()->Branch(branchname,&fPadHits, kBufferSize);
2e5f0f7b 1467 printf("Making Branch %s for PadHits\n",branchname);
fe4da5cc 1468 }
1469
ddae0931 1470// one branch for digits per chamber
1471 Int_t i;
1472
237c933d 1473 for (i=0; i<kNCH ;i++) {
ddae0931 1474 sprintf(branchname,"%sDigits%d",GetName(),i+1);
1475
1476 if (fDchambers && gAlice->TreeD()) {
237c933d 1477 gAlice->TreeD()->Branch(branchname,&((*fDchambers)[i]), kBufferSize);
ddae0931 1478 printf("Making Branch %s for digits in chamber %d\n",branchname,i+1);
1479 }
fe4da5cc 1480 }
2e5f0f7b 1481
1482// one branch for raw clusters per chamber
237c933d 1483 for (i=0; i<kNCH ;i++) {
2e5f0f7b 1484 sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
1485
1486 if (fRawClusters && gAlice->TreeR()) {
237c933d 1487 gAlice->TreeR()->Branch(branchname,&((*fRawClusters)[i]), kBufferSize);
2e5f0f7b 1488 printf("Making Branch %s for raw clusters in chamber %d\n",branchname,i+1);
1489 }
1490 }
1491
1492 // one branch for rec hits per chamber
237c933d 1493 for (i=0; i<kNCH ;i++) {
a4622d0f 1494 sprintf(branchname,"%sRecHits1D%d",GetName(),i+1);
1495
1496 if (fRecHits1D && gAlice->TreeR()) {
1497 gAlice->TreeR()->Branch(branchname,&((*fRecHits1D)[i]), kBufferSize);
1498 printf("Making Branch %s for 1D rec. hits in chamber %d\n",branchname,i+1);
1499 }
1500 }
1501 for (i=0; i<kNCH ;i++) {
1502 sprintf(branchname,"%sRecHits3D%d",GetName(),i+1);
2e5f0f7b 1503
a4622d0f 1504 if (fRecHits3D && gAlice->TreeR()) {
1505 gAlice->TreeR()->Branch(branchname,&((*fRecHits3D)[i]), kBufferSize);
1506 printf("Making Branch %s for 3D rec. hits in chamber %d\n",branchname,i+1);
2e5f0f7b 1507 }
1508 }
a4622d0f 1509
fe4da5cc 1510}
1511
ddae0931 1512//___________________________________________
1513void AliRICH::SetTreeAddress()
fe4da5cc 1514{
237c933d 1515 // Set branch address for the Hits and Digits Tree.
2e5f0f7b 1516 char branchname[20];
1517 Int_t i;
1518
ddae0931 1519 AliDetector::SetTreeAddress();
1520
1521 TBranch *branch;
1522 TTree *treeH = gAlice->TreeH();
1523 TTree *treeD = gAlice->TreeD();
2e5f0f7b 1524 TTree *treeR = gAlice->TreeR();
ddae0931 1525
1526 if (treeH) {
2e5f0f7b 1527 if (fPadHits) {
1528 branch = treeH->GetBranch("RICHPadHits");
1529 if (branch) branch->SetAddress(&fPadHits);
fe4da5cc 1530 }
ddae0931 1531 if (fCerenkovs) {
1532 branch = treeH->GetBranch("RICHCerenkov");
1533 if (branch) branch->SetAddress(&fCerenkovs);
fe4da5cc 1534 }
ddae0931 1535 }
1536
1537 if (treeD) {
237c933d 1538 for (int i=0; i<kNCH; i++) {
ddae0931 1539 sprintf(branchname,"%sDigits%d",GetName(),i+1);
1540 if (fDchambers) {
1541 branch = treeD->GetBranch(branchname);
1542 if (branch) branch->SetAddress(&((*fDchambers)[i]));
1543 }
fe4da5cc 1544 }
fe4da5cc 1545 }
2e5f0f7b 1546 if (treeR) {
237c933d 1547 for (i=0; i<kNCH; i++) {
2e5f0f7b 1548 sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
1549 if (fRawClusters) {
1550 branch = treeR->GetBranch(branchname);
1551 if (branch) branch->SetAddress(&((*fRawClusters)[i]));
1552 }
1553 }
1554
237c933d 1555 for (i=0; i<kNCH; i++) {
a4622d0f 1556 sprintf(branchname,"%sRecHits1D%d",GetName(),i+1);
1557 if (fRecHits1D) {
2e5f0f7b 1558 branch = treeR->GetBranch(branchname);
a4622d0f 1559 if (branch) branch->SetAddress(&((*fRecHits1D)[i]));
2e5f0f7b 1560 }
1561 }
1562
a4622d0f 1563 for (i=0; i<kNCH; i++) {
1564 sprintf(branchname,"%sRecHits3D%d",GetName(),i+1);
1565 if (fRecHits3D) {
1566 branch = treeR->GetBranch(branchname);
1567 if (branch) branch->SetAddress(&((*fRecHits3D)[i]));
1568 }
1569 }
1570
2e5f0f7b 1571 }
ddae0931 1572}
1573//___________________________________________
1574void AliRICH::ResetHits()
1575{
237c933d 1576 // Reset number of clusters and the cluster array for this detector
ddae0931 1577 AliDetector::ResetHits();
2e5f0f7b 1578 fNPadHits = 0;
1579 fNcerenkovs = 0;
1580 if (fPadHits) fPadHits->Clear();
ddae0931 1581 if (fCerenkovs) fCerenkovs->Clear();
fe4da5cc 1582}
1583
2e5f0f7b 1584
ddae0931 1585//____________________________________________
1586void AliRICH::ResetDigits()
1587{
237c933d 1588 //
1589 // Reset number of digits and the digits array for this detector
1590 //
1591 for ( int i=0;i<kNCH;i++ ) {
ddae0931 1592 if ((*fDchambers)[i]) (*fDchambers)[i]->Clear();
1593 if (fNdch) fNdch[i]=0;
1594 }
1595}
2e5f0f7b 1596
ddae0931 1597//____________________________________________
2e5f0f7b 1598void AliRICH::ResetRawClusters()
fe4da5cc 1599{
237c933d 1600 //
1601 // Reset number of raw clusters and the raw clust array for this detector
1602 //
1603 for ( int i=0;i<kNCH;i++ ) {
2e5f0f7b 1604 if ((*fRawClusters)[i]) ((TClonesArray*)(*fRawClusters)[i])->Clear();
1605 if (fNrawch) fNrawch[i]=0;
fe4da5cc 1606 }
ddae0931 1607}
fe4da5cc 1608
2e5f0f7b 1609//____________________________________________
a4622d0f 1610void AliRICH::ResetRecHits1D()
fe4da5cc 1611{
237c933d 1612 //
1613 // Reset number of raw clusters and the raw clust array for this detector
1614 //
2e5f0f7b 1615
237c933d 1616 for ( int i=0;i<kNCH;i++ ) {
a4622d0f 1617 if ((*fRecHits1D)[i]) ((TClonesArray*)(*fRecHits1D)[i])->Clear();
1618 if (fNrechits1D) fNrechits1D[i]=0;
1619 }
1620}
1621
1622//____________________________________________
1623void AliRICH::ResetRecHits3D()
1624{
1625 //
1626 // Reset number of raw clusters and the raw clust array for this detector
1627 //
1628
1629 for ( int i=0;i<kNCH;i++ ) {
1630 if ((*fRecHits3D)[i]) ((TClonesArray*)(*fRecHits3D)[i])->Clear();
1631 if (fNrechits3D) fNrechits3D[i]=0;
2e5f0f7b 1632 }
fe4da5cc 1633}
1634
ddae0931 1635//___________________________________________
2e5f0f7b 1636void AliRICH::SetGeometryModel(Int_t id, AliRICHGeometry *geometry)
fe4da5cc 1637{
237c933d 1638
1639//
1640// Setter for the RICH geometry model
1641//
1642
1643
2e5f0f7b 1644 ((AliRICHChamber*) (*fChambers)[id])->GeometryModel(geometry);
fe4da5cc 1645}
1646
ddae0931 1647//___________________________________________
a2f7eaf6 1648void AliRICH::SetSegmentationModel(Int_t id, AliSegmentation *segmentation)
ddae0931 1649{
237c933d 1650
1651//
1652// Setter for the RICH segmentation model
1653//
1654
a2f7eaf6 1655 ((AliRICHChamber*) (*fChambers)[id])->SetSegmentationModel(segmentation);
ddae0931 1656}
fe4da5cc 1657
ddae0931 1658//___________________________________________
2e5f0f7b 1659void AliRICH::SetResponseModel(Int_t id, AliRICHResponse *response)
fe4da5cc 1660{
237c933d 1661
1662//
1663// Setter for the RICH response model
1664//
1665
2e5f0f7b 1666 ((AliRICHChamber*) (*fChambers)[id])->ResponseModel(response);
fe4da5cc 1667}
1668
2e5f0f7b 1669void AliRICH::SetReconstructionModel(Int_t id, AliRICHClusterFinder *reconst)
ddae0931 1670{
237c933d 1671
1672//
1673// Setter for the RICH reconstruction model (clusters)
1674//
1675
a2f7eaf6 1676 ((AliRICHChamber*) (*fChambers)[id])->SetReconstructionModel(reconst);
fe4da5cc 1677}
1678
ddae0931 1679void AliRICH::SetNsec(Int_t id, Int_t nsec)
1680{
237c933d 1681
1682//
1683// Sets the number of padplanes
1684//
1685
2e5f0f7b 1686 ((AliRICHChamber*) (*fChambers)[id])->SetNsec(nsec);
ddae0931 1687}
1688
1689
1690//___________________________________________
ddae0931 1691void AliRICH::StepManager()
fe4da5cc 1692{
237c933d 1693
452a64c6 1694// Full Step Manager
1695
1696 Int_t copy, id;
1697 static Int_t idvol;
1698 static Int_t vol[2];
1699 Int_t ipart;
1700 static Float_t hits[18];
1701 static Float_t ckovData[19];
1702 TLorentzVector position;
1703 TLorentzVector momentum;
1704 Float_t pos[3];
1705 Float_t mom[4];
1706 Float_t localPos[3];
1707 Float_t localMom[4];
1708 Float_t localTheta,localPhi;
1709 Float_t theta,phi;
1710 Float_t destep, step;
1711 Float_t ranf[2];
1712 Int_t nPads;
1713 Float_t coscerenkov;
1714 static Float_t eloss, xhit, yhit, tlength;
1715 const Float_t kBig=1.e10;
1716
1717 TClonesArray &lhits = *fHits;
452a64c6 1718 TParticle *current = (TParticle*)(*gAlice->Particles())[gAlice->CurrentTrack()];
237c933d 1719
452a64c6 1720 //if (current->Energy()>1)
1721 //{
1722
1723 // Only gas gap inside chamber
1724 // Tag chambers and record hits when track enters
1725
1726 idvol=-1;
1727 id=gMC->CurrentVolID(copy);
1728 Float_t cherenkovLoss=0;
1729 //gAlice->KeepTrack(gAlice->CurrentTrack());
1730
1731 gMC->TrackPosition(position);
1732 pos[0]=position(0);
1733 pos[1]=position(1);
1734 pos[2]=position(2);
8cda28a5 1735 //bzero((char *)ckovData,sizeof(ckovData)*19);
452a64c6 1736 ckovData[1] = pos[0]; // X-position for hit
1737 ckovData[2] = pos[1]; // Y-position for hit
1738 ckovData[3] = pos[2]; // Z-position for hit
3cc8edee 1739 ckovData[6] = 0; // dummy track length
452a64c6 1740 //ckovData[11] = gAlice->CurrentTrack();
8cda28a5 1741
1742 //printf("\n+++++++++++\nTrack: %d\n++++++++++++\n",gAlice->CurrentTrack());
452a64c6 1743
1744 //AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH");
1745
1746 /********************Store production parameters for Cerenkov photons************************/
1747//is it a Cerenkov photon?
8cda28a5 1748 if (gMC->TrackPid() == 50000050) {
452a64c6 1749
1750 //if (gMC->VolId("GAP ")==gMC->CurrentVolID(copy))
1751 //{
1752 Float_t ckovEnergy = current->Energy();
1753 //energy interval for tracking
1754 if (ckovEnergy > 5.6e-09 && ckovEnergy < 7.8e-09 )
1755 //if (ckovEnergy > 0)
1756 {
8cda28a5 1757 if (gMC->IsTrackEntering()){ //is track entering?
1758 //printf("Track entered (1)\n");
452a64c6 1759 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
1760 { //is it in freo?
5cfcdf54 1761 if (gMC->IsNewTrack()){ //is it the first step?
8cda28a5 1762 //printf("I'm in!\n");
452a64c6 1763 Int_t mother = current->GetFirstMother();
1764
1765 //printf("Second Mother:%d\n",current->GetSecondMother());
1766
1767 ckovData[10] = mother;
1768 ckovData[11] = gAlice->CurrentTrack();
1769 ckovData[12] = 1; //Media where photon was produced 1->Freon, 2->Quarz
8cda28a5 1770 //printf("Produced in FREO\n");
452a64c6 1771 fCkovNumber++;
1772 fFreonProd=1;
1773 //printf("Index: %d\n",fCkovNumber);
1774 } //first step question
1775 } //freo question
1776
5cfcdf54 1777 if (gMC->IsNewTrack()){ //is it first step?
452a64c6 1778 if (gMC->VolId("QUAR")==gMC->CurrentVolID(copy)) //is it in quarz?
1779 {
1780 ckovData[12] = 2;
8cda28a5 1781 //printf("Produced in QUAR\n");
452a64c6 1782 } //quarz question
1783 } //first step question
1784
1785 //printf("Before %d\n",fFreonProd);
1786 } //track entering question
1787
1788 if (ckovData[12] == 1) //was it produced in Freon?
1789 //if (fFreonProd == 1)
1790 {
1791 if (gMC->IsTrackEntering()){ //is track entering?
8cda28a5 1792 //printf("Track entered (2)\n");
1793 //printf("Current volume (should be META): %s\n",gMC->CurrentVolName());
1794 //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("META"),gMC->CurrentVolID(copy));
452a64c6 1795 if (gMC->VolId("META")==gMC->CurrentVolID(copy)) //is it in gap?
1796 {
8cda28a5 1797 //printf("Got in META\n");
452a64c6 1798 gMC->TrackMomentum(momentum);
1799 mom[0]=momentum(0);
1800 mom[1]=momentum(1);
1801 mom[2]=momentum(2);
1802 mom[3]=momentum(3);
1803 // Z-position for hit
1804
1805
1806 /**************** Photons lost in second grid have to be calculated by hand************/
1807
1808 Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
1809 Float_t t = (1. - .025 / cophi) * (1. - .05 / cophi);
1810 gMC->Rndm(ranf, 1);
1811 //printf("grid calculation:%f\n",t);
1812 if (ranf[0] > t) {
5cfcdf54 1813 gMC->StopTrack();
452a64c6 1814 ckovData[13] = 5;
1815 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
8cda28a5 1816 //printf("Added One (1)!\n");
452a64c6 1817 //printf("Lost one in grid\n");
1818 }
1819 /**********************************************************************************/
1820 } //gap
1821
8cda28a5 1822 //printf("Current volume (should be CSI) (1): %s\n",gMC->CurrentVolName());
1823 //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
452a64c6 1824 if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy)) //is it in csi?
1825 {
8cda28a5 1826 //printf("Got in CSI\n");
452a64c6 1827 gMC->TrackMomentum(momentum);
1828 mom[0]=momentum(0);
1829 mom[1]=momentum(1);
1830 mom[2]=momentum(2);
1831 mom[3]=momentum(3);
1832
1833 /********* Photons lost by Fresnel reflection have to be calculated by hand********/
1834 /***********************Cerenkov phtons (always polarised)*************************/
1835
1836 Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
1837 Float_t t = Fresnel(ckovEnergy*1e9,cophi,1);
1838 gMC->Rndm(ranf, 1);
1839 if (ranf[0] < t) {
5cfcdf54 1840 gMC->StopTrack();
452a64c6 1841 ckovData[13] = 6;
1842 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
8cda28a5 1843 //printf("Added One (2)!\n");
452a64c6 1844 //printf("Lost by Fresnel\n");
1845 }
1846 /**********************************************************************************/
1847 }
1848 } //track entering?
1849
1850
1851 /********************Evaluation of losses************************/
1852 /******************still in the old fashion**********************/
1853
5cfcdf54 1854 TArrayI procs;
1855 Int_t i1 = gMC->StepProcesses(procs); //number of physics mechanisms acting on the particle
452a64c6 1856 for (Int_t i = 0; i < i1; ++i) {
1857 // Reflection loss
5cfcdf54 1858 if (procs[i] == kPLightReflection) { //was it reflected
452a64c6 1859 ckovData[13]=10;
1860 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
1861 ckovData[13]=1;
1862 if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
1863 ckovData[13]=2;
5cfcdf54 1864 //gMC->StopTrack();
9dda3582 1865 //AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
452a64c6 1866 } //reflection question
8cda28a5 1867
452a64c6 1868 // Absorption loss
5cfcdf54 1869 else if (procs[i] == kPLightAbsorption) { //was it absorbed?
8cda28a5 1870 //printf("Got in absorption\n");
452a64c6 1871 ckovData[13]=20;
1872 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
1873 ckovData[13]=11;
1874 if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
1875 ckovData[13]=12;
1876 if (gMC->CurrentVolID(copy) == gMC->VolId("META"))
1877 ckovData[13]=13;
1878 if (gMC->CurrentVolID(copy) == gMC->VolId("GAP "))
1879 ckovData[13]=13;
1880
1881 if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC"))
1882 ckovData[13]=15;
1883
1884 // CsI inefficiency
1885 if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
1886 ckovData[13]=16;
1887 }
5cfcdf54 1888 gMC->StopTrack();
452a64c6 1889 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
8cda28a5 1890 //printf("Added One (3)!\n");
452a64c6 1891 //printf("Added cerenkov %d\n",fCkovNumber);
1892 } //absorption question
1893
1894
1895 // Photon goes out of tracking scope
5cfcdf54 1896 else if (procs[i] == kPStop) { //is it below energy treshold?
452a64c6 1897 ckovData[13]=21;
5cfcdf54 1898 gMC->StopTrack();
452a64c6 1899 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
8cda28a5 1900 //printf("Added One (4)!\n");
452a64c6 1901 } // energy treshold question
1902 } //number of mechanisms cycle
1903 /**********************End of evaluation************************/
1904 } //freon production question
1905 } //energy interval question
1906 //}//inside the proximity gap question
1907 } //cerenkov photon question
1908
1909 /**************************************End of Production Parameters Storing*********************/
1910
1911
1912 /*******************************Treat photons that hit the CsI (Ckovs and Feedbacks)************/
1913
1914 if (gMC->TrackPid() == 50000050 || gMC->TrackPid() == 50000051) {
1915 //printf("Cerenkov\n");
1916 if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
1917 {
8cda28a5 1918 //printf("Current volume (should be CSI) (2): %s\n",gMC->CurrentVolName());
1919 //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
1920 //printf("Got in CSI\n");
452a64c6 1921 if (gMC->Edep() > 0.){
1922 gMC->TrackPosition(position);
1923 gMC->TrackMomentum(momentum);
1924 pos[0]=position(0);
1925 pos[1]=position(1);
1926 pos[2]=position(2);
1927 mom[0]=momentum(0);
1928 mom[1]=momentum(1);
1929 mom[2]=momentum(2);
1930 mom[3]=momentum(3);
1931 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
1932 Double_t rt = TMath::Sqrt(tc);
1933 theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
1934 phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
1935 gMC->Gmtod(pos,localPos,1);
1936 gMC->Gmtod(mom,localMom,2);
1937
1938 gMC->CurrentVolOffID(2,copy);
1939 vol[0]=copy;
1940 idvol=vol[0]-1;
1941
1942 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
1943 //->Sector(localPos[0], localPos[2]);
1944 //printf("Sector:%d\n",sector);
1945
1946 /*if (gMC->TrackPid() == 50000051){
1947 fFeedbacks++;
1948 printf("Feedbacks:%d\n",fFeedbacks);
1949 }*/
1950
1951 ((AliRICHChamber*) (*fChambers)[idvol])
1952 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
1953 if(idvol<kNCH) {
1954 ckovData[0] = gMC->TrackPid(); // particle type
1955 ckovData[1] = pos[0]; // X-position for hit
1956 ckovData[2] = pos[1]; // Y-position for hit
1957 ckovData[3] = pos[2]; // Z-position for hit
1958 ckovData[4] = theta; // theta angle of incidence
1959 ckovData[5] = phi; // phi angle of incidence
1960 ckovData[8] = (Float_t) fNPadHits; // first padhit
1961 ckovData[9] = -1; // last pad hit
1962 ckovData[13] = 4; // photon was detected
1963 ckovData[14] = mom[0];
1964 ckovData[15] = mom[1];
1965 ckovData[16] = mom[2];
1966
1967 destep = gMC->Edep();
1968 gMC->SetMaxStep(kBig);
1969 cherenkovLoss += destep;
1970 ckovData[7]=cherenkovLoss;
1971
1972 nPads = MakePadHits(localPos[0],localPos[2],cherenkovLoss,idvol,kCerenkov);
1973 if (fNPadHits > (Int_t)ckovData[8]) {
1974 ckovData[8]= ckovData[8]+1;
1975 ckovData[9]= (Float_t) fNPadHits;
1976 }
1977
1978 ckovData[17] = nPads;
1979 //printf("nPads:%d",nPads);
1980
1981 //TClonesArray *Hits = RICH->Hits();
1982 AliRICHHit *mipHit = (AliRICHHit*) (fHits->UncheckedAt(0));
1983 if (mipHit)
1984 {
1985 mom[0] = current->Px();
1986 mom[1] = current->Py();
1987 mom[2] = current->Pz();
1988 Float_t mipPx = mipHit->fMomX;
1989 Float_t mipPy = mipHit->fMomY;
1990 Float_t mipPz = mipHit->fMomZ;
1991
1992 Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
1993 Float_t rt = TMath::Sqrt(r);
1994 Float_t mipR = mipPx*mipPx + mipPy*mipPy + mipPz*mipPz;
1995 Float_t mipRt = TMath::Sqrt(mipR);
1996 if ((rt*mipRt) > 0)
1997 {
1998 coscerenkov = (mom[0]*mipPx + mom[1]*mipPy + mom[2]*mipPz)/(rt*mipRt);
1999 }
2000 else
2001 {
2002 coscerenkov = 0;
2003 }
2004 Float_t cherenkov = TMath::ACos(coscerenkov);
2005 ckovData[18]=cherenkov;
2006 }
2007 //if (sector != -1)
2008 //{
2009 AddHit(gAlice->CurrentTrack(),vol,ckovData);
2010 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
8cda28a5 2011 //printf("Added One (5)!\n");
452a64c6 2012 //}
2013 }
2014 }
2015 }
2016 }
2017
2018 /***********************************************End of photon hits*********************************************/
2019
2020
2021 /**********************************************Charged particles treatment*************************************/
2022
2023 else if (gMC->TrackCharge())
2024 //else if (1 == 1)
2025 {
2026//If MIP
2027 /*if (gMC->IsTrackEntering())
2028 {
2029 hits[13]=20;//is track entering?
2030 }*/
2031 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2032 {
2033 fFreonProd=1;
2034 }
2035
2036 if (gMC->VolId("GAP ")== gMC->CurrentVolID(copy)) {
2037// Get current particle id (ipart), track position (pos) and momentum (mom)
2038
2039 gMC->CurrentVolOffID(3,copy);
2040 vol[0]=copy;
2041 idvol=vol[0]-1;
2042
2043 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
2044 //->Sector(localPos[0], localPos[2]);
2045 //printf("Sector:%d\n",sector);
2046
2047 gMC->TrackPosition(position);
2048 gMC->TrackMomentum(momentum);
2049 pos[0]=position(0);
2050 pos[1]=position(1);
2051 pos[2]=position(2);
2052 mom[0]=momentum(0);
2053 mom[1]=momentum(1);
2054 mom[2]=momentum(2);
2055 mom[3]=momentum(3);
2056 gMC->Gmtod(pos,localPos,1);
2057 gMC->Gmtod(mom,localMom,2);
2058
2059 ipart = gMC->TrackPid();
2060 //
2061 // momentum loss and steplength in last step
2062 destep = gMC->Edep();
2063 step = gMC->TrackStep();
2064
2065 //
2066 // record hits when track enters ...
2067 if( gMC->IsTrackEntering()) {
2068// gMC->SetMaxStep(fMaxStepGas);
2069 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
2070 Double_t rt = TMath::Sqrt(tc);
2071 theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
2072 phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
2073
2074
2075 Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
2076 Double_t localRt = TMath::Sqrt(localTc);
2077 localTheta = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])))*kRaddeg;
2078 localPhi = Float_t(TMath::ATan2(Double_t(localMom[2]),Double_t(localMom[0])))*kRaddeg;
2079
2080 hits[0] = Float_t(ipart); // particle type
2081 hits[1] = localPos[0]; // X-position for hit
2082 hits[2] = localPos[1]; // Y-position for hit
2083 hits[3] = localPos[2]; // Z-position for hit
2084 hits[4] = localTheta; // theta angle of incidence
2085 hits[5] = localPhi; // phi angle of incidence
2086 hits[8] = (Float_t) fNPadHits; // first padhit
2087 hits[9] = -1; // last pad hit
2088 hits[13] = fFreonProd; // did id hit the freon?
2089 hits[14] = mom[0];
2090 hits[15] = mom[1];
2091 hits[16] = mom[2];
2092
2093 tlength = 0;
2094 eloss = 0;
2095 fFreonProd = 0;
2096
2097 Chamber(idvol).LocaltoGlobal(localPos,hits+1);
2098
2099
2100 //To make chamber coordinates x-y had to pass localPos[0], localPos[2]
2101 xhit = localPos[0];
2102 yhit = localPos[2];
2103 // Only if not trigger chamber
2104 if(idvol<kNCH) {
2105 //
2106 // Initialize hit position (cursor) in the segmentation model
2107 ((AliRICHChamber*) (*fChambers)[idvol])
2108 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2109 }
2110 }
2111
2112 //
2113 // Calculate the charge induced on a pad (disintegration) in case
2114 //
2115 // Mip left chamber ...
2116 if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
2117 gMC->SetMaxStep(kBig);
2118 eloss += destep;
2119 tlength += step;
2120
2121
2122 // Only if not trigger chamber
2123 if(idvol<kNCH) {
2124 if (eloss > 0)
2125 {
2126 if(gMC->TrackPid() == kNeutron)
2127 printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
2128 nPads = MakePadHits(xhit,yhit,eloss,idvol,kMip);
2129 hits[17] = nPads;
2130 //printf("nPads:%d",nPads);
2131 }
2132 }
2133
2134 hits[6]=tlength;
2135 hits[7]=eloss;
2136 if (fNPadHits > (Int_t)hits[8]) {
2137 hits[8]= hits[8]+1;
2138 hits[9]= (Float_t) fNPadHits;
2139 }
2140
2141 //if(sector !=-1)
2142 new(lhits[fNhits++]) AliRICHHit(fIshunt,gAlice->CurrentTrack(),vol,hits);
2143 eloss = 0;
2144 //
2145 // Check additional signal generation conditions
2146 // defined by the segmentation
2147 // model (boundary crossing conditions)
2148 } else if
2149 (((AliRICHChamber*) (*fChambers)[idvol])
2150 ->SigGenCond(localPos[0], localPos[2], localPos[1]))
2151 {
2152 ((AliRICHChamber*) (*fChambers)[idvol])
2153 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2154 if (eloss > 0)
2155 {
2156 if(gMC->TrackPid() == kNeutron)
2157 printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
2158 nPads = MakePadHits(xhit,yhit,eloss,idvol,kMip);
2159 hits[17] = nPads;
2160 //printf("Npads:%d",NPads);
2161 }
2162 xhit = localPos[0];
2163 yhit = localPos[2];
2164 eloss = destep;
2165 tlength += step ;
2166 //
2167 // nothing special happened, add up energy loss
2168 } else {
2169 eloss += destep;
2170 tlength += step ;
2171 }
2172 }
2173 }
2174 /*************************************************End of MIP treatment**************************************/
2175 //}
ddae0931 2176}
2e5f0f7b 2177
237c933d 2178void AliRICH::FindClusters(Int_t nev,Int_t lastEntry)
ddae0931 2179{
2e5f0f7b 2180
ddae0931 2181//
2182// Loop on chambers and on cathode planes
2183//
2e5f0f7b 2184 for (Int_t icat=1;icat<2;icat++) {
2185 gAlice->ResetDigits();
2186 gAlice->TreeD()->GetEvent(1); // spurious +1 ...
237c933d 2187 for (Int_t ich=0;ich<kNCH;ich++) {
2e5f0f7b 2188 AliRICHChamber* iChamber=(AliRICHChamber*) (*fChambers)[ich];
237c933d 2189 TClonesArray *pRICHdigits = this->DigitsAddress(ich);
2190 if (pRICHdigits == 0)
2e5f0f7b 2191 continue;
2192 //
2193 // Get ready the current chamber stuff
2194 //
2195 AliRICHResponse* response = iChamber->GetResponseModel();
a2f7eaf6 2196 AliSegmentation* seg = iChamber->GetSegmentationModel();
2e5f0f7b 2197 AliRICHClusterFinder* rec = iChamber->GetReconstructionModel();
2198 if (seg) {
2199 rec->SetSegmentation(seg);
2200 rec->SetResponse(response);
237c933d 2201 rec->SetDigits(pRICHdigits);
2e5f0f7b 2202 rec->SetChamber(ich);
2203 if (nev==0) rec->CalibrateCOG();
2204 rec->FindRawClusters();
2205 }
2206 TClonesArray *fRch;
2207 fRch=RawClustAddress(ich);
2208 fRch->Sort();
2209 } // for ich
2210
2211 gAlice->TreeR()->Fill();
2212 TClonesArray *fRch;
237c933d 2213 for (int i=0;i<kNCH;i++) {
2e5f0f7b 2214 fRch=RawClustAddress(i);
2215 int nraw=fRch->GetEntriesFast();
2216 printf ("Chamber %d, raw clusters %d\n",i,nraw);
2217 }
2218
2219 ResetRawClusters();
2220
2221 } // for icat
2222
2223 char hname[30];
2224 sprintf(hname,"TreeR%d",nev);
33984590 2225 gAlice->TreeR()->Write(hname,kOverwrite,0);
2e5f0f7b 2226 gAlice->TreeR()->Reset();
2227
2228 //gObjectTable->Print();
ddae0931 2229}
2230
2231
2232//______________________________________________________________________________
2233void AliRICH::Streamer(TBuffer &R__b)
2234{
2235 // Stream an object of class AliRICH.
2e5f0f7b 2236 AliRICHChamber *iChamber;
a2f7eaf6 2237 AliSegmentation *segmentation;
2e5f0f7b 2238 AliRICHResponse *response;
ddae0931 2239 TClonesArray *digitsaddress;
2e5f0f7b 2240 TClonesArray *rawcladdress;
a4622d0f 2241 TClonesArray *rechitaddress1D;
2242 TClonesArray *rechitaddress3D;
ddae0931 2243
2244 if (R__b.IsReading()) {
2245 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
2246 AliDetector::Streamer(R__b);
2e5f0f7b 2247 R__b >> fNPadHits;
2248 R__b >> fPadHits; // diff
2249 R__b >> fNcerenkovs;
2250 R__b >> fCerenkovs; // diff
ddae0931 2251 R__b >> fDchambers;
2e5f0f7b 2252 R__b >> fRawClusters;
a4622d0f 2253 R__b >> fRecHits1D; //diff
2254 R__b >> fRecHits3D; //diff
33984590 2255 R__b >> fDebugLevel; //diff
9dda3582 2256 R__b.ReadStaticArray(fNdch);
2257 R__b.ReadStaticArray(fNrawch);
a4622d0f 2258 R__b.ReadStaticArray(fNrechits1D);
2259 R__b.ReadStaticArray(fNrechits3D);
2e5f0f7b 2260//
ddae0931 2261 R__b >> fChambers;
2262// Stream chamber related information
237c933d 2263 for (Int_t i =0; i<kNCH; i++) {
2e5f0f7b 2264 iChamber=(AliRICHChamber*) (*fChambers)[i];
ddae0931 2265 iChamber->Streamer(R__b);
2e5f0f7b 2266 segmentation=iChamber->GetSegmentationModel();
2267 segmentation->Streamer(R__b);
2268 response=iChamber->GetResponseModel();
ddae0931 2269 response->Streamer(R__b);
2e5f0f7b 2270 rawcladdress=(TClonesArray*) (*fRawClusters)[i];
2271 rawcladdress->Streamer(R__b);
a4622d0f 2272 rechitaddress1D=(TClonesArray*) (*fRecHits1D)[i];
2273 rechitaddress1D->Streamer(R__b);
2274 rechitaddress3D=(TClonesArray*) (*fRecHits3D)[i];
2275 rechitaddress3D->Streamer(R__b);
ddae0931 2276 digitsaddress=(TClonesArray*) (*fDchambers)[i];
2277 digitsaddress->Streamer(R__b);
2278 }
452a64c6 2279 R__b >> fDebugLevel;
2280 R__b >> fCkovNumber;
2281 R__b >> fCkovQuarz;
2282 R__b >> fCkovGap;
2283 R__b >> fCkovCsi;
2284 R__b >> fLostRfreo;
2285 R__b >> fLostRquar;
2286 R__b >> fLostAfreo;
2287 R__b >> fLostAquarz;
2288 R__b >> fLostAmeta;
2289 R__b >> fLostCsi;
2290 R__b >> fLostWires;
2291 R__b >> fFreonProd;
2292 R__b >> fMipx;
2293 R__b >> fMipy;
2294 R__b >> fFeedbacks;
2295 R__b >> fLostFresnel;
ddae0931 2296
2297 } else {
2298 R__b.WriteVersion(AliRICH::IsA());
2299 AliDetector::Streamer(R__b);
2e5f0f7b 2300 R__b << fNPadHits;
2301 R__b << fPadHits; // diff
2302 R__b << fNcerenkovs;
2303 R__b << fCerenkovs; // diff
ddae0931 2304 R__b << fDchambers;
2e5f0f7b 2305 R__b << fRawClusters;
a4622d0f 2306 R__b << fRecHits1D; //diff
2307 R__b << fRecHits3D; //diff
33984590 2308 R__b << fDebugLevel; //diff
237c933d 2309 R__b.WriteArray(fNdch, kNCH);
2310 R__b.WriteArray(fNrawch, kNCH);
a4622d0f 2311 R__b.WriteArray(fNrechits1D, kNCH);
2312 R__b.WriteArray(fNrechits3D, kNCH);
2313//
ddae0931 2314 R__b << fChambers;
2315// Stream chamber related information
237c933d 2316 for (Int_t i =0; i<kNCH; i++) {
2e5f0f7b 2317 iChamber=(AliRICHChamber*) (*fChambers)[i];
ddae0931 2318 iChamber->Streamer(R__b);
2e5f0f7b 2319 segmentation=iChamber->GetSegmentationModel();
2320 segmentation->Streamer(R__b);
2321 response=iChamber->GetResponseModel();
ddae0931 2322 response->Streamer(R__b);
2e5f0f7b 2323 rawcladdress=(TClonesArray*) (*fRawClusters)[i];
2324 rawcladdress->Streamer(R__b);
a4622d0f 2325 rechitaddress1D=(TClonesArray*) (*fRecHits1D)[i];
2326 rechitaddress1D->Streamer(R__b);
2327 rechitaddress3D=(TClonesArray*) (*fRecHits3D)[i];
2328 rechitaddress3D->Streamer(R__b);
ddae0931 2329 digitsaddress=(TClonesArray*) (*fDchambers)[i];
2330 digitsaddress->Streamer(R__b);
2331 }
452a64c6 2332 R__b << fDebugLevel;
2333 R__b << fCkovNumber;
2334 R__b << fCkovQuarz;
2335 R__b << fCkovGap;
2336 R__b << fCkovCsi;
2337 R__b << fLostRfreo;
2338 R__b << fLostRquar;
2339 R__b << fLostAfreo;
2340 R__b << fLostAquarz;
2341 R__b << fLostAmeta;
2342 R__b << fLostCsi;
2343 R__b << fLostWires;
2344 R__b << fFreonProd;
2345 R__b << fMipx;
2346 R__b << fMipy;
2347 R__b << fFeedbacks;
2348 R__b << fLostFresnel;
fe4da5cc 2349 }
ddae0931 2350}
2e5f0f7b 2351AliRICHPadHit* AliRICH::FirstPad(AliRICHHit* hit,TClonesArray *clusters )
ddae0931 2352{
2353//
2354 // Initialise the pad iterator
2355 // Return the address of the first padhit for hit
2356 TClonesArray *theClusters = clusters;
2357 Int_t nclust = theClusters->GetEntriesFast();
2358 if (nclust && hit->fPHlast > 0) {
2e5f0f7b 2359 sMaxIterPad=Int_t(hit->fPHlast);
2360 sCurIterPad=Int_t(hit->fPHfirst);
2361 return (AliRICHPadHit*) clusters->UncheckedAt(sCurIterPad-1);
ddae0931 2362 } else {
2363 return 0;
fe4da5cc 2364 }
ddae0931 2365
fe4da5cc 2366}
2367
2e5f0f7b 2368AliRICHPadHit* AliRICH::NextPad(TClonesArray *clusters)
fe4da5cc 2369{
237c933d 2370
2371 // Iterates over pads
2372
ddae0931 2373 sCurIterPad++;
2374 if (sCurIterPad <= sMaxIterPad) {
2e5f0f7b 2375 return (AliRICHPadHit*) clusters->UncheckedAt(sCurIterPad-1);
ddae0931 2376 } else {
2377 return 0;
2378 }
fe4da5cc 2379}
2380
fe4da5cc 2381
d28b5fc2 2382void AliRICH::Digitise(Int_t nev, Int_t flag, Option_t *option,Text_t *filename)
fe4da5cc 2383{
ddae0931 2384 // keep galice.root for signal and name differently the file for
2385 // background when add! otherwise the track info for signal will be lost !
2e5f0f7b 2386
c90dd3e2 2387 static Bool_t first=kTRUE;
237c933d 2388 static TFile *pFile;
2389 char *addBackground = strstr(option,"Add");
53a60579 2390 Int_t particle;
ddae0931 2391
2e5f0f7b 2392 FILE* points; //these will be the digits...
2393
2394 points=fopen("points.dat","w");
2395
2396 AliRICHChamber* iChamber;
a2f7eaf6 2397 AliSegmentation* segmentation;
ddae0931 2398
53a60579 2399 Int_t digitise=0;
ddae0931 2400 Int_t trk[50];
2401 Int_t chtrk[50];
2402 TObjArray *list=new TObjArray;
237c933d 2403 static TClonesArray *pAddress=0;
2404 if(!pAddress) pAddress=new TClonesArray("TVector",1000);
2e5f0f7b 2405 Int_t digits[5];
ddae0931 2406
237c933d 2407 AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
a2f7eaf6 2408 AliHitMap* pHitMap[10];
2e5f0f7b 2409 Int_t i;
237c933d 2410 for (i=0; i<10; i++) {pHitMap[i]=0;}
2411 if (addBackground ) {
ddae0931 2412 if(first) {
2413 fFileName=filename;
2414 cout<<"filename"<<fFileName<<endl;
237c933d 2415 pFile=new TFile(fFileName);
ddae0931 2416 cout<<"I have opened "<<fFileName<<" file "<<endl;
2e5f0f7b 2417 fHits2 = new TClonesArray("AliRICHHit",1000 );
2418 fClusters2 = new TClonesArray("AliRICHPadHit",10000);
c90dd3e2 2419 first=kFALSE;
ddae0931 2420 }
237c933d 2421 pFile->cd();
ddae0931 2422 // Get Hits Tree header from file
2423 if(fHits2) fHits2->Clear();
2424 if(fClusters2) fClusters2->Clear();
2e5f0f7b 2425 if(TrH1) delete TrH1;
2426 TrH1=0;
2427
ddae0931 2428 char treeName[20];
2429 sprintf(treeName,"TreeH%d",nev);
2e5f0f7b 2430 TrH1 = (TTree*)gDirectory->Get(treeName);
2431 if (!TrH1) {
ddae0931 2432 printf("ERROR: cannot find Hits Tree for event:%d\n",nev);
2433 }
2434 // Set branch addresses
2435 TBranch *branch;
2436 char branchname[20];
2437 sprintf(branchname,"%s",GetName());
2e5f0f7b 2438 if (TrH1 && fHits2) {
2439 branch = TrH1->GetBranch(branchname);
ddae0931 2440 if (branch) branch->SetAddress(&fHits2);
2441 }
2e5f0f7b 2442 if (TrH1 && fClusters2) {
2443 branch = TrH1->GetBranch("RICHCluster");
ddae0931 2444 if (branch) branch->SetAddress(&fClusters2);
2445 }
2446 }
33984590 2447
a2f7eaf6 2448 AliHitMap* hm;
2e5f0f7b 2449 Int_t countadr=0;
33984590 2450 Int_t counter=0;
2451 for (i =0; i<kNCH; i++) {
2452 iChamber=(AliRICHChamber*) (*fChambers)[i];
2453 segmentation=iChamber->GetSegmentationModel(1);
2454 pHitMap[i] = new AliRICHHitMapA1(segmentation, list);
2455 }
2456 //
2457 // Loop over tracks
2458 //
2459
2460 TTree *treeH = gAlice->TreeH();
2461 Int_t ntracks =(Int_t) treeH->GetEntries();
2462 for (Int_t track=0; track<ntracks; track++) {
2463 gAlice->ResetHits();
2464 treeH->GetEvent(track);
2465 //
2466 // Loop over hits
2467 for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
2468 mHit;
2469 mHit=(AliRICHHit*)pRICH->NextHit())
2470 {
2471
33984590 2472 Int_t nch = mHit->fChamber-1; // chamber number
53a60579 2473 Int_t index = mHit->Track();
33984590 2474 if (nch >kNCH) continue;
2475 iChamber = &(pRICH->Chamber(nch));
2476
53a60579 2477 TParticle *current = (TParticle*)(*gAlice->Particles())[index];
33984590 2478
53a60579 2479 if (current->GetPdgCode() >= 50000050)
2480 {
2481 TParticle *motherofcurrent = (TParticle*)(*gAlice->Particles())[current->GetFirstMother()];
2482 particle = motherofcurrent->GetPdgCode();
2483 }
2484 else
2485 {
2486 particle = current->GetPdgCode();
2487 }
2488
33984590 2489
2490 //printf("Flag:%d\n",flag);
2491 //printf("Track:%d\n",track);
2492 //printf("Particle:%d\n",particle);
2493
53a60579 2494 digitise=1;
33984590 2495
2496 if (flag == 1)
53a60579 2497 if(TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
2498 digitise=0;
33984590 2499
2500 if (flag == 2)
2501 if(TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
2502 || TMath::Abs(particle)==311)
53a60579 2503 digitise=0;
33984590 2504
2505 if (flag == 3 && TMath::Abs(particle)==2212)
53a60579 2506 digitise=0;
33984590 2507
2508 if (flag == 4 && TMath::Abs(particle)==13)
53a60579 2509 digitise=0;
33984590 2510
2511 if (flag == 5 && TMath::Abs(particle)==11)
53a60579 2512 digitise=0;
33984590 2513
2514 if (flag == 6 && TMath::Abs(particle)==2112)
53a60579 2515 digitise=0;
33984590 2516
2517
53a60579 2518 //printf ("Particle: %d, Flag: %d, Digitise: %d\n",particle,flag,digitise);
33984590 2519
2520
53a60579 2521 if (digitise)
ddae0931 2522 {
d28b5fc2 2523
33984590 2524 //
2525 // Loop over pad hits
2526 for (AliRICHPadHit* mPad=
2527 (AliRICHPadHit*)pRICH->FirstPad(mHit,fPadHits);
2528 mPad;
2529 mPad=(AliRICHPadHit*)pRICH->NextPad(fPadHits))
d28b5fc2 2530 {
33984590 2531 Int_t ipx = mPad->fPadX; // pad number on X
2532 Int_t ipy = mPad->fPadY; // pad number on Y
2533 Int_t iqpad = mPad->fQpad; // charge per pad
2534 //
2535 //
2536 //printf("X:%d, Y:%d, Q:%d\n",ipx,ipy,iqpad);
2537
a2f7eaf6 2538 Float_t thex, they, thez;
176a9917 2539 segmentation=iChamber->GetSegmentationModel(0);
a2f7eaf6 2540 segmentation->GetPadC(ipx,ipy,thex,they,thez);
33984590 2541 new((*pAddress)[countadr++]) TVector(2);
2542 TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]);
2543 trinfo(0)=(Float_t)track;
2544 trinfo(1)=(Float_t)iqpad;
2545
2546 digits[0]=ipx;
2547 digits[1]=ipy;
2548 digits[2]=iqpad;
2549
2550 AliRICHTransientDigit* pdigit;
2551 // build the list of fired pads and update the info
2552 if (!pHitMap[nch]->TestHit(ipx, ipy)) {
2553 list->AddAtAndExpand(new AliRICHTransientDigit(nch,digits),counter);
2554 pHitMap[nch]->SetHit(ipx, ipy, counter);
2555 counter++;
2556 pdigit=(AliRICHTransientDigit*)list->At(list->GetLast());
2557 // list of tracks
2558 TObjArray *trlist=(TObjArray*)pdigit->TrackList();
2559 trlist->Add(&trinfo);
2560 } else {
2561 pdigit=(AliRICHTransientDigit*) pHitMap[nch]->GetHit(ipx, ipy);
2562 // update charge
2563 (*pdigit).fSignal+=iqpad;
2564 // update list of tracks
2565 TObjArray* trlist=(TObjArray*)pdigit->TrackList();
2566 Int_t lastEntry=trlist->GetLast();
2567 TVector *ptrkP=(TVector*)trlist->At(lastEntry);
2568 TVector &ptrk=*ptrkP;
2569 Int_t lastTrack=Int_t(ptrk(0));
2570 Int_t lastCharge=Int_t(ptrk(1));
2571 if (lastTrack==track) {
2572 lastCharge+=iqpad;
2573 trlist->RemoveAt(lastEntry);
2574 trinfo(0)=lastTrack;
2575 trinfo(1)=lastCharge;
2576 trlist->AddAt(&trinfo,lastEntry);
2577 } else {
2578 trlist->Add(&trinfo);
2579 }
2580 // check the track list
2581 Int_t nptracks=trlist->GetEntriesFast();
2582 if (nptracks > 2) {
2583 printf("Attention - tracks: %d (>2)\n",nptracks);
2584 //printf("cat,nch,ix,iy %d %d %d %d \n",icat+1,nch,ipx,ipy);
2585 for (Int_t tr=0;tr<nptracks;tr++) {
2586 TVector *pptrkP=(TVector*)trlist->At(tr);
2587 TVector &pptrk=*pptrkP;
2588 trk[tr]=Int_t(pptrk(0));
2589 chtrk[tr]=Int_t(pptrk(1));
2590 }
2591 } // end if nptracks
2592 } // end if pdigit
2593 } //end loop over clusters
2594 }// track type condition
2595 } // hit loop
2596 } // track loop
2597
2598 // open the file with background
2599
2600 if (addBackground ) {
2601 ntracks =(Int_t)TrH1->GetEntries();
2602 //printf("background - icat,ntracks1 %d %d\n",icat,ntracks);
2603 //printf("background - Start loop over tracks \n");
2604 //
2605 // Loop over tracks
2606 //
2607 for (Int_t trak=0; trak<ntracks; trak++) {
2608 if (fHits2) fHits2->Clear();
2609 if (fClusters2) fClusters2->Clear();
2610 TrH1->GetEvent(trak);
2611 //
2612 // Loop over hits
2613 AliRICHHit* mHit;
2614 for(int j=0;j<fHits2->GetEntriesFast();++j)
2615 {
2616 mHit=(AliRICHHit*) (*fHits2)[j];
2617 Int_t nch = mHit->fChamber-1; // chamber number
2618 if (nch >6) continue;
2619 iChamber = &(pRICH->Chamber(nch));
2620 Int_t rmin = (Int_t)iChamber->RInner();
2621 Int_t rmax = (Int_t)iChamber->ROuter();
2622 //
2623 // Loop over pad hits
2624 for (AliRICHPadHit* mPad=
2625 (AliRICHPadHit*)pRICH->FirstPad(mHit,fClusters2);
2626 mPad;
2627 mPad=(AliRICHPadHit*)pRICH->NextPad(fClusters2))
2628 {
33984590 2629 Int_t ipx = mPad->fPadX; // pad number on X
2630 Int_t ipy = mPad->fPadY; // pad number on Y
2631 Int_t iqpad = mPad->fQpad; // charge per pad
2e5f0f7b 2632
a2f7eaf6 2633 Float_t thex, they, thez;
176a9917 2634 segmentation=iChamber->GetSegmentationModel(0);
a2f7eaf6 2635 segmentation->GetPadC(ipx,ipy,thex,they,thez);
33984590 2636 Float_t rpad=TMath::Sqrt(thex*thex+they*they);
2637 if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
2638 new((*pAddress)[countadr++]) TVector(2);
2639 TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]);
2640 trinfo(0)=-1; // tag background
2641 trinfo(1)=-1;
2642 digits[0]=ipx;
2643 digits[1]=ipy;
2644 digits[2]=iqpad;
2645 if (trak <4 && nch==0)
2646 printf("bgr - pHitMap[nch]->TestHit(ipx, ipy),trak %d %d\n",
2647 pHitMap[nch]->TestHit(ipx, ipy),trak);
2648 AliRICHTransientDigit* pdigit;
2649 // build the list of fired pads and update the info
2650 if (!pHitMap[nch]->TestHit(ipx, ipy)) {
2651 list->AddAtAndExpand(new AliRICHTransientDigit(nch,digits),counter);
2652
2653 pHitMap[nch]->SetHit(ipx, ipy, counter);
2654 counter++;
2655 printf("bgr new elem in list - counter %d\n",counter);
2656
2657 pdigit=(AliRICHTransientDigit*)list->At(list->GetLast());
2658 // list of tracks
2659 TObjArray *trlist=(TObjArray*)pdigit->TrackList();
2660 trlist->Add(&trinfo);
2661 } else {
2662 pdigit=(AliRICHTransientDigit*) pHitMap[nch]->GetHit(ipx, ipy);
2663 // update charge
2664 (*pdigit).fSignal+=iqpad;
2665 // update list of tracks
2666 TObjArray* trlist=(TObjArray*)pdigit->TrackList();
2667 Int_t lastEntry=trlist->GetLast();
2668 TVector *ptrkP=(TVector*)trlist->At(lastEntry);
2669 TVector &ptrk=*ptrkP;
2670 Int_t lastTrack=Int_t(ptrk(0));
2671 if (lastTrack==-1) {
2672 continue;
2673 } else {
2674 trlist->Add(&trinfo);
2675 }
2676 // check the track list
2677 Int_t nptracks=trlist->GetEntriesFast();
2678 if (nptracks > 0) {
2679 for (Int_t tr=0;tr<nptracks;tr++) {
2680 TVector *pptrkP=(TVector*)trlist->At(tr);
2681 TVector &pptrk=*pptrkP;
2682 trk[tr]=Int_t(pptrk(0));
2683 chtrk[tr]=Int_t(pptrk(1));
2684 }
2685 } // end if nptracks
2686 } // end if pdigit
2687 } //end loop over clusters
2688 } // hit loop
2689 } // track loop
ddae0931 2690 TTree *fAli=gAlice->TreeK();
237c933d 2691 if (fAli) pFile =fAli->GetCurrentFile();
2692 pFile->cd();
33984590 2693 } // if Add
2694
2695 Int_t tracks[10];
2696 Int_t charges[10];
2697 //cout<<"Start filling digits \n "<<endl;
2698 Int_t nentries=list->GetEntriesFast();
2699 //printf(" \n \n nentries %d \n",nentries);
2700
2701 // start filling the digits
2702
2703 for (Int_t nent=0;nent<nentries;nent++) {
2704 AliRICHTransientDigit *address=(AliRICHTransientDigit*)list->At(nent);
2705 if (address==0) continue;
2706
2707 Int_t ich=address->fChamber;
2708 Int_t q=address->fSignal;
2709 iChamber=(AliRICHChamber*) (*fChambers)[ich];
2710 AliRICHResponse * response=iChamber->GetResponseModel();
2711 Int_t adcmax= (Int_t) response->MaxAdc();
2712
2713
2714 // add white noise and do zero-suppression and signal truncation (new electronics,old electronics gaus 1.2,0.2)
2715 //printf("Treshold: %d\n",iChamber->fTresh->GetHitIndex(address->fPadX,address->fPadY));
2716 Int_t pedestal = iChamber->fTresh->GetHitIndex(address->fPadX,address->fPadY);
ddae0931 2717
33984590 2718 //printf("Pedestal:%d\n",pedestal);
2719 //Int_t pedestal=0;
9dda3582 2720 Float_t treshold = (pedestal + 4*2.2);
33984590 2721
9dda3582 2722 Float_t meanNoise = gRandom->Gaus(2.2, 0.3);
33984590 2723 Float_t noise = gRandom->Gaus(0, meanNoise);
2724 q+=(Int_t)(noise + pedestal);
2725 //q+=(Int_t)(noise);
2726 // magic number to be parametrised !!!
2727 if ( q <= treshold)
2728 {
2729 q = q - pedestal;
2730 continue;
2e5f0f7b 2731 }
33984590 2732 q = q - pedestal;
2733 if ( q >= adcmax) q=adcmax;
2734 digits[0]=address->fPadX;
2735 digits[1]=address->fPadY;
2736 digits[2]=q;
2737
2738 TObjArray* trlist=(TObjArray*)address->TrackList();
2739 Int_t nptracks=trlist->GetEntriesFast();
2740
2741 // this was changed to accomodate the real number of tracks
2742 if (nptracks > 10) {
2743 cout<<"Attention - tracks > 10 "<<nptracks<<endl;
2744 nptracks=10;
2745 }
2746 if (nptracks > 2) {
2747 printf("Attention - tracks > 2 %d \n",nptracks);
2748 //printf("cat,ich,ix,iy,q %d %d %d %d %d \n",
2749 //icat,ich,digits[0],digits[1],q);
2750 }
2751 for (Int_t tr=0;tr<nptracks;tr++) {
2752 TVector *ppP=(TVector*)trlist->At(tr);
2753 TVector &pp =*ppP;
2754 tracks[tr]=Int_t(pp(0));
2755 charges[tr]=Int_t(pp(1));
2756 } //end loop over list of tracks for one pad
2757 if (nptracks < 10 ) {
2758 for (Int_t t=nptracks; t<10; t++) {
2759 tracks[t]=0;
2760 charges[t]=0;
ddae0931 2761 }
33984590 2762 }
2763 //write file
2764 if (ich==2)
2765 fprintf(points,"%4d, %4d, %4d\n",digits[0],digits[1],digits[2]);
2766
2767 // fill digits
2768 pRICH->AddDigits(ich,tracks,charges,digits);
2769 }
2770 gAlice->TreeD()->Fill();
2771
2772 list->Delete();
2773 for(Int_t ii=0;ii<kNCH;++ii) {
2774 if (pHitMap[ii]) {
2775 hm=pHitMap[ii];
2776 delete hm;
2777 pHitMap[ii]=0;
2778 }
2779 }
2780
2781 //TTree *TD=gAlice->TreeD();
2782 //Stat_t ndig=TD->GetEntries();
2783 //cout<<"number of digits "<<ndig<<endl;
2784 TClonesArray *fDch;
2785 for (int k=0;k<kNCH;k++) {
2786 fDch= pRICH->DigitsAddress(k);
2787 int ndigit=fDch->GetEntriesFast();
2788 printf ("Chamber %d digits %d \n",k,ndigit);
2789 }
2790 pRICH->ResetDigits();
ddae0931 2791 char hname[30];
2792 sprintf(hname,"TreeD%d",nev);
33984590 2793 gAlice->TreeD()->Write(hname,kOverwrite,0);
2794
2795 // reset tree
2796 // gAlice->TreeD()->Reset();
2e5f0f7b 2797 delete list;
237c933d 2798 pAddress->Clear();
33984590 2799 // gObjectTable->Print();
2e5f0f7b 2800}
ddae0931 2801
237c933d 2802AliRICH& AliRICH::operator=(const AliRICH& rhs)
2803{
2804// Assignment operator
2805 return *this;
2806
2807}
ddae0931 2808
237c933d 2809Int_t AliRICH::MakePadHits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
2810{
2811//
2812// Calls the charge disintegration method of the current chamber and adds
2813// the simulated cluster to the root treee
2814//
176a9917 2815 Int_t clhits[5];
2816 Float_t newclust[4][500];
237c933d 2817 Int_t nnew;
2818
2819//
2820// Integrated pulse height on chamber
2821
2822 clhits[0]=fNhits+1;
2823
2824 ((AliRICHChamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust, res);
2825 Int_t ic=0;
2826
2827//
2828// Add new clusters
2829 for (Int_t i=0; i<nnew; i++) {
176a9917 2830 if (Int_t(newclust[0][i]) > 0) {
237c933d 2831 ic++;
237c933d 2832// Cluster Charge
176a9917 2833 clhits[1] = Int_t(newclust[0][i]);
237c933d 2834// Pad: ix
176a9917 2835 clhits[2] = Int_t(newclust[1][i]);
237c933d 2836// Pad: iy
176a9917 2837 clhits[3] = Int_t(newclust[2][i]);
237c933d 2838// Pad: chamber sector
176a9917 2839 clhits[4] = Int_t(newclust[3][i]);
237c933d 2840
2841 AddPadHit(clhits);
2842 }
2843 }
2844return nnew;
2845}
ddae0931 2846