removed #include <netinet/in.h>
[u/mrichter/AliRoot.git] / PHOS / AliPHOSRaw2Digits.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 // Class designed to convert raw data to aliroot format. We assume, that
20 // prototype is situated in the center of 3 PHOS module and convert prototype
21 // outpur to AliPHOSDigits. In addition, we fill branch of TreeE with 
22 // AliPHOSBeamTestEvent, contaning description of event(triggers etc).
23 // Note, that one byte per channel in raw data is transvormed to class 
24 // AliPHOSDigit, so finale zise increase ~ 100 times. So output can be split 
25 // into peases of reasonable size: each file can not contain more than 
26 // fMaxPerFile: if there are still events in raw file, then new directory 
27 // is created and header+digits files are written to it.
28 // 
29 // Use Case:
30 //   AliPHOSRaw2Digits * r = new AliPHOSRaw2Digits("path/input.file") ;
31 //                   //note, that it can be gzipped file!
32 //   //Set position of the target in the given run.
33 //   //Z axis along beam direction, from target to prototype (0-surface of prototype)
34 //   //X axis along columns of prototype (0-center of prototype)
35 //   //Y axis along raws of prototype    (0-center of prototype)
36 //   Double_t pos[3]={0,0,-452.} ;
37 //   r->SetTargetPosition(pos) ;
38 //   //Read/create connection Table:
39 //   TFile f("ConTableDB.root") ;
40 //   AliPHOSConTableDB * cdb = f.Get("AliPHOSConTableDB") ;
41 //   f.Close() ;
42 //   r->SetConTableDB(cdb) ;
43 //   r->ExecuteTask() ;
44 //
45 // As a result files galice.root and PHOS.Digits.root should be produced in 
46 // current dir, and, possibly, dirs 1,2,3... each with galice.root and PHOS.Digits.root,
47 // where the rest of data are written. 
48 //
49 /*-- Author: Maxim Volkov (RRC KI)
50              Dmitri Peressounko (RRC KI & SUBATECH)
51              Yuri Kharlov (IHEP & SUBATECH)     */
52
53 //////////////////////////////////////////////////////////////////////////////
54
55 // --- ROOT system ---
56 #include "TClonesArray.h"
57 #include "TFile.h"
58 #include "TTree.h"
59 #include "TSystem.h"
60 #include "Bytes.h"
61
62 // --- Standard library ---
63 #include <sys/stat.h>
64 #include <fcntl.h>
65 #include <unistd.h>
66 #include <ctype.h>
67 #include <stdio.h>
68 #include <stdlib.h>
69 #include <unistd.h>
70 #include <string.h>
71
72 // --- AliRoot header files ---
73 #include "AliPHOSDigit.h"
74 #include "AliPHOSConTableDB.h"
75 #include "AliPHOSBeamTestEvent.h"
76 #include "AliPHOSRaw2Digits.h"
77 #include "AliPHOSv1.h"
78 #include "../EVGEN/AliGenBox.h"
79 #include "AliRun.h"
80
81 ClassImp(AliPHOSRaw2Digits)
82   
83   
84 //____________________________________________________________________________ 
85   AliPHOSRaw2Digits::AliPHOSRaw2Digits():TTask() 
86 {
87   fInName="";  
88   fMK1 = 0x0123CDEF ;
89   fMK2 = 0x80708070 ;
90   fMK3 = 0x4321ABCD ;
91   fMK4 = 0x80618061 ;
92   fCKW = 0x4640E400 ;
93   fDebug = kFALSE;             //  Debug flag
94   fToSplit = kFALSE ;
95   fIsInitialized = kFALSE ;
96   fTarget[0] = 0 ;
97   fTarget[1] = 0 ;
98   fTarget[2] = 0 ;
99   fDigits = 0 ;
100   fHeaderFile = 0 ;
101   fDigitsFile = 0 ;
102   fPHOSHeader =0 ;
103   fMaxPerFile = 20000 ;
104   fEvent = 0 ;
105   fctdb = 0;
106 }
107 //____________________________________________________________________________ 
108   AliPHOSRaw2Digits::AliPHOSRaw2Digits(const char * filename,Bool_t toSplit):TTask("Default","") 
109 {
110   fInName=filename;
111   fToSplit = toSplit ;
112   TString outname("") ;
113   if(fToSplit)
114     outname = "galice.root" ;
115   else{
116     outname =fInName ;
117     outname.ToLower() ;
118     outname.ReplaceAll(".fz",".root") ;
119     outname.ReplaceAll(".gz","") ;
120   }
121   SetTitle(outname) ;
122
123   fMK1 = 0x0123CDEF ;
124   fMK2 = 0x80708070 ;
125   fMK3 = 0x4321ABCD ;
126   fMK4 = 0x80618061 ;
127   fCKW = 0x4640E400 ;
128   fDebug = kFALSE;             //  Debug flag
129   fIsInitialized = kFALSE ;
130   fTarget[0] = 0 ;
131   fTarget[1] = 0 ;
132   fTarget[2] = 0 ;
133   fDigits = 0 ;
134   fPHOSHeader =0 ;
135   fHeaderFile = 0 ;
136   fDigitsFile = 0 ;
137   fMaxPerFile = 20000 ;
138   fEvent = 0 ;
139   fctdb = 0;
140 }
141
142 //____________________________________________________________________________ 
143 AliPHOSRaw2Digits::~AliPHOSRaw2Digits()
144 {
145   if(fPHOSHeader)
146     fPHOSHeader->Delete() ;
147   if(fDigits){
148     fDigits->Delete() ;
149     delete fDigits ;
150   }
151   
152 }
153 //____________________________________________________________________________ 
154 void AliPHOSRaw2Digits::Exec(Option_t * option){
155   //This is steering method performing all the conversion
156
157   if(!fIsInitialized) //need initialization
158     if(!Init())       //failed to initialize
159       return ;
160
161   ProcessRawFile() ;
162
163   FinishRun() ;
164
165 //____________________________________________________________________________ 
166 Bool_t AliPHOSRaw2Digits::Init(void){
167   //Makes initialization of contaniers
168
169   if(fIsInitialized)
170     return kTRUE;
171
172   //Make container for digits
173   fDigits = new TClonesArray("AliPHOSDigit",1000) ;
174   fPHOSHeader = new  AliPHOSBeamTestEvent() ;
175   fIsInitialized = kTRUE ;
176   return StartRootFiles() ;
177
178 }
179 //____________________________________________________________________________ 
180 Bool_t AliPHOSRaw2Digits::StartRootFiles(void ){
181   //Create PHOS geometry, sets magnetic field to zero, 
182   //create Generator - to store target position, 
183   //opens out file, creates TreeE 
184
185   //create gAlice if nececcary
186   if(!gAlice)
187     new AliRun("gAlice","The ALICE Off-line Simulation Framework") ;
188
189   //Create PHOS
190   if(!gAlice->GetModule("PHOS"))
191     new AliPHOSv1("PHOS","GPS2") ;
192
193   //Set Magnetic field
194   gAlice->SetField(0,2);  
195
196   //Set positin of the virtex
197   AliGenerator * gener = gAlice->Generator() ; 
198   if(!gener)    
199     gener = new AliGenBox(1);
200   Float_t ox = fTarget[1]; 
201   Float_t oy = fTarget[2]-460.; 
202   Float_t oz = fTarget[0];
203   gener->SetOrigin(ox, oy, oz);
204
205   //make directory 
206   Int_t nRootFile = (fEvent+1)/fMaxPerFile ;    
207   if(nRootFile){
208     char dname[20];
209     sprintf(dname,"%d",nRootFile) ;
210     if(gSystem->AccessPathName(dname)) //strange return: 0 if exists
211       if(gSystem->MakeDirectory(dname)!=0)
212         Fatal("StartRootFiles","Can not make directory %s \n",dname) ;
213     
214     if(!gSystem->ChangeDirectory(dname))
215       Fatal("StartRootFiles","Can not cd to %s\n",dname) ;
216   }
217
218   //  Create the output file
219   TString outname("") ;
220   if(strstr(GetTitle(),"root")){
221     outname=GetTitle();
222   }
223   else{
224     outname = fInName ;
225     outname.ToLower() ;
226     outname.ReplaceAll(".fz",".root") ;
227   }
228
229   fHeaderFile = new TFile(outname,"recreate");
230   fHeaderFile->SetCompressionLevel(2);
231   
232   // Create the Root Trees
233   gAlice->MakeTree("E") ;
234   
235   //Fill now TreeE
236   Int_t splitlevel = 0 ;
237   Int_t bufferSize = 32000 ;    
238   TBranch * headerBranch = gAlice->TreeE()->Branch("AliPHOSBeamTestEvent", 
239                                                    "AliPHOSBeamTestEvent", 
240                                                    &fPHOSHeader,bufferSize,splitlevel);
241   headerBranch->SetName("AliPHOSBeamTestEvent") ;
242
243   if(fToSplit){
244     fDigitsFile = new TFile("PHOS.Digits.root","recreate") ;
245     fDigitsFile->SetCompressionLevel(2) ;
246   }
247   return kTRUE ;
248 }
249 //____________________________________________________________________________ 
250 Bool_t AliPHOSRaw2Digits::CloseRootFiles(void ){
251   //cleans everething to start next root file
252
253   delete gAlice ;
254   
255   if(fHeaderFile){
256     fHeaderFile->Close() ;
257     delete fHeaderFile ;
258     fHeaderFile = 0;
259   }   
260         
261   if(fDigitsFile){
262     fDigitsFile->Close() ;
263     delete fDigitsFile ;
264     fDigitsFile = 0 ;
265   }
266   
267   Int_t nRootFile = (fEvent-1)/fMaxPerFile ;    
268   if(nRootFile){
269     if(!gSystem->ChangeDirectory("../")){
270      Fatal("CloseRootFile","Can not return to initial dir \n") ;      
271      return kFALSE ;
272     }
273   }
274   return kTRUE ;
275 }
276 //____________________________________________________________________________ 
277 Bool_t AliPHOSRaw2Digits::ProcessRawFile(){
278   
279   //Method opens zebra file and reads successively bytes from it,
280   //filling corresponding fields in header and digits.
281   
282   
283   fStatus= -3 ;
284   //First of all, open file and check if it is a zebra file
285   
286   char command[256];
287   sprintf(command,"zcat %s",fInName.Data());
288   FILE *dataFile = popen(command, "r");
289   if (dataFile == NULL) {
290     Warning("ProcessRawFile", " Cannot open file %s\n", fInName.Data() ) ;
291     perror(fInName.Data()) ;
292     fStatus = -1 ;
293     return kFALSE ;
294   }
295
296   // Check if byte ordering is little-endian 
297   UInt_t w = 0x12345678;
298   Int_t swapo = memcmp(&w, "\x78\x56\x34\x12", sizeof(UInt_t)) == 0;
299   if(fDebug)
300     Info("ProcessRawFile", "swapo=%f\n", swapo ) ;
301   
302   
303   UInt_t recBuf[300] ;
304
305   // Read physical record control words 
306   UInt_t nb = 8*sizeof(UInt_t);
307   Int_t n = fread(recBuf, nb, 1, dataFile);
308   if(static_cast<UInt_t>(n) != 1) {
309     if (n < 0 )
310       perror(fInName.Data());
311     else
312       Error("ProcessRawFile", "Could not read physical record control words" ) ;
313     fStatus = -2 ;
314     return kFALSE;
315   }
316   
317   if(fDebug)
318     Info("ProcessRawFile", "recbuf[0] = %d\n", recBuf[0] );
319   
320   // Check if it is a ZEBRA file and if the words are swapped 
321   UInt_t swapi = 0 ;
322   if (recBuf[0] != fMK1) {
323     Swab4(recBuf, &w, 1);
324     if (w != fMK1) {
325       Error("ProcessRawFile", "Does not look like a ZEBRA file\n" ) ;
326       pclose(dataFile) ;
327       fStatus = -2 ;
328       return kFALSE;
329     }
330     swapi=1 ;
331   }
332   
333   if(fDebug){
334     TString message ; 
335     message  = "        w = %f\n" ; 
336     message += "    swapi = %f\n" ; 
337     Info("ProcessRawFile", message.Data(), w, swapi ) ; 
338   }
339   
340   // Get number of words in physical record 
341   UInt_t  nwphr ;
342   if (swapi)
343     Swab4(&recBuf[4],&nwphr,1);
344   else 
345     nwphr = recBuf[4];
346   nwphr*=2; // 1998 -- Now we have 2 records 150 words each 
347
348
349   //start loop over data  
350   // Read rest of record 
351   nb = (nwphr-8)*sizeof(UInt_t);
352   n = fread(&recBuf[8], nb, 1, dataFile) ;
353   if (static_cast<UInt_t>(n) != 1) {
354     if (n < 0 ){
355       perror(fInName.Data());
356       fStatus = -2 ;
357       return kFALSE;
358     }
359   }
360   nb = nwphr *sizeof(UInt_t);
361
362   UInt_t userVector[16] ;
363   UInt_t zheader[12];    
364   UShort_t pattern ;
365   UShort_t scanning[32] ;
366   UShort_t charge[12];
367   UInt_t scaler[12]; 
368   UShort_t tdc2228[32];
369   
370   //read untill the end of file
371   fEvent=0 ;
372   while(1){
373
374     //    StartNewEvent() ;
375     fDigits->Delete() ;
376     if((fEvent%fMaxPerFile == 0) && fEvent ){
377       CloseRootFiles() ;
378       StartRootFiles() ;
379     }
380     gAlice->SetEvent(fEvent%fMaxPerFile) ;
381
382     //Set Beam Energy
383     fPHOSHeader->SetBeamEnergy(fBeamEnergy) ;
384           
385     Int_t i ;
386     for(i=0;i<16;i++)
387       userVector[i]=*(recBuf+21+i);
388     if(!swapi)
389       Swab4(userVector, userVector, 16);     
390     fPHOSHeader->SetUserVector(userVector) ;
391     
392     
393     // ZEBRA event header
394     for(i=0;i<12;i++)
395       zheader[i]=*(recBuf+47+i);
396     if(swapi)
397       Swab4(zheader, zheader, 12);
398     fPHOSHeader->SetHeader(zheader) ;
399     
400     // Swap input 
401     if (swapi)
402       Swab4(recBuf, recBuf, nwphr);
403     
404     /* Physical record control words */
405     UInt_t * recptr = recBuf;  //Pointer to current position
406
407     if(recptr[7] != 1) {
408       Error("ProcessRawFile", "Cannot handle fast blocks" ) ; 
409       fStatus = -2 ;
410       return kFALSE;
411     }    
412     recptr += 8;
413     
414     // Logical record control words   
415     UInt_t lrtyp = recptr[1];
416     if (lrtyp != 3) {
417       Error("ProcessRawFile", "Can not handle logical record type %d", lrtyp ) ;
418       fStatus = -2 ;
419       return kFALSE;
420     }
421     
422     recptr += 2;
423     if (recptr[0] != fCKW) {
424       Error("ProcessRawFile", "Bad check word" ) ;
425       fStatus = -2 ;
426       return kFALSE;
427     }
428     
429     UInt_t  nwuh = recptr[9];
430     recptr += 10+nwuh;
431     
432     // Bank system words 
433     UInt_t nd = recptr[8];            /* Number of data words */
434     recptr += 10;
435     
436     // Data words 
437     UInt_t evtno = recptr[2];                   /* Event number */
438     
439     if(fDebug)
440        Info("ProcessRawFile", "evtno= %d", evtno);
441     
442     UInt_t nh = recptr[4];                   /* Number of header words in data bank */
443     recptr += nh;
444     
445     // Unswap data from VME 
446     if (swapi)
447       Swab4(recptr, recptr, nd-nh-3);
448     
449     // Give buffer to monitor program 
450     //  UInt_t esize = nd-nh-3;
451     //    if (swapo)
452     //       Swab2(recptr, recptr, esize); 
453     // Two byte data are correct after this operation. 
454     //But we're in trouble if the data array contains 4 byte data!
455     
456     // From now on deal with VME data (MSB first, or network byte order).
457     
458     
459     // Fill the event with data from ADCs
460     UChar_t *byteptr=(UChar_t*)recptr;
461     
462     // Trigger bit register  
463     pattern=net2host(*(UShort_t*)byteptr);
464     fPHOSHeader->SetPattern(pattern) ;
465     byteptr+=sizeof(UShort_t);
466     
467     // Either peak ADCs, 10 modulesX8=80 channels, 
468     //or Kurchatov 64+2 channel ADC 
469     //(the rest of the channels padded with 0xffff) 
470     for(i=0;i<80;i++){
471       Int_t peak = static_cast<Int_t>(net2host(*(UShort_t*)byteptr));
472       //make digit
473       Int_t absID = fctdb->Raw2AbsId(i) ;
474       if(absID > 0)
475         new((*fDigits)[i])AliPHOSDigit(-1,absID,peak,0.,i) ;
476       if(fDebug){
477         if(peak>(UShort_t)1000)
478           Info("ProcessRawFile", "event= %d peak[%d] = %f", fEvent, i, peak);
479       }
480       byteptr+=sizeof(UShort_t);
481     }
482     
483     // Scanning ADCs, 4 modulesX8=32 channels
484     for(i=0;i<32;i++){
485       scanning[i]=net2host(*(UShort_t*)byteptr);
486       byteptr+=sizeof(UShort_t);
487     }
488     fPHOSHeader->SetScanning(scanning) ;
489     
490     // Charge ADCs, 1 moduleX12=12 channels
491     for(i=0;i<12;i++){
492       charge[i]=net2host(*(UShort_t*)byteptr);
493       byteptr+=sizeof(UShort_t);
494     }
495     fPHOSHeader->SetCharge(charge) ;
496     
497     // Scalers, 1 moduleX12=12 (4 byte) channels
498     for(i=0;i<12;i++){
499       scaler[i]=net2host(*(UInt_t*)byteptr);
500       byteptr+=sizeof(UInt_t);
501     }
502     fPHOSHeader->SetScaler(scaler) ;
503     
504     // LeCroy TDC 2228A, 4 moduleX8=32 channels
505     for(i=0;i<8;i++){
506       tdc2228[i]=net2host(*(UShort_t*)byteptr);
507       byteptr+=sizeof(UShort_t);
508     }
509     fPHOSHeader->SetTDC(tdc2228) ;
510
511     WriteDigits() ;
512     if(fDebug)
513       Info("ProcessRawFile", "event= %d written", fEvent) ;
514  
515     // Read next record 
516     UInt_t nb = nwphr *sizeof(UInt_t);
517     n = fread( recBuf, nb,1,dataFile);
518     if (n < 0 ){
519       perror(fInName);
520       fStatus = -2 ;
521       return kFALSE;
522     }
523     if (static_cast<UInt_t>(n) != 1) {
524       pclose(dataFile) ;
525       fStatus = 1 ;
526       return kTRUE ; //all read
527     }
528     fEvent++ ;
529   }
530   
531   fStatus = 1 ;  
532   return kTRUE ;  
533 }
534 //____________________________________________________________________________ 
535 void AliPHOSRaw2Digits::Swab4(void *from, void *to, size_t nwords){
536   // The function swaps 4 bytes: byte#3<-->byte#0, byte#2<-->byte#1 
537   register char *pf=static_cast<char*>(from) ;
538   register char *pt=static_cast<char*>(to) ;
539   register char c;
540   while (nwords-- > 0 ) {
541     c = pf[0];
542     pt[0] = pf[3];
543     pt[3] = c;
544     c = pf[1];
545     pt[1] = pf[2];
546     pt[2] = c;
547     pf += 4;
548     pt += 4;
549   }
550 }
551
552 //____________________________________________________________________________ 
553 void AliPHOSRaw2Digits::Swab2(void *from, void *to, size_t nwords)
554 { //The function swaps 2x2 bytes: byte#0<-->byte#1, byte#2<-->byte#3 
555   register char *pf=static_cast<char*>(from) ;
556   register char *pt=static_cast<char*>(to);
557   register char c;   
558   while (nwords-- > 0 ) {
559     c = pf[0];
560     pt[0] = pf[1];
561     pt[1] = c;
562     c = pf[2];
563     pt[2] = pf[3];
564     pt[3] = c;
565     pf += 4;
566     pt += 4;
567   }
568 }
569
570 //____________________________________________________________________________ 
571 void AliPHOSRaw2Digits::FinishRun(){
572   //Write geometry and header tree
573   gAlice->Write(0,TObject::kOverwrite);
574   gAlice->TreeE()->Write(0,TObject::kOverwrite);
575   
576 }
577 //____________________________________________________________________________ 
578 void AliPHOSRaw2Digits::WriteDigits(void){
579   //In this method we create TreeD, write digits and Raw2Digits to it
580   // and write Header to TreeE. Finally we write TreeD to root file 
581   
582   //Start from Digits
583   fDigits->Sort() ;
584   fDigits->Expand(fDigits->GetEntriesFast()) ;
585   for(Int_t i=0;i<fDigits->GetEntriesFast(); i++)
586     static_cast<AliPHOSDigit*>(fDigits->At(i))->SetIndexInList(i) ;
587
588   char hname[30];
589   sprintf(hname,"TreeD%d",fEvent%fMaxPerFile);
590   TTree * treeD = new TTree(hname,"Digits");
591   //treeD->Write(0,TObject::kOverwrite);
592   
593   // -- create Digits branch
594   Int_t bufferSize = 32000 ;    
595   TBranch * digitsBranch = treeD->Branch("PHOS",&fDigits,bufferSize);
596   digitsBranch->SetTitle("Default");
597   
598   // -- Create Digitizer branch
599   Int_t splitlevel = 0 ;
600   const AliPHOSRaw2Digits * d = this ;
601   TBranch * digitizerBranch = treeD->Branch("AliPHOSRaw2Digits", 
602                                             "AliPHOSRaw2Digits", &d,bufferSize,splitlevel); 
603   digitizerBranch->SetTitle("Default");
604
605   if(fDigitsFile)
606     fDigitsFile->cd() ;
607   digitsBranch->Fill() ;
608   digitizerBranch->Fill() ; 
609   treeD->Write(0,TObject::kOverwrite);
610  
611   delete treeD ;
612
613   //Write header
614   fHeaderFile->cd() ;
615   gAlice->TreeE()->Fill();
616 }
617 //____________________________________________________________________________ 
618 void AliPHOSRaw2Digits::Print(Option_t * option)const{
619   
620
621   printf("----------AliPHOSRaw2Digits---------- \n") ;
622   printf("Current input  File: %s\n",fInName.Data()) ; 
623   printf("Current output File: %s\n", GetTitle()); 
624   printf("Events processes in the last file %d\n",fEvent) ; 
625   printf("Input file status\n") ;
626   switch (fStatus){
627   case 0: printf("`Have not processed yet'\n") ;
628     break ;
629   case 1: printf("`Processed normally'\n") ;
630     break ;
631   case -1: printf("`File not found'\n") ;
632     break ;
633   case -2: printf("`Error in reading'\n") ;
634     break ;
635   case -3: printf("'Interupted'\n") ;
636   default: ;
637   }
638   printf("Connection table: " );
639   if(fctdb)
640     printf("%s %s \n",fctdb->GetName(), fctdb->GetTitle() ) ; 
641   else
642     printf(" no DB \n" );
643
644 }