New classes for beam test analysis
[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 // Short description  
20 //
21 /*-- Author: Maxim Volkov (RRC KI)
22              Dmitri Peressounko (RRC KI & SUBATECH)
23              Yuri Kharlov (IHEP & SUBATECH)     */
24
25 //////////////////////////////////////////////////////////////////////////////
26
27 // --- ROOT system ---
28 #include "TClonesArray.h"
29 #include "TFile.h"
30 #include "TTree.h"
31
32 // --- Standard library ---
33 #include <iostream.h>
34 //#include <sys/types.h>
35 #include <sys/stat.h>
36 #include <fcntl.h>
37 #include <unistd.h>
38 #include <ctype.h>
39 #include <stdio.h>
40 #include <stdlib.h>
41 #include <unistd.h>
42 #include <string.h>
43 #include <netinet/in.h>
44
45 // --- AliRoot header files ---
46 #include "AliPHOSDigit.h"
47 #include "AliPHOSConTableDB.h"
48 #include "AliPHOSBeamTestEvent.h"
49 #include "AliPHOSRaw2Digits.h"
50 #include "AliPHOSv1.h"
51 #include "../EVGEN/AliGenBox.h"
52 #include "AliRun.h"
53
54 ClassImp(AliPHOSRaw2Digits)
55   
56   
57 //____________________________________________________________________________ 
58   AliPHOSRaw2Digits::AliPHOSRaw2Digits():TTask() 
59 {
60   fInName="";  
61   fMK1 = 0x0123CDEF ;
62   fMK2 = 0x80708070 ;
63   fMK3 = 0x4321ABCD ;
64   fMK4 = 0x80618061 ;
65   fCKW = 0x4640E400 ;
66   fDebug = kFALSE;             //  Debug flag
67   fIsInitialized = kFALSE ;
68   fTarget[0] = 0 ;
69   fTarget[1] = 0 ;
70   fTarget[2] = 0 ;
71   fDigits = 0 ;
72   fPHOSHeader =0 ;
73   fEvent = 0 ;
74   fctdb = 0;
75 }
76 //____________________________________________________________________________ 
77   AliPHOSRaw2Digits::AliPHOSRaw2Digits(const char * filename):TTask("Default","") 
78 {
79   fInName=filename;
80   TString outname(fInName) ;
81   outname.ToLower() ;
82   outname.ReplaceAll(".fz",".root") ;
83   outname.ReplaceAll(".gz","") ;
84   SetTitle(outname) ;
85
86   fMK1 = 0x0123CDEF ;
87   fMK2 = 0x80708070 ;
88   fMK3 = 0x4321ABCD ;
89   fMK4 = 0x80618061 ;
90   fCKW = 0x4640E400 ;
91   fDebug = kFALSE;             //  Debug flag
92   fIsInitialized = kFALSE ;
93   fTarget[0] = 0 ;
94   fTarget[1] = 0 ;
95   fTarget[2] = 0 ;
96   fDigits = 0 ;
97   fPHOSHeader =0 ;
98   fEvent = 0 ;
99   fctdb = 0;
100 }
101
102 //____________________________________________________________________________ 
103 AliPHOSRaw2Digits::~AliPHOSRaw2Digits()
104 {
105   if(fPHOSHeader)
106     fPHOSHeader->Delete() ;
107   if(fDigits){
108     fDigits->Delete() ;
109     delete fDigits ;
110   }
111   
112 }
113 //____________________________________________________________________________ 
114 void AliPHOSRaw2Digits::Exec(Option_t * option){
115   //This is steering method performing all the conversion
116
117   if(!fIsInitialized) //need initialization
118     if(!Init())       //failed to initialize
119       return ;
120
121   ProcessRawFile() ;
122
123   FinishRun() ;
124
125 //____________________________________________________________________________ 
126 Bool_t AliPHOSRaw2Digits::Init(void){
127   //Create PHOS geometry, sets magnetic field to zero, 
128   //create Generator - to store target position, 
129   //opens out file, creates TreeE and make initialization of contaniers
130
131
132   if(fIsInitialized)
133     return kTRUE;
134
135   //Create PHOS
136   new AliPHOSv1("PHOS","GPS2") ;
137
138   //Set Magnetic field
139   gAlice->SetField(0,2);  
140
141   //Set positin of the virtex
142   AliGenBox *gener = new AliGenBox(1);
143   Float_t ox = fTarget[1]; 
144   Float_t oy = fTarget[2]-460.; 
145   Float_t oz = fTarget[0];
146   gener->SetOrigin(ox, oy, oz);
147
148   //  Create the output file
149   TString outname("") ;
150   if(strstr(GetTitle(),"root")){
151     outname=GetTitle();
152   }
153   else{
154     outname = fInName ;
155     outname.ToLower() ;
156     outname.ReplaceAll(".fz",".root") ;
157   }
158   TFile *rootfile = new TFile(outname,"recreate");
159   rootfile->SetCompressionLevel(2);
160
161   // Create the Root Trees
162   gAlice->MakeTree("E") ;
163
164   //Make container for digits
165   fDigits = new TClonesArray("AliPHOSDigit",1000) ;
166
167   //Fill now TreeE
168   fPHOSHeader = new  AliPHOSBeamTestEvent() ;
169   Int_t splitlevel = 0 ;
170   Int_t bufferSize = 32000 ;    
171   TBranch * headerBranch = gAlice->TreeE()->Branch("AliPHOSBeamTestEvent", 
172                                                    "AliPHOSBeamTestEvent", 
173                                                    &fPHOSHeader,bufferSize,splitlevel);
174   headerBranch->SetName("AliPHOSBeamTestEvent") ;
175
176   fIsInitialized = kTRUE ;
177   return kTRUE ;
178 }
179 //____________________________________________________________________________ 
180 Bool_t AliPHOSRaw2Digits::ProcessRawFile(){
181
182   //Method opens zebra file and reads successively bytes from it,
183   //filling corresponding fields in header and digits.
184
185
186   fStatus= -3 ;
187   //First of all, open file and check if it is a zebra file
188
189   char command[256];
190   sprintf(command,"zcat %s",fInName.Data());
191   FILE *dataFile = popen(command, "r");
192   if (dataFile == NULL) {
193     cout << " Can not open file " << fInName.Data() << endl ;
194     perror(fInName.Data()) ;
195     fStatus = -1 ;
196     return kFALSE ;
197   }
198   printf("Open pipe: %s\n",command);
199
200   // Check if byte ordering is little-endian 
201   UInt_t w = 0x12345678;
202   Int_t swapo = memcmp(&w, "\x78\x56\x34\x12", sizeof(UInt_t)) == 0;
203   if(fDebug)
204     cout << "swapo=" << swapo << endl;
205   
206   
207   UInt_t recBuf[300] ;
208
209   // Read physical record control words 
210   UInt_t nb = 8*sizeof(UInt_t);
211   Int_t n = fread(recBuf, nb, 1, dataFile);
212   if(static_cast<UInt_t>(n) != 1) {
213     if (n < 0 )
214       perror(fInName.Data());
215     else
216       cout << "Could not read physical record control words" << endl;
217     fStatus = -2 ;
218     return kFALSE;
219   }
220   
221   if(fDebug)
222     cout << "recbuf[0] = " << recBuf[0] << endl ;
223   
224   // Check if it is a ZEBRA file and if the words are swapped 
225   UInt_t swapi = 0 ;
226   if (recBuf[0] != fMK1) {
227     Swab4(recBuf, &w, 1);
228     if (w != fMK1) {
229       cout << "Does not look like a ZEBRA file\n";
230       pclose(dataFile) ;
231       fStatus = -2 ;
232       return kFALSE;
233     }
234     swapi=1 ;
235   }
236   
237   if(fDebug){
238     cout << "        w = " << w << endl ;;
239     cout << "    swapi = " << swapi << endl ;;
240   }
241   
242   // Get number of words in physical record 
243   UInt_t  nwphr ;
244   if (swapi)
245     Swab4(&recBuf[4],&nwphr,1);
246   else 
247     nwphr = recBuf[4];
248   nwphr*=2; // 1998 -- Now we have 2 records 150 words each 
249
250
251   //start loop over data  
252   // Read rest of record 
253   nb = (nwphr-8)*sizeof(UInt_t);
254   n = fread(&recBuf[8], nb, 1, dataFile) ;
255   if (static_cast<UInt_t>(n) != 1) {
256     if (n < 0 ){
257       perror(fInName.Data());
258       fStatus = -2 ;
259       return kFALSE;
260     }
261   }
262   nb = nwphr *sizeof(UInt_t);
263
264   UInt_t userVector[16] ;
265   UInt_t zheader[12];    
266   UShort_t pattern ;
267   UShort_t scanning[32] ;
268   UShort_t charge[12];
269   UInt_t scaler[12]; 
270   UShort_t tdc2228[32];
271   
272   //read untill the end of file
273   fEvent=0 ;
274   while(1){
275
276     //    StartNewEvent() ;
277     fDigits->Delete() ;
278     gAlice->SetEvent(fEvent) ;
279           
280     Int_t i ;
281     for(i=0;i<16;i++)
282       userVector[i]=*(recBuf+21+i);
283     if(!swapi)
284       Swab4(userVector, userVector, 16);     
285     fPHOSHeader->SetUserVector(userVector) ;
286     
287     
288     // ZEBRA event header
289     for(i=0;i<12;i++)
290       zheader[i]=*(recBuf+47+i);
291     if(swapi)
292       Swab4(zheader, zheader, 12);
293     fPHOSHeader->SetHeader(zheader) ;
294     
295     // Swap input 
296     if (swapi)
297       Swab4(recBuf, recBuf, nwphr);
298     
299     /* Physical record control words */
300     UInt_t * recptr = recBuf;  //Pointer to current position
301
302     if(recptr[7] != 1) {
303       cout << "Can not handle fast blocks" << endl;
304       fStatus = -2 ;
305       return kFALSE;
306     }    
307     recptr += 8;
308     
309     // Logical record control words   
310     UInt_t lrtyp = recptr[1];
311     if (lrtyp != 3) {
312       cout << "Can not handle logical record type " <<  lrtyp << endl ;
313       fStatus = -2 ;
314       return kFALSE;
315     }
316     
317     recptr += 2;
318     if (recptr[0] != fCKW) {
319       cout <<  "Bad check word" << endl;
320       fStatus = -2 ;
321       return kFALSE;
322     }
323     
324     UInt_t  nwuh = recptr[9];
325     recptr += 10+nwuh;
326     
327     // Bank system words 
328     UInt_t nd = recptr[8];            /* Number of data words */
329     recptr += 10;
330     
331     // Data words 
332     UInt_t evtno = recptr[2];                   /* Event number */
333     
334     if(fDebug)
335       cout << "evtno=" << evtno << endl;
336     
337     UInt_t nh = recptr[4];                   /* Number of header words in data bank */
338     recptr += nh;
339     
340     // Unswap data from VME 
341     if (swapi)
342       Swab4(recptr, recptr, nd-nh-3);
343     
344     // Give buffer to monitor program 
345     //  UInt_t esize = nd-nh-3;
346     //    if (swapo)
347     //       Swab2(recptr, recptr, esize); 
348     // Two byte data are correct after this operation. 
349     //But we're in trouble if the data array contains 4 byte data!
350     
351     // From now on deal with VME data (MSB first, or network byte order).
352     
353     
354     // Fill the event with data from ADCs
355     UChar_t *byteptr=(UChar_t*)recptr;
356     
357     // Trigger bit register  
358     pattern=ntohs(*(UShort_t*)byteptr);
359     fPHOSHeader->SetPattern(pattern) ;
360     byteptr+=sizeof(UShort_t);
361     
362     // Either peak ADCs, 10 modulesX8=80 channels, 
363     //or Kurchatov 64+2 channel ADC 
364     //(the rest of the channels padded with 0xffff) 
365     for(i=0;i<80;i++){
366       Int_t peak = static_cast<Int_t>(ntohs(*(UShort_t*)byteptr));
367       //make digit
368       Int_t absID = fctdb->Raw2AbsId(i) ;
369       if(absID > 0)
370         new((*fDigits)[i])AliPHOSDigit(-1,absID,peak,0.,i) ;
371       if(fDebug){
372         if(peak>(UShort_t)1000)
373           cout << "event=" << fEvent << " peak[" << i << "] = "<<peak << endl;
374       }
375       byteptr+=sizeof(UShort_t);
376     }
377     
378     // Scanning ADCs, 4 modulesX8=32 channels
379     for(i=0;i<32;i++){
380       scanning[i]=ntohs(*(UShort_t*)byteptr);
381       byteptr+=sizeof(UShort_t);
382     }
383     fPHOSHeader->SetScanning(scanning) ;
384     
385     // Charge ADCs, 1 moduleX12=12 channels
386     for(i=0;i<12;i++){
387       charge[i]=ntohs(*(UShort_t*)byteptr);
388       byteptr+=sizeof(UShort_t);
389     }
390     fPHOSHeader->SetCharge(charge) ;
391     
392     // Scalers, 1 moduleX12=12 (4 byte) channels
393     for(i=0;i<12;i++){
394       scaler[i]=ntohl(*(UInt_t*)byteptr);
395       byteptr+=sizeof(UInt_t);
396     }
397     fPHOSHeader->SetScaler(scaler) ;
398     
399     // LeCroy TDC 2228A, 4 moduleX8=32 channels
400     for(i=0;i<8;i++){
401       tdc2228[i]=ntohs(*(UShort_t*)byteptr);
402       byteptr+=sizeof(UShort_t);
403     }
404     fPHOSHeader->SetTDC(tdc2228) ;
405
406     WriteDigits() ;
407     if(fDebug)
408       cout << "event=" << fEvent << " written " << endl;
409  
410     // Read next record 
411     UInt_t nb = nwphr *sizeof(UInt_t);
412     n = fread( recBuf, nb,1,dataFile);
413     if (n < 0 ){
414       perror(fInName);
415       fStatus = -2 ;
416       return kFALSE;
417     }
418     if (static_cast<UInt_t>(n) != 1) {
419       pclose(dataFile) ;
420       fStatus = 1 ;
421       return kTRUE ; //all read
422     }
423     fEvent++ ;
424   }
425   
426   fStatus = 1 ;  
427   return kTRUE ;  
428 }
429 //____________________________________________________________________________ 
430 void AliPHOSRaw2Digits::Swab4(void *from, void *to, size_t nwords){
431   // The function swaps 4 bytes: byte#3<-->byte#0, byte#2<-->byte#1 
432   register char *pf=static_cast<char*>(from) ;
433   register char *pt=static_cast<char*>(to) ;
434   register char c;
435   while (nwords-- > 0 ) {
436     c = pf[0];
437     pt[0] = pf[3];
438     pt[3] = c;
439     c = pf[1];
440     pt[1] = pf[2];
441     pt[2] = c;
442     pf += 4;
443     pt += 4;
444   }
445 }
446
447 //____________________________________________________________________________ 
448 void AliPHOSRaw2Digits::Swab2(void *from, void *to, size_t nwords)
449 { //The function swaps 2x2 bytes: byte#0<-->byte#1, byte#2<-->byte#3 
450   register char *pf=static_cast<char*>(from) ;
451   register char *pt=static_cast<char*>(to);
452   register char c;   
453   while (nwords-- > 0 ) {
454     c = pf[0];
455     pt[0] = pf[1];
456     pt[1] = c;
457     c = pf[2];
458     pt[2] = pf[3];
459     pt[3] = c;
460     pf += 4;
461     pt += 4;
462   }
463 }
464
465 //____________________________________________________________________________ 
466 void AliPHOSRaw2Digits::FinishRun(){
467   //Write geometry and header tree
468   gAlice->Write(0,TObject::kOverwrite);
469   gAlice->TreeE()->Write(0,TObject::kOverwrite);
470   
471 }
472 //____________________________________________________________________________ 
473 void AliPHOSRaw2Digits::WriteDigits(void){
474   //In this method we create TreeD, write digits and Raw2Digits to it
475   // and write Header to TreeE. Finally we write TreeD to root file 
476   
477   //Start from Digits
478   fDigits->Sort() ;
479   fDigits->Expand(fDigits->GetEntriesFast()) ;
480   for(Int_t i=0;i<fDigits->GetEntriesFast(); i++)
481     static_cast<AliPHOSDigit*>(fDigits->At(i))->SetIndexInList(i) ;
482
483   char hname[30];
484   sprintf(hname,"TreeD%d",fEvent);
485   TTree * treeD = new TTree(hname,"Digits");
486   //treeD->Write(0,TObject::kOverwrite);
487   
488   // -- create Digits branch
489   Int_t bufferSize = 32000 ;    
490   TBranch * digitsBranch = treeD->Branch("PHOS",&fDigits,bufferSize);
491   digitsBranch->SetTitle("Default");
492   
493   // -- Create Digitizer branch
494   Int_t splitlevel = 0 ;
495   const AliPHOSRaw2Digits * d = this ;
496   TBranch * digitizerBranch = treeD->Branch("AliPHOSRaw2Digits", 
497                                             "AliPHOSRaw2Digits", &d,bufferSize,splitlevel); 
498   digitizerBranch->SetTitle("Default");
499   
500   digitsBranch->Fill() ;
501   digitizerBranch->Fill() ; 
502   treeD->Write(0,TObject::kOverwrite);
503  
504   delete treeD ;
505
506   //Write header
507   gAlice->TreeE()->Fill();
508 }
509 //____________________________________________________________________________ 
510 void AliPHOSRaw2Digits::Print(Option_t * option)const{
511   
512   cout << "----------AliPHOSRaw2Digits----------" << endl ;
513   cout << "Input stream " << endl ;
514   cout << "Current input  File: " << fInName.Data() << endl ;
515   cout << "Current output File: " << GetTitle() << endl ;
516   cout << "Events processes in the last file " << fEvent << endl ;
517   cout << "Input file status " ;
518   switch (fStatus){
519   case 0: cout << "`Have not processed yet' " << endl ;
520     break ;
521   case 1: cout << "`Processed normally' " << endl ;
522     break ;
523   case -1: cout << "`File not found'" << endl ;
524     break ;
525   case -2: cout << "`Error in reading' " << endl ;
526     break ;
527   case -3: cout << "'Interupted'" << endl ;
528   default: ;
529   }
530   cout << "Connection table: " ;
531   if(fctdb)
532     cout << fctdb->GetName() << "  " << fctdb->GetTitle() << endl ;
533   else
534     cout << " no DB " << endl ;
535   
536 }