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