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