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