Update to AD geometry:
[u/mrichter/AliRoot.git] / AD / ADsim / AliAD.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 /* $Id: AliAD.cxx  $ */
17
18 ///////////////////////////////////////////////////////////////////////////
19 //                                                                       //
20 //                  AD (ALICE Diffractive)  Detector                     //
21 //                                                                       //
22 //  This class contains the base procedures for the AD  detector         //
23 //  All comments should be sent to :                                     //
24 //                                                                       //
25 //                                                                       //
26 ///////////////////////////////////////////////////////////////////////////
27
28
29 // --- Standard libraries ---
30 #include <Riostream.h>
31 #include <stdlib.h>
32
33 // --- ROOT libraries ---
34 #include <TNamed.h>
35 #include "TROOT.h"
36 #include "TFile.h"
37 #include "TNetFile.h"
38 #include "TRandom.h"
39 #include "TTree.h"
40 #include "TBranch.h"
41 #include "TClonesArray.h"
42 #include "TGeoGlobalMagField.h"
43 #include "AliMagF.h"
44 #include "TStopwatch.h"
45 #include "TParameter.h"
46 #include "TF1.h"
47
48 // --- AliRoot header files ---
49 #include "AliRun.h"
50 #include "AliMC.h"
51 #include "AliAD.h"
52 #include "AliADhit.h"
53 #include "AliADLoader.h"
54 #include "AliADDigitizer.h"
55 #include "AliDigitizationInput.h"
56 #include "AliADdigit.h"
57 #include "AliADSDigit.h"
58 #include "AliADBuffer.h"
59 #include "AliDAQ.h"
60 #include "AliRawReader.h"
61 #include "AliCDBManager.h"
62 #include "AliCDBEntry.h"
63 #include "AliADReconstructor.h"
64
65 ClassImp(AliAD)
66  //__________________________________________________________________
67 AliAD::AliAD()
68    : AliDetector(),
69      fSetADATwoInstalled(0),
70      fSetADCTwoInstalled(0)
71 {
72         /// Default Constructor
73     
74
75 }
76
77 //_____________________________________________________________________________
78 AliAD::AliAD(const char *name, const char *title)
79    : AliDetector(name,title),
80      fSetADATwoInstalled(kTRUE),
81      fSetADCTwoInstalled(kTRUE)
82
83 {
84   
85    // Standard constructor for AD Detector
86   
87
88 }
89
90 //_____________________________________________________________________________
91 AliAD::~AliAD()
92 {
93    //
94    // Default destructor for AD Detector
95    //
96   
97 }
98
99 //_____________________________________________________________________________
100 void AliAD::CreateMaterials()
101 {
102   //
103   // MATERIALS FOR ADC AND ADA
104   //
105   
106   // Parameters for simulation scope for ADA and ADC (stolen from AliVZEROv7::CreateMaterials )
107   Int_t    fieldType       = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Integ();     // Field type 
108   Double_t maxField        = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Max();       // Field max.
109   // Int_t   isxfld1 = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Integ();
110   // Float_t sxmgmx  = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Max();
111   Int_t    isxfld2         = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->PrecInteg();
112   Double_t maxBending      = 10;    // Max Angle
113   Double_t maxStepSize     = 0.01;  // Max step size 
114   Double_t maxEnergyLoss   = 1;     // Max Delta E
115   Double_t precision       = 0.003; // Precision
116   Double_t minStepSize     = 0.003; // Minimum step size 
117   Float_t  density,  as[11], zs[11], ws[11];
118   Double_t radLength, absLength, a_ad, z_ad;
119   Int_t    id;
120   
121   //
122   // Air
123   //
124   Float_t aAir[4]={12.0107,14.0067,15.9994,39.948};
125   Float_t zAir[4]={6.,7.,8.,18.};
126   Float_t wAir[4]={0.000124,0.755267,0.231781,0.012827};
127   Float_t dAir = 1.20479E-3;
128   Float_t dAir1 = 1.20479E-10;
129   // Steel  
130   Float_t asteel[4] = { 55.847,51.9961, 58.6934, 28.0855 };
131   Float_t zsteel[4] = {    26.,    24.,     28., 14.     };
132   Float_t wsteel[4] = {   .715,    .18,      .1, .005    };
133   //
134   // Cast iron
135   //
136   Float_t acasti[4] = {55.847,12.011,28.085,54.938};
137   Float_t zcasti[4] = {26.,6.,14.,25.};
138   Float_t wcasti[4] = {0.929,0.035,0.031,0.005};
139   // --- Define the various materials for GEANT --- 
140   //     Aluminum 
141   AliMaterial(9,  "ALU1      ", 26.98, 13., 2.7, 8.9, 37.2);
142   AliMaterial(29, "ALU2      ", 26.98, 13., 2.7, 8.9, 37.2);
143   AliMaterial(49, "ALU3      ", 26.98, 13., 2.7, 8.9, 37.2);
144   
145   //     Iron 
146   AliMaterial(10, "IRON1     ", 55.85, 26., 7.87, 1.76, 17.1);
147   AliMaterial(30, "IRON2     ", 55.85, 26., 7.87, 1.76, 17.1);
148   AliMaterial(50, "IRON3     ", 55.85, 26., 7.87, 1.76, 17.1);
149
150   //
151   //     Copper
152   AliMaterial(11, "COPPER1   ", 63.55, 29., 8.96, 1.43, 15.1);
153   AliMaterial(31, "COPPER2   ", 63.55, 29., 8.96, 1.43, 15.1);
154   AliMaterial(51, "COPPER3   ", 63.55, 29., 8.96, 1.43, 15.1);
155   
156   //     Vacuum 
157   AliMixture(16, "VACUUM1 ", aAir, zAir, dAir1, 4, wAir);
158   AliMixture(36, "VACUUM2 ", aAir, zAir, dAir1, 4, wAir);
159   AliMixture(56, "VACUUM3 ", aAir, zAir, dAir1, 4, wAir);
160   
161   //     Stainless Steel 
162   AliMixture(19, "STAINLESS STEEL1", asteel, zsteel, 7.88, 4, wsteel);
163   AliMixture(39, "STAINLESS STEEL2", asteel, zsteel, 7.88, 4, wsteel);
164   AliMixture(59, "STAINLESS STEEL3", asteel, zsteel, 7.88, 4, wsteel);
165   
166   //     Cast iron
167   AliMixture(18, "CAST IRON1", acasti, zcasti, 7.2, 4, wcasti);
168   AliMixture(38, "CAST IRON2", acasti, zcasti, 7.2, 4, wcasti);
169   AliMixture(58, "CAST IRON3", acasti, zcasti, 7.2, 4, wcasti);
170   // **************** 
171   //     Defines tracking media parameters. 
172   //     Les valeurs sont commentees pour laisser le defaut 
173   //     a GEANT (version 3-21, page CONS200), f.m. 
174   Float_t epsil, stmin, tmaxfd, deemax, stemax;
175   epsil  = .001;  // Tracking precision, 
176   stemax = -1.;   // Maximum displacement for multiple scat 
177   tmaxfd = -20.;  // Maximum angle due to field deflection 
178   deemax = -.3;   // Maximum fractional energy loss, DLS 
179   stmin  = -.8;
180   // *************** 
181   
182   //    Aluminum 
183   AliMedium(9,  "ALU_C0          ", 9, 0,  fieldType, maxField, tmaxfd, stemax, deemax, epsil, stmin);
184   AliMedium(29, "ALU_C1          ", 29, 0, fieldType, maxField, tmaxfd, stemax, deemax, epsil, stmin);
185   AliMedium(49, "ALU_C2          ", 49, 0, fieldType, maxField, tmaxfd, stemax, deemax, epsil, stmin);
186   
187   //    Iron 
188   AliMedium(10, "FE_C0           ", 10, 0, fieldType, maxField, tmaxfd, stemax, deemax, epsil, stmin);
189   AliMedium(30, "FE_C1           ", 30, 0, fieldType, maxField, tmaxfd, stemax, deemax, epsil, stmin);
190   AliMedium(50, "FE_C2           ", 50, 0, fieldType, maxField, tmaxfd, stemax, deemax, epsil, stmin);
191
192   //    Copper 
193   AliMedium(11, "Cu_C0           ", 11, 0, fieldType, maxField, tmaxfd, stemax, deemax, epsil, stmin);
194   AliMedium(31, "Cu_C1           ", 31, 0, fieldType, maxField, tmaxfd, stemax, deemax, epsil, stmin);
195   AliMedium(51, "Cu_C2           ", 51, 0, fieldType, maxField, tmaxfd, stemax, deemax, epsil, stmin);
196   
197   //    Vacuum 
198   AliMedium(16, "VA_C0           ", 16, 0, fieldType, maxField, tmaxfd, stemax, deemax, epsil, stmin);
199   AliMedium(36, "VA_C1           ", 36, 0, fieldType, maxField, tmaxfd, stemax, deemax, epsil, stmin);
200   AliMedium(56, "VA_C2           ", 56, 0, fieldType, maxField, tmaxfd, stemax, deemax, epsil, stmin);
201   
202   //    Steel 
203   AliMedium(19, "ST_C0           ", 19, 0, fieldType, maxField, tmaxfd, stemax, deemax, epsil, stmin);
204   AliMedium(39, "ST_C1           ", 39, 0, fieldType, maxField, tmaxfd, stemax, deemax, epsil, stmin);
205   AliMedium(59, "ST_C3           ", 59, 0, fieldType, maxField, tmaxfd, stemax, deemax, epsil, stmin);
206
207    //
208    // Parameters  for AD scintillator: NE-102 (BC400)
209    //
210    // NE-102, has the following properties : (from internet, need better reference)
211    //    Density : ca. 1.03 g/cm3
212    //    Electrons/cm3: 3.39 x 10^23
213    //    H atoms/cm3: 5.28 x 10^22
214    //    C atoms/cm3: 4.78 x 10^22
215    //    Ratio of H to C : 1.104 .
216    //    wavelength of emission : ~4.23 nm.
217    //    Decay time : ~2.4 ns.
218    //    Luminescent efficiency : typically 18% of NaI(Tl)
219    //    Photons/MeV: 2.5 x 10^4 
220    //
221    // H                // C 
222    // as[0] = 1.00794;    as[1] = 12.011;
223    // zs[0] = 1.;         zs[1] = 6.;
224    // ws[0] = 5.23;       ws[1] = 4.74;
225    // density = 1.032;
226    // id      = 1;
227    // AliMixture( id, "NE102", as, zs, density, -2, ws );
228    // AliMedium( id, "NE102", id, 1, fieldType, maxField, maxBending, maxStepSize,
229    //            maxEnergyLoss, precision, minStepSize );
230
231    // ecalvovi@cern.ch
232    // Parameters  for AD scintillator: BC404
233    //
234    // NE-102, has the following properties : (from internet, need better reference)
235    //    Density : ca. 1.032 g/cm3
236    //    Electrons/cm3: 3.37 x 10^23
237    //    H atoms/cm3: 5.21 x 10^22
238    //    C atoms/cm3: 4.74 x 10^22
239    //    Ratio of H to C : 1.100 
240    //    wavelength of emission : 408 nm.
241    //    Decay time : 1.8 ns.
242    //    Luminescent efficiency : typically 18% of NaI(Tl)
243    //    Photons/MeV: ??
244    //
245    // H                // C 
246    as[0] = 1.00794;    as[1] = 12.011;
247    zs[0] = 1.;         zs[1] = 6.;
248    ws[0] = 5.21;       ws[1] = 4.74;
249    density = 1.032;
250    id      = 1;
251    AliMixture( id, "BC404", as, zs, density, -2, ws );
252    AliMedium ( id, "BC404", id, 1, fieldType, maxField, maxBending, maxStepSize,
253               maxEnergyLoss, precision, minStepSize );
254    // parameters AliMedium: numed  name   nmat   isvol  ifield fieldm tmaxfd stemax deemax epsil  stmin  
255    // ... 
256    // isvol       sensitive volume if isvol!=0
257    // ifield      magnetic field flag (see below)
258    // fieldm      maximum magnetic field
259    // ...
260    // ifield =  0       no magnetic field
261    //        = -1       user decision in guswim
262    //        =  1       tracking performed with Runge Kutta
263    //        =  2       tracking performed with helix
264    //        =  3       constant magnetic field along z
265
266
267    //
268    // Parameters for lightGuide:  
269    //     TODO check material 
270    // Should be Poly(methyl methacrylate) (PMMA) acrylic 
271    // (C5O2H8)n 
272    // Density  1.18 g/cm3
273    // Mixture PMMA    Aeff=12.3994 Zeff=6.23653 rho=1.18 radlen=34.0677 intlen=63.3073
274    // Element #0 : C  Z=  6.00 A= 12.01 w= 0.600 natoms=5
275    // Element #1 : H  Z=  1.00 A=  1.01 w= 0.081 natoms=8
276    // Element #2 : O  Z=  8.00 A= 16.00 w= 0.320 natoms=2
277
278    // Carbon          Hydrogen          Oxygen
279    as[0] = 12.0107;   as[1] = 1.00794;  as[2] = 15.9994;
280    zs[0] = 6.;        zs[1] = 1.;       zs[2] = 8.;
281    ws[0] = 0.60;      ws[1] = 0.081;    ws[2] = 0.32;
282    density = 1.18;
283    id      = 2;
284    AliMixture( id, "PMMA", as, zs, density, 3, ws );
285    AliMedium( id,"PMMA", id, 1, fieldType, maxField, maxBending, maxStepSize,
286               maxEnergyLoss, precision, minStepSize );
287
288    
289    // mu-metal
290    // Niquel          Iron              Molybdenum        Manganese
291    as[0] = 58.6934;   as[1] = 55.845;   as[2] = 95.94;    as[3] = 54.9380;  
292    zs[0] = 28.;       zs[1] = 26.;      zs[2] = 42.;      zs[3] = 25.;   
293    ws[0] = 0.802;     ws[1] = 0.14079;  ws[2] = 0.0485;   ws[3] = 0.005;
294    // Silicon         Chromium          Cobalt            Aluminium
295    as[4] = 28.0855;   as[5] = 51.9961;  as[6] = 58.9332;  as[7] = 26.981539;   
296    zs[4] = 14.;       zs[5] = 24.;      zs[6] = 27.;      zs[7] = 13.;   
297    ws[4] = 0.003;     ws[5] = 0.0002;   ws[6] = 0.0002;   ws[7] = 0.0001;
298    // Carbon          Phosphorus        Sulfur
299    as[8] = 12.0107;   as[9] = 30.97376; as[10] = 32.066;
300    zs[8] = 6.;        zs[9] = 15.;      zs[10] = 16.;
301    ws[8] = 0.00015;   ws[9] = 0.00005;  ws[10] = 0.00001;
302    density = 8.25;
303    id      = 3;
304    AliMixture( id, "MuMetal", as, zs, density, 11, ws );
305    AliMedium( id,"MuMetal", id, 1, fieldType, maxField, maxBending, maxStepSize,
306               maxEnergyLoss, precision, minStepSize );
307
308    // Parameters for ADCPMA: Aluminium
309    a_ad = 26.98; 
310    z_ad = 13.00;
311    density   = 2.7;
312    radLength = 8.9;
313    absLength = 37.2;
314    id = 4;
315    AliMaterial (id, "Alum",  a_ad, z_ad, density, radLength, absLength, 0, 0 );
316    AliMedium( id, "Alum", id, 1, fieldType, maxField, maxBending, maxStepSize,
317               maxEnergyLoss, precision, minStepSize );
318
319    // Parameters for ADCPMG: Glass for the simulation Aluminium 
320    // TODO fix material
321    a_ad = 26.98; 
322    z_ad = 13.00;
323    density   = 2.7;
324    radLength = 8.9;
325    absLength = 37.2;
326    id = 5;
327    AliMaterial( id, "Glass",  a_ad, z_ad, density, radLength, absLength, 0, 0 );
328    AliMedium( id, "Glass", id, 1, fieldType, maxField, maxBending, maxStepSize,
329               maxEnergyLoss, precision, minStepSize );
330
331
332 }
333 //_____________________________________________________________________________
334 void AliAD::SetTreeAddress()
335 {
336    //
337    // Sets tree address for hits.
338    //
339
340         TBranch *branch;
341         char branchname[20];
342         snprintf(branchname,19,"%s",GetName());
343         // Branch address for hit tree
344         TTree *treeH = fLoader->TreeH();
345         if (treeH ) 
346         {
347                 branch = treeH->GetBranch(branchname);
348                 if (branch) branch->SetAddress(&fHits);
349         }
350 }
351
352
353 //_____________________________________________________________________________
354 void AliAD::MakeBranch(Option_t* opt)
355 {
356         const char* oH = strstr(opt,"H");
357         if (fLoader->TreeH() && oH && (fHits==0x0))
358         {
359                 fHits = new TClonesArray("AliADhit",1000);
360                 fNhits = 0;
361         }
362         AliDetector::MakeBranch(opt);
363 }
364 //_____________________________________________________________________________
365 AliLoader* AliAD::MakeLoader(const char* topfoldername)
366
367  
368    AliDebug(1,Form("Creating AliADLoader, Top folder is %s ",topfoldername));
369    fLoader = new AliADLoader(GetName(),topfoldername);
370    return fLoader;
371 }
372
373 //_____________________________________________________________________________
374 AliDigitizer* AliAD::CreateDigitizer(AliDigitizationInput* digInput) const
375 {
376    //
377    // Creates a digitizer for AD
378    //
379    return new AliADDigitizer(digInput);
380 }
381 //_____________________________________________________________________________
382 void AliAD::Hits2Digits(){
383   //
384   // Converts hits to digits
385   //
386   // Creates the AD digitizer 
387   AliADDigitizer* dig = new AliADDigitizer(this,AliADDigitizer::kHits2Digits);
388
389   // Creates the digits
390   dig->Digitize("");
391
392   // deletes the digitizer
393   delete dig;
394 }
395
396 //_____________________________________________________________________________
397 void AliAD::Hits2SDigits(){
398   //
399   // Converts hits to summable digits
400   //
401   // Creates the AD digitizer 
402   AliADDigitizer* dig = new AliADDigitizer(this,AliADDigitizer::kHits2SDigits);
403
404   // Creates the sdigits
405   dig->Digitize("");
406
407   // deletes the digitizer
408   delete dig;
409 }
410
411 //_____________________________________________________________________________
412
413 void AliAD::Digits2Raw()
414 {
415            //
416    //  Converts digits of the current event to raw data
417    //
418    AliAD *fAD = (AliAD*)gAlice->GetDetector("AD");
419    fLoader->LoadDigits();
420    TTree* digits = fLoader->TreeD();
421    if (!digits) {
422       Error("Digits2Raw", "no digits tree");
423       return;
424    }
425    TClonesArray * ADdigits = new TClonesArray("AliADdigit",1000);
426    fAD->SetTreeAddress();               
427    digits->GetBranch("ADDigit")->SetAddress(&ADdigits); 
428   
429    const char *fileName    = AliDAQ::DdlFileName("AD",0);
430    AliADBuffer* buffer  = new AliADBuffer(fileName);
431    
432    // Get Trigger information first
433    // Read trigger inputs from trigger-detector object
434    AliDataLoader * dataLoader = fLoader->GetDigitsDataLoader();
435    if( !dataLoader->IsFileOpen() ) 
436         dataLoader->OpenFile("READ");
437    AliTriggerDetector* trgdet = (AliTriggerDetector*)dataLoader->GetDirectory()->Get("Trigger");
438    UInt_t triggerInfo = 0;
439    if(trgdet) {
440       triggerInfo = trgdet->GetMask() & 0xffff;
441    }
442    else {
443       AliError(Form("There is no trigger object for %s",fLoader->GetName()));
444    }
445
446    buffer->WriteTriggerInfo((UInt_t)triggerInfo); 
447    buffer->WriteTriggerScalers(); 
448    buffer->WriteBunchNumbers(); 
449   
450    // Now retrieve the channel information: charge smaples + time and 
451    // dump it into ADC and Time arrays
452    Int_t nEntries = Int_t(digits->GetEntries());
453    Short_t aADC[16][kNClocks];
454    Float_t aTime[16];
455    Float_t aWidth[16];
456    Bool_t  aIntegrator[16];
457    Bool_t  aBBflag[16];
458    Bool_t  aBGflag[16];
459   
460    for (Int_t i = 0; i < nEntries; i++) {
461      fAD->ResetDigits();
462      digits->GetEvent(i);
463      Int_t ndig = ADdigits->GetEntriesFast(); 
464    
465      if(ndig == 0) continue;
466      for(Int_t k=0; k<ndig; k++){
467          AliADdigit* fADDigit = (AliADdigit*) ADdigits->At(k);
468
469          Int_t iChannel       = fADDigit->PMNumber();
470          
471          for(Int_t iClock = 0; iClock < kNClocks; ++iClock) aADC[iChannel][iClock] = fADDigit->ChargeADC(iClock);
472          aTime[iChannel]      = fADDigit->Time(); //divide by resolution
473          aWidth[iChannel]     = fADDigit->Width(); //divide by resolution
474          aIntegrator[iChannel]= fADDigit->Integrator();
475          aBBflag[iChannel]    = fADDigit->GetBBflag();
476          aBGflag[iChannel]    = fADDigit->GetBGflag();
477          
478          //AliDebug(1,Form("DDL: %s\tdigit number: %d\tPM number: %d\tADC: %d\tTime: %f",
479                          //fileName,k,fADDigit->PMNumber(),aADC[iChannel][AliADdigit::kNClocks/2],aTime[iChannel])); 
480      }        
481    }
482
483    // Now fill raw data 
484    Int_t iCIU=0;
485    for (Int_t  iV0CIU = 0; iV0CIU < 8; iV0CIU++) {
486    
487       if(iV0CIU != 2 && iV0CIU != 5) {
488         buffer->WriteEmptyCIU();
489         continue;
490        }
491       // decoding of one Channel Interface Unit numbered iCIU - there are 8 channels per CIU (and 2 CIUs) : 
492       for(Int_t iChannel_Offset = iCIU*8; iChannel_Offset < (iCIU*8)+8; iChannel_Offset=iChannel_Offset+4) { 
493          for(Int_t iChannel = iChannel_Offset; iChannel < iChannel_Offset+4; iChannel++) {
494              buffer->WriteChannel(iChannel, aADC[iChannel], aIntegrator[iChannel]);       
495          }
496          buffer->WriteBeamFlags(&aBBflag[iChannel_Offset],&aBGflag[iChannel_Offset]); 
497          buffer->WriteMBInfo(); 
498          buffer->WriteMBFlags();   
499          buffer->WriteBeamScalers(); 
500       } 
501
502       for(Int_t iChannel = iCIU*8 + 7; iChannel >= iCIU*8; iChannel--) {
503           buffer->WriteTiming(aTime[iChannel], aWidth[iChannel]); 
504       } // End of decoding of one CIU card
505       iCIU++;    
506   } // end of decoding the eight CIUs
507
508      
509   delete buffer;
510   fLoader->UnloadDigits();  
511 }
512
513 //_____________________________________________________________________________
514
515 Bool_t AliAD::Raw2SDigits(AliRawReader* rawReader)
516 {
517         // reads raw data to produce sdigits
518         // for AD not implemented yet 
519         return kTRUE;
520 }