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