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