ROOTSYS for AMORE
[u/mrichter/AliRoot.git] / FIT / FITsim / AliFITv2.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16
17 /////////////////////////////////////////////////////////////////////
18 //                                                                 //
19 // FIT detector full geometry  version 1                             //
20 //
21 //Begin Html       
22 /*
23 <img src="gif/AliFITv2Class.gif">
24 */
25 //End Html
26 //                                                                  //
27 //                                                                  //
28 //////////////////////////////////////////////////////////////////////
29
30 #include <Riostream.h>
31 #include <stdlib.h>
32
33 #include "TGeoCompositeShape.h"
34 #include "TGeoManager.h"
35 #include "TGeoMatrix.h"
36 #include "TGeoVolume.h"
37 #include "TGeoTube.h"
38 #include "TGeoBBox.h"
39 #include "TGeoNode.h"
40
41
42 #include <TGeoGlobalMagField.h>
43 #include <TGraph.h>
44 #include <TLorentzVector.h>
45 #include <TMath.h>
46 #include <TVirtualMC.h>
47 #include <TString.h>
48
49 #include "AliLog.h"
50 #include "AliMagF.h"
51 #include "AliRun.h"
52
53 #include "AliFITHits.h"
54 #include "AliFITv2.h"
55
56 #include "AliMC.h"
57 #include "AliCDBLocal.h"
58 #include "AliCDBStorage.h"
59 #include "AliCDBManager.h"
60 #include "AliCDBEntry.h"
61 #include "AliTrackReference.h"
62
63 ClassImp(AliFITv2)
64
65
66 //--------------------------------------------------------------------
67 AliFITv2::AliFITv2():  AliFIT(),
68                      fIdSens1(0),
69                      fPMTeff(0x0)
70
71 {
72   //
73   // Standart constructor for T0 Detector version 0
74 }
75 //--------------------------------------------------------------------
76 AliFITv2::AliFITv2(const char *name, const char *title):
77   AliFIT(name,title),
78   fIdSens1(0),
79   fPMTeff(0x0)
80
81 {
82   //
83   // Standart constructor for T0 Detector version 0
84   //
85   fIshunt = 2; 
86   SetPMTeff();
87 }
88 //_____________________________________________________________________________
89
90 AliFITv2::~AliFITv2() 
91 {
92   // desctructor  
93 }
94
95 //-------------------------------------------------------------------------
96 void AliFITv2::CreateGeometry()
97 {
98   //
99   // Create the geometry of FIT Detector version 1 full geometry
100   //
101   // begin Html
102   //
103
104   Int_t *idtmed = fIdtmed->GetArray();
105   Float_t zdetC = 85;
106   Float_t zdetA = 335;
107   
108   Int_t idrotm[999];
109   Double_t x,y,z;
110   Float_t pstartC[3] = {6., 20 ,5};
111   Float_t pstartA[3] = {2.55, 20 ,5};
112   // Float_t pinstart[3] = {3.2,3.2,3.9};
113   Float_t pinstart[3] = {3.2,3.2,4.};
114   Float_t  pmcp[3] = {3.19, 3.19, 2.8}; //MCP
115   Float_t ptop[3] = {1.324, 1.324, 1.};//cherenkov radiator
116   Float_t preg[3] = {1.324, 1.324, 0.05};//photcathode 
117
118   Float_t zV0A = 329.;
119   Float_t pV0Amother[3] = {4.25, 41.25, 0.6};
120   Float_t pV0A[3] = {4.3, 41.2, 0.5};
121
122   AliMatrix(idrotm[901], 90., 0., 90., 90., 180., 0.);
123   
124   //-------------------------------------------------------------------
125   //  T0 volume 
126   //-------------------------------------------------------------------
127   //C side
128
129    Float_t zc[36] = {2,0, 2,0, 2,0, 2,0, 2,2, 2,2, 2,2, 2,2, 0,0,0,0};
130    Float_t za[25] = {2,0, 2,0, 2,2, 0,2, 0,2, 2,0, 2,0, 0,0, 0,0,0,0};
131   
132    
133   TGeoVolumeAssembly * stlinA = new TGeoVolumeAssembly("0STL");  // A side mother 
134   TGeoVolumeAssembly * stlinC = new TGeoVolumeAssembly("0STR");  // C side mother 
135  //T0 interior
136   TVirtualMC::GetMC()->Gsvolu("0INS","BOX",idtmed[kAir],pinstart,3);
137   TGeoVolume *ins = gGeoManager->GetVolume("0INS");
138  // 
139   TGeoTranslation *tr[40];
140   TString nameTr;
141   //C side
142   Float_t xa=-12.8;
143   Float_t ya=-12.8;
144   Int_t itr=0;
145   Int_t itrHole=0;
146   for (Int_t itrx=0; itrx<5; itrx++) {
147     for (Int_t itry=0; itry<5; itry++) {
148       nameTr = Form("0TR%i",itr+1);
149       z=-pstartA[2]+pinstart[2]/*+za[itr]*/;
150       if(itr !=12){
151         if(TMath::Abs(xa)<10 && TMath::Abs(ya)<10) z= z-2;
152         tr[itr] = new TGeoTranslation(nameTr.Data(),xa,ya, z );
153         tr[itr]->RegisterYourself();
154         stlinA->AddNode(ins,itr,tr[itr]);
155       }
156       printf(" A %i %f %f %f \n",itr, xa, ya, z+zdetA);
157       itr++;
158       ya+=6.4;
159     }
160     ya=-12.8;
161     xa+=6.4;
162   }
163   Float_t xc=-16;
164   Float_t yc=-16;
165   for (Int_t itrx=0; itrx<6; itrx++) {
166     for (Int_t itry=0; itry<6; itry++) {
167       nameTr = Form("0TR%i",itr+1);
168       z=-pstartC[2]+pinstart[2]/*+zc[itr]*/;
169       if (itr!=39 && itr!=40 && itr!=45 &&itr!=46) {
170         if( TMath::Abs(xc)<10 &&  TMath::Abs(yc)<10) z= z+2;
171         tr[itr] = new TGeoTranslation(nameTr.Data(),xc,yc, z );
172         tr[itr]->RegisterYourself();
173         stlinC->AddNode(ins,itr,tr[itr]);
174       }
175       printf(" C %i %f %f %f \n",itr,xc, yc, z+zdetC);
176       itr++;
177       yc+=6.4;
178     }
179     yc=-16;
180     xc+=6.4;
181   }
182     TGeoVolume *alice = gGeoManager->GetVolume("ALIC");
183   alice->AddNode(stlinA,1,new TGeoTranslation(0,0, zdetA ) );
184   //  alice->AddNode(stlinC,1,new TGeoTranslation(0,0, zdetC ) );
185   TGeoRotation * rotC = new TGeoRotation( "rotC",90., 0., 90., 90., 180., 0.);
186   alice->AddNode(stlinC,1, new TGeoCombiTrans(0., 0., -zdetC , rotC) );
187   
188   x=0;
189   y=0;
190    
191   // Entry window (glass)
192   TVirtualMC::GetMC()->Gsvolu("0TOP","BOX",idtmed[kOpGlass],ptop,3); //glass
193   TGeoVolume *top = gGeoManager->GetVolume("0TOP");
194   TVirtualMC::GetMC()->Gsvolu ("0REG", "BOX", idtmed[kOpGlassCathode], preg, 3); 
195   TGeoVolume *cat = gGeoManager->GetVolume("0REG");
196   TVirtualMC::GetMC()->Gsvolu("0MCP","BOX",idtmed[kGlass],pmcp,3); //glass
197   TGeoVolume *mcp = gGeoManager->GetVolume("0MCP");
198  
199   Int_t ntops=0;
200    Float_t xin=0, yin=0;
201    for (Int_t ix=0; ix<2; ix++) {
202      xin = - pinstart[0] + 0.55 + (ix+0.5)*2*ptop[0] ;
203      for (Int_t iy=0; iy<2 ; iy++) {
204        z = - pinstart[2]+ptop[2];
205        yin = - pinstart[1] + 0.55 + (iy+0.5)*2*ptop[1];
206        ntops++;
207        ins->AddNode(top, ntops, new TGeoTranslation(xin,yin,z) );
208        // printf(" 0TOP  full x %f y %f z %f \n", xin, yin, z);
209        z = -pinstart[2] + 2 * ptop[2] + preg[2];
210        ins->AddNode(cat, ntops, new TGeoTranslation(xin,yin,z) );
211
212        // printf(" GEOGEO  %i %i %i %f %f %f %f %f %f \n", ntops, ix, iy,
213        //  xin,yin,x1[ntops],y1[ntops],x1[ntops]+xin,y1[ntops]+yin);
214      }
215    }
216 // MCP
217    z=-pinstart[2] + 2*ptop[2] + 2*preg[2] + pmcp[2];
218   ins->AddNode(mcp, 1 , new TGeoTranslation(0,0,z) );
219 /*
220   //V0A 
221    TVirtualMC::GetMC()->Gsvolu("0V0AM","TUBE",idtmed[kAir],pV0Amother,3);
222    TVirtualMC::GetMC()->Gspos ("0V0AM",1, "ALIC", 0,0,zV0A , 0, "ONLY");
223    TVirtualMC::GetMC()->Gsvolu("0V0A","TUBE",idtmed[kSensAir],pV0A,3);
224    TVirtualMC::GetMC()->Gspos ("0V0A",1, "0V0AM", 0, 0, 0, 0, "ONLY");
225 */
226
227  
228 }    
229 //------------------------------------------------------------------------
230 void AliFITv2::AddAlignableVolumes() const
231 {
232   //
233   // Create entries for alignable volumes associating the symbolic volume
234   // name with the corresponding volume path. Needs to be syncronized with
235   // eventual changes in the geometry.
236   //
237   TString volPath;
238   TString symName, sn;
239   TString vpAalign = "/ALIC_1/0STL_1";
240   TString vpCalign = "/ALIC_1/0STR_1";
241   for (Int_t imod=0; imod<2; imod++)  {
242     if (imod==0) {volPath  = vpCalign; symName="/ALIC_1/0STL"; }
243     if (imod==1) {volPath  = vpAalign; symName="/ALIC_1/0STR"; }
244     
245     AliDebug(2,"--------------------------------------------");
246     AliDebug(2,Form("volPath=%s\n",volPath.Data()));
247     AliDebug(2,Form("symName=%s\n",symName.Data()));
248     AliDebug(2,"--------------------------------------------");
249     if(!gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data()))
250       AliFatal(Form("Alignable entry %s not created. Volume path %s not valid",
251 symName.Data(),volPath.Data()));
252    }
253 }   
254 //------------------------------------------------------------------------
255 void AliFITv2::CreateMaterials()
256 {
257
258    Int_t isxfld   = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Integ();
259    Float_t sxmgmx = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Max();
260    //   Float_t a,z,d,radl,absl,buf[1];
261    // Int_t nbuf;
262 // AIR
263                                                                                 
264    Float_t aAir[4]={12.0107,14.0067,15.9994,39.948};
265    Float_t zAir[4]={6.,7.,8.,18.};
266    Float_t wAir[4]={0.000124,0.755267,0.231781,0.012827};
267    Float_t dAir = 1.20479E-3;
268    Float_t dAir1 = 1.20479E-11;
269  // Radiator  glass SiO2
270    Float_t aglass[2]={28.0855,15.9994};
271    Float_t zglass[2]={14.,8.};
272    Float_t wglass[2]={1.,2.};
273    Float_t dglass=2.65;
274  // MCP glass SiO2
275    Float_t dglass_mcp=1.3;
276 //*** Definition Of avaible T0 materials ***
277    AliMixture(1, "Vacuum$", aAir, zAir, dAir1,4,wAir);
278    AliMixture(2, "Air$", aAir, zAir, dAir,4,wAir);
279    AliMixture( 4, "MCP glass   $",aglass,zglass,dglass_mcp,-2,wglass);
280    AliMixture( 24, "Radiator Optical glass$",aglass,zglass,dglass,-2,wglass);
281    
282    AliMedium(1, "Air$", 2, 0, isxfld, sxmgmx, 10., .1, 1., .003, .003);
283    AliMedium(3, "Vacuum$", 1, 0, isxfld, sxmgmx, 10., .01, .1, .003, .003);
284    AliMedium(6, "Glass$", 4, 0, isxfld, sxmgmx, 10., .01, .1, .003, .003);
285   
286    AliMedium(16, "OpticalGlass$", 24, 1, isxfld, sxmgmx, 10., .01, .1, .003, .003);
287     AliMedium(19, "OpticalGlassCathode$", 24, 1, isxfld, sxmgmx, 10., .01, .1, .003, .003);
288    AliMedium(22, "SensAir$", 2, 1, isxfld, sxmgmx, 10., .1, 1., .003, .003);
289    
290    AliDebugClass(1,": ++++++++++++++Medium set++++++++++");
291    
292    
293 }
294
295 //-------------------------------------------------------------------
296 void AliFITv2::DefineOpticalProperties()
297 {
298
299
300 // Optical properties definition.
301    Int_t *idtmed = fIdtmed->GetArray();
302 // Definition Cherenkov parameters
303    int i;
304    const Int_t kNbins=31;
305    
306    Float_t rindexSiO2[kNbins],  efficAll[kNbins], rindexAir[kNbins], absorAir[kNbins],rindexCathodeNext[kNbins], absorbCathodeNext[kNbins];
307    Double_t efficMet[kNbins], aReflMet[kNbins];
308            
309    // quartz 20mm
310    Float_t aAbsSiO2[kNbins]={29.0, 28.6, 28.3, 27.7, 27.3, 26.7, 26.4, 
311                              25.9, 25.3, 24.9, 24.5, 23.7, 
312                              23.2, 22.8, 22.4, 21.8, 21.3,
313                              22.8, 22.1, 21.7, 21.2, 20.5, 
314                              19.9, 19.3, 18.7, 18.0, 17.1, 
315                              16.3, 15.3, 14.3, 14.3   };
316           
317    Float_t aPckov[kNbins]  ={3.87, 3.94, 4.02, 4.11, 4.19, 4.29, 4.38,
318                              4.48, 4.58, 4.69, 4.81, 4.93, 
319                              5.05, 5.19, 5.33, 5.48, 5.63,
320                              5.8,  5.97, 6.16, 6.36, 6.57, 
321                              6.8,  7.04, 7.3,  7.58, 7.89, 
322                              8.22, 8.57, 8.97, 9.39 };  
323   Double_t dPckov[kNbins]  ={3.87, 3.94, 4.02, 4.11, 4.19, 4.29, 4.38,
324                              4.48, 4.58, 4.69, 4.81, 4.93,
325                              5.05, 5.19, 5.33, 5.48, 5.63,
326                              5.8,  5.97, 6.16, 6.36, 6.57,
327                              6.8,  7.04, 7.3,  7.58, 7.89,
328                              8.22, 8.57, 8.97, 9.39 };
329
330    /*     
331    Float_t effCathode[kNbins]={0.11, 0.13, 0.15, 0.16, 0.18, 0.19, 0.20,
332                               0.21, 0.22, 0.23, 0.24, 0.26, 
333                               0.27, 0.29, 0.30, 0.29, 0.29, 
334                               0.28, 0.28, 0.27, 0.26, 0.25, 
335                               0.25, 0.23, 0.20, 0.19, 0.17,
336                               0.17, 0.17, 0.2, 0.23};
337    */     
338    //  Float_t aAbsSiO2[kNbins]; //quartz 30mm
339  for(i=0;i<kNbins;i++)
340
341     {
342       aPckov[i]=aPckov[i]*1e-9;//Photons energy bins 4 eV - 8.5 eV step 0.1 eV
343       dPckov[i]=dPckov[i]*1e-9;//Photons energy bins 4 eV - 8.5 eV step 0.1 eV 
344       //      rindexAir[i]=0.0001;
345       rindexAir[i] = 1.;
346       rindexSiO2[i]=1.458; //refractive index for qwarts
347       rindexCathodeNext[i]=0;
348       efficAll[i]=1.;
349       efficMet[i]=0.;
350       aReflMet[i]=1.;
351       //      aAbsSiO2[i]=28.5; //quartz 30mm
352       absorAir[i]=0.3;      
353       absorbCathodeNext[i]=0;
354
355     }
356   
357   TVirtualMC::GetMC()->SetCerenkov (idtmed[kOpGlass], kNbins, aPckov, aAbsSiO2, efficAll, rindexSiO2 );
358    // TVirtualMC::GetMC()->SetCerenkov (idtmed[kOpGlassCathode], kNbins, aPckov, aAbsSiO2, effCathode, rindexSiO2 );
359    TVirtualMC::GetMC()->SetCerenkov (idtmed[kOpGlassCathode], kNbins, aPckov, aAbsSiO2,efficAll , rindexSiO2 );
360    //  TVirtualMC::GetMC()->SetCerenkov (idtmed[kOpAir], kNbins, aPckov,absorAir , efficAll,rindexAir );
361    //   TVirtualMC::GetMC()->SetCerenkov (idtmed[kOpAirNext], kNbins, aPckov,absorbCathodeNext , efficAll, rindexCathodeNext);
362
363    //Define a boarder for radiator optical properties
364    TVirtualMC::GetMC()->DefineOpSurface("surfRd", kUnified /*kGlisur*/,kDielectric_metal,kPolished, 0.);
365    TVirtualMC::GetMC()->SetMaterialProperty("surfRd", "EFFICIENCY", kNbins, dPckov, efficMet);
366    TVirtualMC::GetMC()->SetMaterialProperty("surfRd", "REFLECTIVITY", kNbins, dPckov, aReflMet);
367
368
369 }
370
371 //-------------------------------------------------------------------
372 void AliFITv2::Init()
373 {
374 // Initialises version 0 of the Forward Multiplicity Detector
375 //
376   AliFIT::Init();
377   fIdSens1=TVirtualMC::GetMC()->VolId("0REG");
378   fIdSens2=TVirtualMC::GetMC()->VolId("0V0A");
379
380    AliDebug(1,Form("%s: *** FIT version 1 initialized ***\n",ClassName()));
381 }
382
383 //-------------------------------------------------------------------
384
385 void AliFITv2::StepManager()
386 {
387   //
388   // Called for every step in the T0 Detector
389   //
390   Int_t id,copy,copy1;
391   static Float_t hits[6];
392   static Int_t vol[3];
393   TLorentzVector pos;
394   TLorentzVector mom;
395   
396   //   TClonesArray &lhits = *fHits;
397   
398   if(!TVirtualMC::GetMC()->IsTrackAlive()) return; // particle has disappeared
399   
400   id=TVirtualMC::GetMC()->CurrentVolID(copy);  
401   // Check the sensetive volume
402   if(id==fIdSens1 ) { 
403     if(TVirtualMC::GetMC()->IsTrackEntering()) {
404       TVirtualMC::GetMC()->CurrentVolOffID(1,copy1);
405       vol[1] = copy1;
406       vol[0]=copy;
407       TVirtualMC::GetMC()->TrackPosition(pos);
408       hits[0] = pos[0];
409       hits[1] = pos[1];
410       hits[2] = pos[2];
411       if(pos[2]<0) vol[2] = 0;
412       else vol[2] = 1 ;
413       //      printf(" volumes pmt %i mcp %i side %i x %f y %f z %f\n",  vol[0], vol[1],  vol[2], hits[0], hits[1], hits[2] );
414       
415       Float_t etot=TVirtualMC::GetMC()->Etot();
416       hits[3]=etot;
417       Int_t iPart= TVirtualMC::GetMC()->TrackPid();
418       Int_t partID=TVirtualMC::GetMC()->IdFromPDG(iPart);
419       hits[4]=partID;
420       Float_t ttime=TVirtualMC::GetMC()->TrackTime();
421       hits[5]=ttime*1e12;
422       if (TVirtualMC::GetMC()->TrackPid() == 50000050)   // If particles is photon then ...
423         {
424           //      if(RegisterPhotoE(vol[1]-1,hits[3])) {
425           if(RegisterPhotoE(hits[3])) {
426             AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),vol,hits);
427             // Create a track reference at the exit of photocatode
428           }         
429         }
430       
431       //charge particle 
432       if ( TVirtualMC::GetMC()->TrackCharge() )
433         AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kFIT);
434       
435      }// trck entering          
436   } //sensitive
437   //V0A
438   if(id==fIdSens2 ) { 
439     if ( TVirtualMC::GetMC()->TrackCharge()  ) {
440       if(TVirtualMC::GetMC()->IsTrackEntering()) {
441         TVirtualMC::GetMC()->TrackPosition(pos);
442         hits[0] = pos[0];
443         hits[1] = pos[1];
444         hits[2] = pos[2];
445         vol[0]=0;
446         vol[1]=0;
447         vol[2]=2;
448         
449         Float_t etot=TVirtualMC::GetMC()->Etot();
450         hits[3]=etot;
451         Int_t iPart= TVirtualMC::GetMC()->TrackPid();
452         Int_t partID=TVirtualMC::GetMC()->IdFromPDG(iPart);
453         hits[4]=partID;
454         Float_t ttime=TVirtualMC::GetMC()->TrackTime();
455         hits[5]=ttime*1e12;
456         AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),vol,hits);
457         //      printf(" volumes pmt %i mcp %i vol %i x %f y %f z %f particle %i all \n",  vol[0], vol[1],  vol[2], hits[0], hits[1], hits[2], hits[4]);
458       }
459     }    
460   }
461
462 }
463
464
465 //------------------------------------------------------------------------
466 Bool_t AliFITv2::RegisterPhotoE(Double_t energy)
467 {
468   
469   
470   //  Float_t hc=197.326960*1.e6; //mev*nm
471   Double_t hc=1.973*1.e-6; //gev*nm
472   Float_t lambda=hc/energy;
473   Float_t eff = fPMTeff->Eval(lambda);
474   Double_t  p = gRandom->Rndm();
475   
476   if (p > eff)
477     return kFALSE;
478   
479   return kTRUE;
480 }
481
482 //----------------------------------------------------------------------------
483
484 void AliFITv2::SetPMTeff()
485 {
486   Float_t lambda[50];
487   Float_t eff[50 ] = {0,        0,       0.23619,  0.202909, 0.177913, 
488                     0.175667, 0.17856, 0.190769, 0.206667, 0.230286,
489                     0.252276, 0.256267,0.26,     0.27125,  0.281818,
490                     0.288118, 0.294057,0.296222, 0.301622, 0.290421, 
491                     0.276615, 0.2666,  0.248,    0.23619,  0.227814, 
492                     0.219818, 0.206667,0.194087, 0.184681, 0.167917, 
493                     0.154367, 0.1364,  0.109412, 0.0834615,0.0725283, 
494                     0.0642963,0.05861, 0.0465,   0.0413333,0.032069, 
495                     0.0252203,0.02066, 0.016262, 0.012,    0.00590476,
496                     0.003875, 0.00190, 0,        0,        0          } ;
497   for (Int_t i=0; i<50; i++) lambda[i]=200+10*i; 
498
499   fPMTeff = new TGraph(50,lambda,eff);
500  
501 }
502
503
504
505
506