Fixing a bug
[u/mrichter/AliRoot.git] / ZDC / ZDCsim / AliZDC.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 //                      Zero Degree Calorimeter                              //
21 //           This class contains the basic functions for the ZDCs;           //
22 //            functions specific to one particular geometry are              //
23 //                      contained in the derived classes                     //
24 //                                                                           //
25 ///////////////////////////////////////////////////////////////////////////////
26
27 // --- ROOT system
28 #include <TClonesArray.h>
29 #include <TTree.h>
30 #include <TFile.h>
31 #include <TSystem.h>
32 #include <TRandom.h>
33 #include <TParticle.h>
34
35 // --- AliRoot header files
36 #include "AliDetector.h"
37 #include "AliRawDataHeaderSim.h"
38 #include "AliRawReader.h"
39 #include "AliLoader.h"
40 #include "AliRun.h"
41 #include "AliMC.h"
42 #include "AliLog.h"
43 #include "AliDAQ.h"
44 #include "AliZDC.h"
45 #include "AliZDCHit.h"
46 #include "AliZDCSDigit.h"
47 #include "AliZDCDigit.h"
48 #include "AliZDCDigitizer.h"
49 #include "AliZDCRawStream.h"
50 #include "AliZDCPedestals.h"
51 #include "AliZDCEnCalib.h"
52 #include "AliZDCTowerCalib.h"
53 #include "AliFstream.h"
54
55  
56 ClassImp(AliZDC)
57
58 //_____________________________________________________________________________
59 AliZDC::AliZDC() :
60   AliDetector(),
61   fNoShower(0),
62   fPedCalib(0),
63   fEnCalibData(0),
64   fTowCalibData(0),
65   fZDCCalibFName(""),
66   fSpectatorTracked(1),
67   fBeamEnergy(0.),
68   fIspASystem(kFALSE),
69   fIsRELDISgen(kFALSE),
70   fOnlyZEM(kFALSE),
71   fFindMother(kFALSE)
72 {
73   //
74   // Default constructor for the Zero Degree Calorimeter base class
75   //
76   
77   fIshunt = 1;
78   fNhits  = 0;
79   fHits = 0;
80   fDigits = 0;
81   fNdigits = 0;
82   
83 }
84  
85 //_____________________________________________________________________________
86 AliZDC::AliZDC(const char *name, const char *title) : 
87   AliDetector(name,title),
88   fNoShower(0),
89   fPedCalib(0),
90   fEnCalibData(0),
91   fTowCalibData(0),
92   fZDCCalibFName(""),
93   fSpectatorTracked(1),
94   fBeamEnergy(0.),
95   fIspASystem(kFALSE),
96   fIsRELDISgen(kFALSE),
97   fOnlyZEM(kFALSE),
98   fFindMother(kFALSE)
99 {
100   //
101   // Standard constructor for the Zero Degree Calorimeter base class
102   //
103     
104   fIshunt = 1;
105   fNhits  = 0;
106   fDigits = 0;
107   fNdigits = 0;
108  
109   fHits = new TClonesArray("AliZDCHit",1000);
110   gAlice->GetMCApp()->AddHitList(fHits);
111   
112   SetName("ZDC"); SetTitle("ZDC");
113
114 }
115
116 //____________________________________________________________________________ 
117 AliZDC::~AliZDC()
118 {
119   //
120   // ZDC destructor
121   //
122
123   fIshunt = 0;
124   if(fPedCalib) delete fPedCalib;
125   if(fEnCalibData) delete fEnCalibData;
126   if(fEnCalibData) delete fEnCalibData;
127
128 }
129
130 //_____________________________________________________________________________
131 AliZDC::AliZDC(const AliZDC& ZDC) :
132 AliDetector("ZDC","ZDC"),
133 fNoShower(ZDC.fNoShower),
134 fPedCalib(ZDC.fPedCalib),
135 fEnCalibData(ZDC.fEnCalibData),
136 fTowCalibData(ZDC.fTowCalibData),
137 fZDCCalibFName(ZDC.fZDCCalibFName),
138 fSpectatorTracked(ZDC.fSpectatorTracked),
139 fBeamEnergy(ZDC.fBeamEnergy),
140 fIspASystem(ZDC.fIspASystem),
141 fIsRELDISgen(ZDC.fIsRELDISgen),
142 fOnlyZEM(ZDC.fOnlyZEM),
143 fFindMother(ZDC.fFindMother)
144 {
145   // copy constructor
146 }
147
148 //_____________________________________________________________________________
149 AliZDC& AliZDC::operator=(const AliZDC& ZDC)
150 {
151   // assignement operator
152   if(this!=&ZDC){
153     fNoShower = ZDC.fNoShower;
154     fPedCalib = ZDC.fPedCalib;
155     fEnCalibData = ZDC.fEnCalibData;
156     fTowCalibData = ZDC.fTowCalibData;
157     fZDCCalibFName = ZDC.fZDCCalibFName;
158     fBeamEnergy = ZDC.fBeamEnergy;
159     fIspASystem = ZDC.fIspASystem;
160     fIsRELDISgen = ZDC.fIsRELDISgen;
161     fOnlyZEM = ZDC.fOnlyZEM;
162     fFindMother = ZDC.fFindMother;
163   } return *this;
164 }
165
166 //_____________________________________________________________________________
167 void AliZDC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
168 {
169   //
170   //            Add a ZDC hit to the hit list.
171   
172   static Float_t trackTime=0., trackEta=0., primKinEn=0., xImpact=0., yImpact=0., sFlag=0.;
173   static Int_t   pcPDGcode=0, motPDGcode=0;
174
175   AliZDCHit *newquad, *curprimquad;
176   newquad = new AliZDCHit(fIshunt, track, vol, hits);
177   TClonesArray &lhits = *fHits;
178   
179   if(fNhits==0){
180       // First hit -> setting flag for primary or secondary particle
181       TParticle * p = gAlice->GetMCApp()->Particle(track);
182       Int_t imo = p->GetFirstMother();
183       //
184       if(track != imo){
185         newquad->SetSFlag(1);  // SECONDARY particle entering the ZDC
186       }
187       else if(track == imo){
188         newquad->SetSFlag(0);  // PRIMARY particle entering the ZDC
189       }
190       //  
191       sFlag      = newquad->GetSFlag();
192       primKinEn  = newquad->GetPrimKinEn();
193       xImpact    = newquad->GetXImpact();
194       yImpact    = newquad->GetYImpact();
195       pcPDGcode  = newquad->GetPDGCode();
196       motPDGcode = newquad->GetMotherPDGCode();
197       trackTime  = newquad->GetTrackTOF();
198       trackEta   = newquad->GetTrackEta();
199
200    }
201    else{       
202       newquad->SetPrimKinEn(primKinEn);
203       newquad->SetXImpact(xImpact);
204       newquad->SetYImpact(yImpact);
205       newquad->SetSFlag(sFlag);
206       newquad->SetPDGCode(pcPDGcode);
207       newquad->SetMotherPDGCode(motPDGcode);
208       newquad->SetTrackTOF(trackTime);
209       newquad->SetTrackEta(trackEta);
210    }
211   
212   Int_t j;
213   for(j=0; j<fNhits; j++){
214     // If hits are equal (same track, same volume), sum them.
215      curprimquad = (AliZDCHit*) lhits[j];
216      if(*curprimquad == *newquad){
217         *curprimquad = *curprimquad+*newquad;
218         // Ch. debug
219         //printf("\n\t Summing hits **************** \n", fNhits);
220         //curprimquad->Print("");
221         //
222         delete newquad;
223         return;
224      } 
225   }
226   //printf( "PDG from hits[10] = %f\n", hits[10]);
227
228     //Otherwise create a new hit
229     new(lhits[fNhits]) AliZDCHit(*newquad);
230     fNhits++;
231     // Ch. debug
232     //printf("\n\t New ZDC hit added! fNhits = %d\n", fNhits);
233     //newquad->Print("");
234     
235     delete newquad;
236 }
237
238 //____________________________________________________________________________
239 Float_t AliZDC::ZMin(void) const
240 {
241   // Minimum dimension of the ZDC module in z
242   return -11600.;
243 }
244
245 //____________________________________________________________________________
246 Float_t AliZDC::ZMax(void) const
247 {
248   // Maximum dimension of the ZDC module in z
249   return 11750.;
250 }
251   
252
253 //_____________________________________________________________________________
254 void AliZDC::MakeBranch(Option_t *opt)
255 {
256   //
257   // Create Tree branches for the ZDC
258   //
259
260   char branchname[10];
261   snprintf(branchname, 10, "%s", GetName());
262
263   const char *cH = strstr(opt,"H");
264   
265   if(cH && fLoader->TreeH()) {
266     if (fHits) {
267       fHits->Clear();
268       fNhits = 0;
269     }
270     else {
271       fHits   = new TClonesArray("AliZDCHit",1000); 
272       if (gAlice && gAlice->GetMCApp())
273         gAlice->GetMCApp()->AddHitList(fHits);
274     }
275   }
276   
277   AliDetector::MakeBranch(opt);
278 }
279
280 //_____________________________________________________________________________
281 void AliZDC::Hits2SDigits()
282 {
283   // Create summable digits from hits
284   
285   AliDebug(1,"\n        AliZDC::Hits2SDigits() ");
286   
287   fLoader->LoadHits("read");
288   fLoader->LoadSDigits("recreate");
289   AliRunLoader* runLoader = fLoader->GetRunLoader();
290   AliZDCSDigit sdigit;
291   AliZDCSDigit* psdigit = &sdigit;
292
293   // Event loop
294   for(Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
295     Float_t pmZNC[5], pmZPC[5], pmZNA[5], pmZPA[5], pmZEM1=0., pmZEM2=0.;
296     for(Int_t i=0; i<5; i++) pmZNC[i] = pmZPC[i] =  pmZNA[i] = pmZPA[i] = 0;
297
298     runLoader->GetEvent(iEvent);
299     TTree* treeH = fLoader->TreeH();
300     Int_t ntracks = (Int_t) treeH->GetEntries();
301     ResetHits();
302
303     // Tracks loop
304     Int_t sector[2]; Float_t trackTime = 0.;
305     for(Int_t itrack = 0; itrack < ntracks; itrack++) {
306       treeH->GetEntry(itrack);
307       for(AliZDCHit* zdcHit = (AliZDCHit*)FirstHit(-1); zdcHit;
308           zdcHit = (AliZDCHit*)NextHit()) { 
309                       
310         sector[0] = zdcHit->GetVolume(0);
311         sector[1] = zdcHit->GetVolume(1);
312         if((sector[1] < 1) || (sector[1]>5)) {
313           Error("Hits2SDigits", "sector[0] = %d, sector[1] = %d", sector[0], sector[1]);
314           continue;
315         }
316         Float_t lightQ = zdcHit->GetLightPMQ();
317         Float_t lightC = zdcHit->GetLightPMC();
318         trackTime = zdcHit->GetTrackTOF();
319         // Signals from ZEM are delayed to arrive in time with ZDC signals
320         if(sector[0] == 3)  trackTime += 320;
321         // Ch. debug
322         //printf("\t det %d vol %d trackTOF %f lightQ %1.0f lightC %1.0f\n",
323         //      sector[0], sector[1], trackTime, lightQ, lightC);
324      
325         if(sector[0] == 1) { //ZNC 
326           pmZNC[0] += lightC;
327           pmZNC[sector[1]] += lightQ;
328         } 
329         else if(sector[0] == 2) { //ZPC 
330           pmZPC[0] += lightC;
331           pmZPC[sector[1]] += lightQ;
332         } 
333         else if(sector[0] == 3) { //ZEM 
334           if(sector[1] == 1) pmZEM1 += lightC;
335           else pmZEM2 += lightQ;
336         }
337         if(sector[0] == 4) { //ZNA 
338           pmZNA[0] += lightC;
339           pmZNA[sector[1]] += lightQ;
340         } 
341         else if(sector[0] == 5) { //ZPA 
342           pmZPA[0] += lightC;
343           pmZPA[sector[1]] += lightQ;
344         } 
345       }//Hits loop
346     }//Tracks loop
347
348     // create the output tree
349     fLoader->MakeTree("S");
350     TTree* treeS = fLoader->TreeS();
351     const Int_t kBufferSize = 4000;
352     treeS->Branch(GetName(), "AliZDCSDigit", &psdigit, kBufferSize);
353
354     // Create sdigits for ZNC
355     sector[0] = 1; // Detector = ZNC
356     for(Int_t j = 0; j < 5; j++) {
357       sector[1] = j; 
358       if(pmZNC[j]>0){
359         new(psdigit) AliZDCSDigit(sector, pmZNC[j], trackTime);
360         treeS->Fill();
361         // Ch. debug
362         //printf("\t SDigit created: det %d quad %d pmZNC[%d] %1.0f trackTOF %f\n",
363         //      sector[0], sector[1], j, pmZNC[j], trackTime);
364       }
365     }
366   
367     // Create sdigits for ZPC
368     sector[0] = 2; // Detector = ZPC
369     for(Int_t j = 0; j < 5; j++) {
370       sector[1] = j; // Towers PM ADCs
371       if(pmZPC[j]>0){
372         new(psdigit) AliZDCSDigit(sector, pmZPC[j], trackTime);
373         treeS->Fill();
374         // Ch. debug
375         //printf("\t SDigit created: det %d quad %d pmZPC[%d] %1.0f trackTOF %f\n",
376         //      sector[0], sector[1], j, pmZPC[j], trackTime);
377       }
378     }
379
380     // Create sdigits for ZEM
381     sector[0] = 3; 
382     sector[1] = 1; // Detector = ZEM1
383     if(pmZEM1>0){ 
384       new(psdigit) AliZDCSDigit(sector, pmZEM1, trackTime);
385       treeS->Fill();
386       // Ch. debug
387       //printf("\t SDigit created: det %d quad %d pmZEM1 %1.0f trackTOF %f\n",
388       //        sector[0], sector[1], pmZEM1, trackTime);
389     }
390     sector[1] = 2; // Detector = ZEM2
391     if(pmZEM2>0){
392       new(psdigit) AliZDCSDigit(sector, pmZEM2, trackTime);
393       treeS->Fill();
394       // Ch. debug
395       //printf("\t SDigit created: det %d quad %d pmZEM2 %1.0f trackTOF %f\n",
396       //        sector[0], sector[1], pmZEM2, trackTime);
397     }
398
399     // Create sdigits for ZNA
400     sector[0] = 4; // Detector = ZNA
401     for(Int_t j = 0; j < 5; j++) {
402       sector[1] = j; // Towers PM ADCs
403       if(pmZNA[j]>0){
404         new(psdigit) AliZDCSDigit(sector, pmZNA[j], trackTime);
405         treeS->Fill();
406         // Ch. debug
407         //printf("\t SDigit created: det %d quad %d pmZNA[%d] %1.0f trackTOF %f\n",
408         //      sector[0], sector[1], j, pmZNA[j], trackTime);
409       }
410     }
411   
412     // Create sdigits for ZPA
413     sector[0] = 5; // Detector = ZPA
414     sector[1] = 0; // Common PM ADC
415     for(Int_t j = 0; j < 5; j++) {
416       sector[1] = j; // Towers PM ADCs
417       if(pmZPA[j]>0){
418         new(psdigit) AliZDCSDigit(sector, pmZPA[j], trackTime);
419         treeS->Fill();
420         // Ch. debug
421         //printf("\t SDigit created: det %d quad %d pmZPA[%d] %1.0f trackTOF %f\n",
422         //      sector[0], sector[1], j, pmZPA[j], trackTime);
423       }
424     }
425
426     // write the output tree
427     fLoader->WriteSDigits("OVERWRITE");
428   }
429
430   fLoader->UnloadHits();
431   fLoader->UnloadSDigits();
432 }
433
434 //_____________________________________________________________________________
435 AliDigitizer* AliZDC::CreateDigitizer(AliDigitizationInput* digInput) const{
436   // Create the digitizer for ZDC
437   AliZDCDigitizer *zdcDigitizer = new AliZDCDigitizer(digInput);
438   if(fSpectatorTracked==0) zdcDigitizer->SetSpectators2Track();
439   if(fBeamEnergy>0.01) zdcDigitizer->SetBeamEnergy(fBeamEnergy);
440   if(fIspASystem==kTRUE) zdcDigitizer->SetpAsystem();
441   if(fIsRELDISgen==kTRUE) zdcDigitizer->SetRELDISGenerator();
442   //if(fIspASystem==kTRUE) printf("\n **** ZDC digitizer initialized for p-A collisions\n\n");
443   return zdcDigitizer;
444 }
445
446 //_____________________________________________________________________________
447 void AliZDC::Digits2Raw()
448 {
449   // Convert ZDC digits to raw data
450
451   // Format: 24 int values -> ZN1(C+Q1-4), ZP1(C+Q1-4), ZEM1, ZEM2, ZN(C+Q1-4), ZP2(C+Q1-4), 2 Ref PMs
452   //         + 24 int values for the corresponding out of time channels
453   // For the CAEN module V965 we have an Header, the Data Words and an End Of Block
454   //    12 channels x 2 gain chains read from 1st ADC module
455   //    12 channels x 2 gain chains read from 2nd ADC module
456   //    12 channels x 2 gain chains read from 3rd ADC module (o.o.t.)
457   //    12 channels x 2 gain chains read from 4rth ADC module (o.o.t.)
458   //
459   const int knADCData1=12, knADCData2=12; 
460   const int knADCData3=12, knADCData4=12; 
461   //
462   UInt_t lADCHeader1; 
463   UInt_t lADCHeader2; 
464   UInt_t lADCHeader3; 
465   UInt_t lADCHeader4; 
466   //
467   UInt_t lADCData1[2*knADCData1];
468   UInt_t lADCData2[2*knADCData2];
469   UInt_t lADCData3[2*knADCData3];
470   UInt_t lADCData4[2*knADCData4];
471   //
472   UInt_t lADCEndBlock;
473
474   // load the digits
475   fLoader->LoadDigits("read");
476   AliZDCDigit digit;
477   AliZDCDigit* pdigit = &digit;
478   TTree* treeD = fLoader->TreeD();
479   if(!treeD) return;
480   treeD->SetBranchAddress("ZDC", &pdigit);
481   //printf("\t AliZDC::Digits2Raw -> TreeD has %d entries\n",(Int_t) treeD->GetEntries());
482
483   // Reading channel map
484   //printf("\n\t Reading ADC mapping from OCDB\n");
485   AliZDCChMap * chMap = GetChMap();
486   const int nCh = knADCData1+knADCData2+knADCData3+knADCData4;
487   Int_t  mapADC[nCh][4]; 
488   for(Int_t i=0; i<nCh; i++){
489     mapADC[i][0] = chMap->GetADCModule(i);
490     mapADC[i][1] = chMap->GetADCChannel(i);
491     mapADC[i][2] = chMap->GetDetector(i);
492     mapADC[i][3] = chMap->GetSector(i);
493     // Ch. debug
494     //printf("  mapADC[%d] = (%d %d %d %d)\n", i,
495     //  mapADC[i][0],mapADC[i][1],mapADC[i][2],mapADC[i][3]);
496   }
497
498   // *** Fill data array
499   // ** ADC header
500   UInt_t lADCHeaderGEO1 = 0;
501   UInt_t lADCHeaderGEO2 = 1;
502   UInt_t lADCHeaderGEO3 = 2;
503   UInt_t lADCHeaderGEO4 = 3;
504   UInt_t lADCHeaderCRATE = 0;
505   UInt_t lADCHeaderCNT1 = knADCData1;
506   UInt_t lADCHeaderCNT2 = knADCData2;
507   UInt_t lADCHeaderCNT3 = knADCData3;
508   UInt_t lADCHeaderCNT4 = knADCData4;
509     
510   lADCHeader1 = lADCHeaderGEO1 << 27 | 0x1 << 25 | lADCHeaderCRATE << 16 |
511                lADCHeaderCNT1 << 8 ;
512   lADCHeader2 = lADCHeaderGEO2 << 27 | 0x1 << 25 | lADCHeaderCRATE << 16 |
513                lADCHeaderCNT2 << 8 ;
514   lADCHeader3 = lADCHeaderGEO3 << 27 | 0x1 << 25 | lADCHeaderCRATE << 16 |
515                lADCHeaderCNT3 << 8 ;
516   lADCHeader4 = lADCHeaderGEO4 << 27 | 0x1 << 25 | lADCHeaderCRATE << 16 |
517                lADCHeaderCNT4 << 8 ;
518       
519   // ** ADC data word
520   UInt_t lADCDataGEO = 0;
521   //
522   UInt_t lADCDataValue1[2*knADCData1];
523   UInt_t lADCDataValue2[2*knADCData2];
524   UInt_t lADCDataValue3[2*knADCData3];
525   UInt_t lADCDataValue4[2*knADCData4];
526   //
527   UInt_t lADCDataOvFlwHG = 0;
528   UInt_t lADCDataOvFlwLG = 0;
529   //
530   for(Int_t i=0; i<2*knADCData1 ; i++) lADCDataValue1[i] = 0;
531   for(Int_t i=0; i<2*knADCData2 ; i++) lADCDataValue2[i] = 0;
532   for(Int_t i=0; i<2*knADCData3 ; i++) lADCDataValue3[i] = 0;
533   for(Int_t i=0; i<2*knADCData4 ; i++) lADCDataValue4[i] = 0;
534   //
535   UInt_t lADCDataChannel = 0;
536   
537   Int_t indADC0=0, indADC1=0, indADC2=0, indADC3=0;
538   
539   // loop over digits
540   for(Int_t iDigit=0; iDigit<(Int_t) (treeD->GetEntries()); iDigit++){
541     treeD->GetEntry(iDigit);
542     if(!pdigit) continue;
543     //digit.Print("");
544    
545     // *** ADC data
546     // Scan of the map to assign the correct ADC module-channel
547     for(Int_t k=0; k<nCh; k++){
548       if(iDigit<knADCData1+knADCData2){ 
549        if(digit.GetSector(0)==mapADC[k][2] && digit.GetSector(1)==mapADC[k][3]){
550          lADCDataGEO = (UInt_t) mapADC[k][0];
551          lADCDataChannel = (UInt_t) mapADC[k][1];
552          break;
553        } 
554       }
555       else{
556        if(digit.GetSector(0)==mapADC[k][2] && digit.GetSector(1)==mapADC[k][3]){
557          lADCDataGEO = (UInt_t) mapADC[k][0];
558          lADCDataChannel = (UInt_t) mapADC[k][1];
559          if(k>knADCData1+knADCData2) break;
560        } 
561       }
562     }
563     // Ch. debug
564     //printf("iDigit %d det %d sec %d -> lADCDataGEO %d  lADCDataChannel %d\n",
565     //  iDigit,digit.GetSector(0),digit.GetSector(1),lADCDataGEO,lADCDataChannel);
566      
567     if(lADCDataGEO==0){ 
568       if(indADC0>=knADCData1){
569         AliWarning(Form(" Problem with digit index %d for ADC0\n", indADC0));
570         return;
571       }
572       Int_t indLG = indADC0+knADCData1;
573       // High gain ADC ch.       
574       if(digit.GetADCValue(0) > 2047) lADCDataOvFlwHG = 1; 
575       lADCDataValue1[indADC0] = digit.GetADCValue(0);    
576       lADCData1[indADC0] = lADCDataGEO << 27 |  lADCDataChannel << 17 | 
577                     lADCDataOvFlwHG << 12 | (lADCDataValue1[indADC0] & 0xfff); 
578       // Low gain ADC ch.
579       if(digit.GetADCValue(1) > 2047) lADCDataOvFlwLG = 1; 
580       lADCDataValue1[indLG] = digit.GetADCValue(1); 
581       lADCData1[indLG] = lADCDataGEO << 27 |  lADCDataChannel << 17 | 0x1 << 16 |
582                     lADCDataOvFlwLG << 12 | (lADCDataValue1[indLG] & 0xfff);  
583       // Ch. debug
584       //printf(" lADCDataGEO %d  ADCdataHG[%d] %d  ADCdataLG[%d] %d\n", 
585       //  lADCDataGEO,indADC0,lADCDataValue1[indADC0],indLG,lADCDataValue1[indLG]);
586                     
587       indADC0++;
588     }
589     else if(lADCDataGEO==1){ 
590       if(indADC1>=knADCData2){
591          AliWarning(Form(" Problem with digit index %d for ADC1\n", indADC1));
592          return;
593       }
594       Int_t indLG = indADC1+knADCData2;
595       // High gain ADC ch.       
596       if(digit.GetADCValue(0) > 2047) lADCDataOvFlwHG = 1; 
597       lADCDataValue2[indADC1] = digit.GetADCValue(0);    
598       lADCData2[indADC1] = lADCDataGEO << 27 | lADCDataChannel << 17 | 
599                     lADCDataOvFlwHG << 12 | (lADCDataValue2[indADC1] & 0xfff); 
600       // Low gain ADC ch.
601       if(digit.GetADCValue(1) > 2047) lADCDataOvFlwLG = 1; 
602       lADCDataValue2[indLG] = digit.GetADCValue(1); 
603       lADCData2[indLG] = lADCDataGEO << 27 |  lADCDataChannel << 17 | 0x1 << 16 |
604                     lADCDataOvFlwLG << 12 | (lADCDataValue2[indLG] & 0xfff);  
605       // Ch. debug
606       //printf(" lADCDataGEO %d  ADCdataHG[%d] %d  ADCdataLG[%d] %d\n", 
607       //  lADCDataGEO,indADC1,lADCDataValue2[indADC1],indLG,lADCDataValue2[indLG]);
608                   
609       indADC1++;
610     }
611     else if(lADCDataGEO==2){ 
612       if(indADC2>=knADCData3){
613         AliWarning(Form(" Problem with digit index %d for ADC2\n", indADC2));
614         return;
615       }
616       Int_t indLG = indADC2+knADCData3;
617       // High gain ADC ch.       
618       if(digit.GetADCValue(0) > 2047) lADCDataOvFlwHG = 1; 
619       lADCDataValue3[indADC1] = digit.GetADCValue(0);    
620       lADCData3[indADC1] = lADCDataGEO << 27 | lADCDataChannel << 17 | 
621                     lADCDataOvFlwHG << 12 | (lADCDataValue3[indADC2] & 0xfff); 
622       // Low gain ADC ch.
623       if(digit.GetADCValue(1) > 2047) lADCDataOvFlwLG = 1; 
624       lADCDataValue3[indLG] = digit.GetADCValue(1); 
625       lADCData3[indLG] = lADCDataGEO << 27 |  lADCDataChannel << 17 | 0x1 << 16 |
626                     lADCDataOvFlwLG << 12 | (lADCDataValue3[indLG] & 0xfff);  
627       // Ch. debug
628       //printf(" lADCDataGEO %d  ADCdataHG[%d] %d  ADCdataLG[%d] %d\n", 
629       //  lADCDataGEO,indADC2,lADCDataValue3[indADC2],indLG,lADCDataValue3[indLG]);
630                   
631       indADC2++;
632     }
633     else if(lADCDataGEO==3){ 
634       if(indADC3>=knADCData4){
635          AliWarning(Form(" Problem with digit index %d for ADC2\n", indADC3));
636          return;
637       }
638       Int_t indLG = indADC3+knADCData4;
639       // High gain ADC ch.       
640       if(digit.GetADCValue(0) > 2047) lADCDataOvFlwHG = 1; 
641       lADCDataValue4[indADC3] = digit.GetADCValue(0);    
642       lADCData4[indADC3] = lADCDataGEO << 27 | lADCDataChannel << 17 | 
643                     lADCDataOvFlwHG << 12 | (lADCDataValue4[indADC3] & 0xfff); 
644       // Low gain ADC ch.
645       if(digit.GetADCValue(1) > 2047) lADCDataOvFlwLG = 1; 
646       lADCDataValue4[indLG] = digit.GetADCValue(1); 
647       lADCData4[indLG] = lADCDataGEO << 27 |  lADCDataChannel << 17 | 0x1 << 16 |
648                     lADCDataOvFlwLG << 12 | (lADCDataValue4[indLG] & 0xfff);  
649       // Ch. debug
650       //printf(" lADCDataGEO %d  ADCdataHG[%d] %d  ADCdataLG[%d] %d\n", 
651       //  lADCDataGEO,indADC3,lADCDataValue4[indADC3],indLG,lADCDataValue4[indLG]);
652                   
653       indADC3++;
654     }             
655
656   }
657   //
658   /*for(Int_t i=0;i<2*knADCData1;i++) printf("\t ADCData1[%d] = %x\n",i,lADCData1[i]);
659   for(Int_t i=0;i<2*knADCData2;i++) printf("\t ADCData2[%d] = %x\n",i,lADCData2[i]);
660   for(Int_t i=0;i<2*knADCData3;i++) printf("\t ADCData3[%d] = %x\n",i,lADCData3[i]);
661   for(Int_t i=0;i<2*knADCData4;i++) printf("\t ADCData4[%d] = %x\n",i,lADCData4[i]);*/
662    
663   // End of Block
664   UInt_t lADCEndBlockGEO = 0;
665   // Event counter in ADC EOB -> getting no. of events in run from AliRunLoader
666   // get run loader
667   AliRunLoader* runLoader = fLoader->GetRunLoader(); 
668   UInt_t lADCEndBlockEvCount = runLoader->GetEventNumber();
669   //  
670   lADCEndBlock = lADCEndBlockGEO << 27 | 0x1 << 26 | lADCEndBlockEvCount;
671   //printf("\t AliZDC::Digits2Raw -> ADCEndBlock = %d\n",lADCEndBlock);
672
673   // open the output file
674   TString fileName;
675   fileName.Form("%s",AliDAQ::DdlFileName("ZDC",0)); 
676
677   AliFstream* file = new AliFstream(fileName.Data());
678
679   // write the DDL data header
680   AliRawDataHeaderSim header;
681   header.fSize = sizeof(header) + 
682                  sizeof(lADCHeader1) + sizeof(lADCData1) + sizeof(lADCEndBlock) +
683                  sizeof(lADCHeader2) + sizeof(lADCData2) + sizeof(lADCEndBlock) +
684                  sizeof(lADCHeader3) + sizeof(lADCData3) + sizeof(lADCEndBlock) +
685                  sizeof(lADCHeader4) + sizeof(lADCData4) + sizeof(lADCEndBlock);
686   //
687   /*printf("sizeof header = %d, ADCHeader1 = %d, ADCData1 = %d, ADCEndBlock = %d\n",
688           sizeof(header),sizeof(lADCHeader1),sizeof(lADCData1),sizeof(lADCEndBlock));
689   printf("sizeof header = %d, ADCHeader2 = %d, ADCData2 = %d, ADCEndBlock = %d\n",
690           sizeof(header),sizeof(lADCHeader2),sizeof(lADCData2),sizeof(lADCEndBlock));
691   printf("sizeof header = %d, ADCHeader3 = %d, ADCData3 = %d, ADCEndBlock = %d\n",
692           sizeof(header),sizeof(lADCHeader1),sizeof(lADCData1),sizeof(lADCEndBlock));
693   printf("sizeof header = %d, ADCHeader4 = %d, ADCData4 = %d, ADCEndBlock = %d\n",
694           sizeof(header),sizeof(lADCHeader2),sizeof(lADCData2),sizeof(lADCEndBlock));*/
695   
696   header.SetAttribute(0);  // valid data
697   file->WriteBuffer((char*)(&header), sizeof(header));
698   // write the raw data and close the file
699   file->WriteBuffer((char*) &lADCHeader1,  sizeof (lADCHeader1));
700   file->WriteBuffer((char*) &lADCData1,   sizeof(lADCData1));
701   file->WriteBuffer((char*) &lADCEndBlock, sizeof(lADCEndBlock));
702   file->WriteBuffer((char*) &lADCHeader2,  sizeof (lADCHeader2));
703   file->WriteBuffer((char*) (lADCData2),   sizeof(lADCData2));
704   file->WriteBuffer((char*) &lADCEndBlock, sizeof(lADCEndBlock));
705   file->WriteBuffer((char*) &lADCHeader3,  sizeof (lADCHeader3));
706   file->WriteBuffer((char*) (lADCData3),   sizeof(lADCData3));
707   file->WriteBuffer((char*) &lADCEndBlock, sizeof(lADCEndBlock));
708   file->WriteBuffer((char*) &lADCHeader4,  sizeof (lADCHeader4));
709   file->WriteBuffer((char*) (lADCData4),   sizeof(lADCData4));
710   file->WriteBuffer((char*) &lADCEndBlock, sizeof(lADCEndBlock));
711   delete file;
712
713   // unload the digits
714   fLoader->UnloadDigits();
715 }
716
717 //_____________________________________________________________________________
718 Bool_t AliZDC::Raw2SDigits(AliRawReader* rawReader)
719 {
720   // Convert ZDC raw data to Sdigits
721   const int kNch = 48;
722   AliLoader* loader = (AliRunLoader::Instance())->GetLoader("ZDCLoader");
723   if(!loader) {
724     AliError("no ZDC loader found");
725     return kFALSE;
726   }
727
728   // Create the output digit tree
729   TTree* treeS = loader->TreeS();
730   if(!treeS){
731     loader->MakeTree("S");
732     treeS = loader->TreeS();
733   }
734   //
735   AliZDCSDigit sdigit;
736   AliZDCSDigit* psdigit = &sdigit;
737   const Int_t kBufferSize = 4000;
738   treeS->Branch("ZDC", "AliZDCSDigit",  &psdigit, kBufferSize);
739   //
740   AliZDCRawStream rawStream(rawReader);
741   Int_t sector[2], resADC, rawADC, corrADC, nPheVal;
742   Int_t jcount = 0;
743   while(rawStream.Next()){
744     if(rawStream.IsADCDataWord()){
745       //For the moment only in-time SDigits are foreseen (1st 48 raw values)
746       if(jcount < kNch){ 
747         for(Int_t j=0; j<2; j++) sector[j] = rawStream.GetSector(j);
748         rawADC = rawStream.GetADCValue();
749         resADC = rawStream.GetADCGain();
750         //printf("\t RAw2SDigits raw%d ->  RawADC[%d, %d, %d] read\n",
751         //    jcount, sector[0], sector[1], rawADC);
752         //
753         corrADC = rawADC - Pedestal(sector[0], sector[1], resADC);
754         if(corrADC<0) corrADC=0;
755         nPheVal = ADCch2Phe(sector[0], sector[1], corrADC, resADC);
756         //
757         //printf("\t \t ->  SDigit[%d, %d, %d] created\n",
758         //    sector[0], sector[1], nPheVal);
759         //
760         new(psdigit) AliZDCSDigit(sector, (Float_t) nPheVal, 0.);
761         treeS->Fill();
762         jcount++;
763       }
764     }//IsADCDataWord
765   }//rawStream.Next
766   // write the output tree
767   fLoader->WriteSDigits("OVERWRITE");
768   fLoader->UnloadSDigits();
769    
770   return kTRUE;
771 }
772
773 //_____________________________________________________________________________
774 Int_t AliZDC::Pedestal(Int_t Det, Int_t Quad, Int_t Res) const
775 {
776   // Returns a pedestal for detector det, PM quad, channel with res.
777   //
778   // Getting calibration object for ZDC set
779   Float_t pedValue = 0.;
780   AliCDBManager *man = AliCDBManager::Instance();
781   AliCDBEntry  *entry = man->Get("ZDC/Calib/Pedestals");
782   AliZDCPedestals *calibPed = 0x0;
783   if(!entry) AliFatal("No calibration data loaded!");  
784   else{
785     calibPed = (AliZDCPedestals*) entry->GetObject();
786     //
787     if(!calibPed){
788       printf("\t No calibration object found for ZDC!");
789       return -1;
790     }
791     //
792     Int_t index=0, kNch=24;
793     if(Quad!=5){
794       if(Det==1)        index = Quad+kNch*Res;   // ZN1
795       else if(Det==2)   index = Quad+5+kNch*Res;         // ZP1
796       else if(Det==3)   index = Quad+9+kNch*Res; // ZEM
797       else if(Det==4)   index = Quad+12+kNch*Res; // ZN2
798       else if(Det==5)   index = Quad+17+kNch*Res; // ZP2
799     }
800     else index = (Det-1)/3+22+kNch*Res; // Reference PMs
801     //
802     //
803     Float_t meanPed = calibPed->GetMeanPed(index);
804     Float_t pedWidth = calibPed->GetMeanPedWidth(index);
805     pedValue = gRandom->Gaus(meanPed,pedWidth);
806     //
807     //printf("\t AliZDC::Pedestal - det(%d, %d) - Ped[%d] = %d\n",Det, Quad, index,(Int_t) pedValue); // Chiara debugging!
808   }
809   
810   return (Int_t) pedValue;
811 }
812
813
814 //_____________________________________________________________________________
815 Int_t AliZDC::ADCch2Phe(Int_t Det, Int_t Quad, Int_t ADCVal, Int_t Res) const
816 {
817   // Evaluation of the no. of phe produced
818   Float_t pmGain[6][5];
819   Float_t resADC[2];
820   for(Int_t j = 0; j < 5; j++){
821     pmGain[0][j] = 50000.;
822     pmGain[1][j] = 100000.;
823     pmGain[2][j] = 100000.;
824     pmGain[3][j] = 50000.;
825     pmGain[4][j] = 100000.;
826     pmGain[5][j] = 100000.;
827   }
828   // ADC Caen V965
829   resADC[0] = 0.0000008; // ADC Resolution high gain: 200 fC/adcCh
830   resADC[1] = 0.0000064; // ADC Resolution low gain:  25  fC/adcCh
831   //
832   Int_t nPhe = (Int_t) (ADCVal / (pmGain[Det-1][Quad] * resADC[Res]));
833   //
834   //printf("\t AliZDC::ADCch2Phe -> det(%d, %d) - ADC %d  phe %d\n",Det,Quad,ADCVal,nPhe);
835
836   return nPhe;
837 }
838
839 //______________________________________________________________________
840 void AliZDC::SetTreeAddress(){
841
842   // Set branch address for the Trees.
843   if(fLoader->TreeH() && (fHits == 0x0))
844     fHits   = new TClonesArray("AliZDCHit",1000);
845       
846   AliDetector::SetTreeAddress();
847 }
848
849 //_____________________________________________________________________________
850 AliZDCChMap* AliZDC::GetChMap() const
851 {
852
853   // Getting calibration object for ZDC
854   AliZDCChMap *calibdata = 0x0;
855   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/ChMap");
856   if(!entry) AliFatal("No calibration data loaded!");  
857   else{
858     calibdata = dynamic_cast<AliZDCChMap*> (entry->GetObject());
859     if(!calibdata) AliFatal("Wrong calibration object in calibration  file!");
860   }
861   return calibdata;
862 }