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