Splitted libraries
[u/mrichter/AliRoot.git] / ZDC / AliZDCDigitizer.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$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 //                      ZDC digitizer class                                  //
21 //                                                                           //
22 ///////////////////////////////////////////////////////////////////////////////
23
24 #include <stdlib.h>
25
26 // --- ROOT system
27 #include <TTree.h>
28 #include <TFile.h>
29 #include <TNtuple.h>
30 #include <TRandom.h>
31
32 // --- AliRoot header files
33 #include "AliLog.h"
34 #include "AliRun.h"
35 #include "AliHeader.h"
36 #include "AliGenHijingEventHeader.h"
37 #include "AliRunDigitizer.h"
38 #include "AliRunLoader.h"
39 #include "AliZDCSDigit.h"
40 #include "AliZDCDigit.h"
41 #include "AliZDCFragment.h"
42 #include "AliZDCDigitizer.h"
43
44 ClassImp(AliZDCDigitizer)
45
46
47 //____________________________________________________________________________
48 AliZDCDigitizer::AliZDCDigitizer()
49 {
50 // Default constructor    
51
52 }
53
54 //____________________________________________________________________________
55 AliZDCDigitizer::AliZDCDigitizer(AliRunDigitizer* manager):
56   AliDigitizer(manager)
57 {
58 // Constructor    
59
60 }
61
62 //____________________________________________________________________________
63 AliZDCDigitizer::~AliZDCDigitizer()
64 {
65 // Destructor
66
67 }
68
69
70 //____________________________________________________________________________
71 Bool_t AliZDCDigitizer::Init()
72 {
73 // Initialize the digitizer
74
75   for (Int_t i = 0; i < 3; i++) {
76     for(Int_t j = 0; j < 5; j++) {
77       fPMGain[i][j]   = 100000.;
78     }
79   }
80
81 // ADC Caen V965
82   fADCRes[0]   = 0.0000008; // ADC Resolution high gain: 200 fC/adcCh
83   fADCRes[1]   = 0.0000064; // ADC Resolution low gain:  25  fC/adcCh
84
85   return kTRUE;
86 }
87
88 //____________________________________________________________________________
89 void AliZDCDigitizer::Exec(Option_t* /*option*/)
90 {
91 // Execute digitization
92
93   Float_t pm[3][5];
94   for (Int_t iSector1 = 0; iSector1 < 3; iSector1++)
95     for (Int_t iSector2 = 0; iSector2 < 5; iSector2++) {
96       pm[iSector1][iSector2] = 0;
97     }
98   // pm[0][...] = light in ZN [C, Q1, Q2, Q3, Q4]
99   // pm[1][...] = light in ZP [C, Q1, Q2, Q3, Q4]
100   // pm[2][...] = light in ZEM [x, 1, 2, x, x]
101
102   // impact parameter and number of spectators
103   Float_t impPar = -1;
104   Int_t specN = 0;
105   Int_t specP = 0;
106
107   // loop over input streams
108   for (Int_t iInput = 0; iInput < fManager->GetNinputs(); iInput++) {
109
110     // get run loader and ZDC loader
111     AliRunLoader* runLoader = 
112       AliRunLoader::GetRunLoader(fManager->GetInputFolderName(iInput));
113     AliLoader* loader = runLoader->GetLoader("ZDCLoader");
114     if (!loader) continue;
115
116     // load sdigits
117     loader->LoadSDigits();
118     TTree* treeS = loader->TreeS();
119     if (!treeS) continue;
120     AliZDCSDigit sdigit;
121     AliZDCSDigit* psdigit = &sdigit;
122     treeS->SetBranchAddress("ZDC", &psdigit);
123
124     // loop over sdigits
125     for (Int_t iSDigit = 0; iSDigit < treeS->GetEntries(); iSDigit++) {
126       treeS->GetEntry(iSDigit);
127       if (!psdigit) continue;
128
129       if ((sdigit.GetSector(1) < 0) || (sdigit.GetSector(1) > 4)) {
130         AliError(Form("\nsector[0] = %d, sector[1] = %d\n", 
131                       sdigit.GetSector(0), sdigit.GetSector(1)));
132         continue;
133       }
134       pm[sdigit.GetSector(0)-1][sdigit.GetSector(1)] += sdigit.GetLightPM();
135       //printf("\n\t sector[0] = %d, sector[1] = %d, pm[%d][%d] = %f \n",
136       //          sdigit.GetSector(0), sdigit.GetSector(1),sdigit.GetSector(0)-1,
137       //          sdigit.GetSector(1), pm[sdigit.GetSector(0)-1][sdigit.GetSector(1)]); // Chiara debugging!
138     }
139
140     // unload sdigits
141     loader->UnloadSDigits();
142
143     // get the impact parameter and the number of spectators in case of hijing
144     if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
145     AliHeader* header = runLoader->GetAliRun()->GetHeader();
146     if (!header) continue;
147     AliGenEventHeader* genHeader = header->GenEventHeader();
148     if (!genHeader) continue;
149     if (!genHeader->InheritsFrom(AliGenHijingEventHeader::Class())) continue;
150     impPar = ((AliGenHijingEventHeader*) genHeader)->ImpactParameter();
151     // Until there is only 1 ZDC set the # of spectators must be divided by 2!!!
152     specN = ((AliGenHijingEventHeader*) genHeader)->Spectatorsn() / 2;
153     specP = ((AliGenHijingEventHeader*) genHeader)->Spectatorsp() / 2;
154     AliDebug(2, Form("\n b = %f fm, Nspecn = %d, Nspecp = %d\n",
155                      impPar, specN, specP));
156   }
157
158   // add spectators
159   if (impPar >= 0) {
160     Int_t freeSpecN, freeSpecP;
161     Fragmentation(impPar, specN, specP, freeSpecN, freeSpecP);
162     SpectatorSignal(1, freeSpecN, pm);
163     SpectatorSignal(2, freeSpecP, pm);
164   }
165
166   // get the output run loader and loader
167   AliRunLoader* runLoader = 
168     AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
169   AliLoader* loader = runLoader->GetLoader("ZDCLoader");
170   if (!loader) {
171     AliError("no ZDC loader found");
172     return;
173   }
174
175   // create the output tree
176   const char* mode = "update";
177   if (runLoader->GetEventNumber() == 0) mode = "recreate";
178   loader->LoadDigits(mode);
179   loader->MakeTree("D");
180   TTree* treeD = loader->TreeD();
181   AliZDCDigit digit;
182   AliZDCDigit* pdigit = &digit;
183   const Int_t kBufferSize = 4000;
184   treeD->Branch("ZDC", "AliZDCDigit", &pdigit, kBufferSize);
185
186   // Create digits
187   Int_t sector[2];
188   Int_t digi[2];
189   for (sector[0] = 1; sector[0] <= 3; sector[0]++)
190     for (sector[1] = 0; sector[1] < 5; sector[1]++)  {
191         if ((sector[0] == 3) && ((sector[1] < 1) || (sector[1] > 2))) continue;
192         for (Int_t res = 0; res < 2; res++){
193           digi[res] = Phe2ADCch(1, 0, pm[sector[0]-1][sector[1]], res) 
194                     + Pedestal();
195           //printf("\n\t sector[0] = %d, sector[1] = %d, digi[%d] = %d\n",
196           //    sector[0], sector[1], res, digi[res]); // Chiara debugging!
197         }
198         new(pdigit) AliZDCDigit(sector, digi);
199         treeD->Fill();
200     }
201
202   // write the output tree
203   loader->WriteDigits("OVERWRITE");
204   loader->UnloadDigits();
205 }
206
207
208 //_____________________________________________________________________________
209 void AliZDCDigitizer::Fragmentation(Float_t impPar, Int_t specN, Int_t specP,
210                                     Int_t &freeSpecN, Int_t &freeSpecP) const
211 {
212 // simulate fragmentation of spectators
213
214   Int_t zz[100], nn[100];
215   AliZDCFragment frag(impPar);
216   for (Int_t j=0; j<=99; j++){
217      zz[j] =0;
218      nn[j] =0;
219   }
220
221   // Fragments generation
222   Int_t nAlpha;
223   frag.GenerateIMF(zz, nAlpha);
224
225   // Attach neutrons
226   Int_t ztot=0;
227   Int_t ntot=0;
228   frag.AttachNeutrons(zz, nn, ztot, ntot);
229   freeSpecN = specN-ntot-2*nAlpha;
230   freeSpecP = specP-ztot-2*nAlpha;
231   if(freeSpecN<0) freeSpecN=0;
232   if(freeSpecP<0) freeSpecP=0;
233   AliDebug(2, Form("FreeSpn = %d, FreeSpp = %d", freeSpecN, freeSpecP));
234 }
235
236 //_____________________________________________________________________________
237 void AliZDCDigitizer::SpectatorSignal(Int_t SpecType, Int_t numEvents, 
238                                       Float_t pm[3][5]) const
239 {
240 // add signal of the spectators
241  
242   TFile* file = NULL;
243   if (SpecType == 1) {          // --- Signal for spectator neutrons
244     file = TFile::Open("$ALICE/$ALICE_LEVEL/ZDC/ZNsignalntu.root");
245   } else if (SpecType == 2) {   // --- Signal for spectator protons
246     file = TFile::Open("$ALICE/$ALICE_LEVEL/ZDC/ZPsignalntu.root");
247   }
248   if (!file || !file->IsOpen()) {
249     AliError("Opening of file failed");
250     return;
251   }
252
253   TNtuple* zdcSignal = (TNtuple*) file->Get("ZDCSignal");
254   Int_t nentries = (Int_t) zdcSignal->GetEntries();
255   
256   Float_t *entry, hitsSpec[7];
257   Int_t pl, i, j, k, iev=0, rnd[125], volume[2];
258   for(pl=0;pl<125;pl++){
259      rnd[pl] = 0;
260   }
261   if (numEvents > 125) {
262     AliWarning(Form("numEvents (%d) is larger than 125", numEvents));
263     numEvents = 125;
264   }
265   for(pl=0;pl<numEvents;pl++){
266      rnd[pl] = (Int_t) (9999*gRandom->Rndm());
267      if(rnd[pl] >= 9998) rnd[pl] = 9997;
268      //printf(" rnd[%d] = %d\n",pl,rnd[pl]);     
269   }
270   // Sorting vector in ascending order with C function QSORT 
271   qsort((void*)rnd,numEvents,sizeof(Int_t),comp);
272   do{
273      for(i=0; i<nentries; i++){  
274         zdcSignal->GetEvent(i);
275         entry = zdcSignal->GetArgs();
276         if(entry[0] == rnd[iev]){
277           for(k=0; k<2; k++) volume[k] = (Int_t) entry[k+1];
278           for(j=0; j<7; j++) hitsSpec[j] = entry[j+3];
279
280           Float_t lightQ = hitsSpec[4];
281           Float_t lightC = hitsSpec[5];
282           AliDebug(3, Form("volume = (%d, %d), lightQ = %.3f, lightC = %.3f",
283                            volume[0], volume[1], lightQ, lightC));
284           if (volume[0] < 3) {  // ZN or ZP
285             pm[volume[0]-1][0] += lightC;
286             pm[volume[0]-1][volume[1]] += lightQ;
287           } else {
288             if (volume[1] == 1) pm[2][1] += lightC;
289             else                pm[2][2] += lightQ;
290         }
291
292
293         }
294         else if(entry[0] > rnd[iev]){
295           iev++;
296           continue;
297         }
298      }
299   }while(iev<numEvents);
300   
301   file->Close();
302   delete file;
303 }
304
305
306 //_____________________________________________________________________________
307 Int_t AliZDCDigitizer::Phe2ADCch(Int_t Det, Int_t Quad, Float_t Light, 
308                                  Int_t Res) const
309 {
310 // Evaluation of the ADC channel corresponding to the light yield Light
311
312   return (Int_t) (Light * fPMGain[Det-1][Quad] * fADCRes[Res]);
313 }
314
315 //_____________________________________________________________________________
316 Int_t AliZDCDigitizer::Pedestal() const
317 {
318   // --- Pedestal value -> extracted from a gaussian distribution
319   // obtained from the beam test on the ZEM prototype (Aug. 2000)
320
321   Float_t pedMean  = 50.;
322   Float_t pedWidth = 5.;
323
324   return (Int_t) gRandom->Gaus(pedMean, pedWidth);
325 }