cal/src/PamVMCCaloDig.cxx_my

#include "PamVMCCaloDig.h"
#include <TTree.h>
#include "CalibCalPedEvent.h"

ClassImp(PamVMCCaloDig)

extern "C"{
  short crc(short, short);
};


void PamVMCCaloDig::LoadCalib(){

  cout<<"Loading Tracker Calibrations..."<<endl; 


  fhc =(TClonesArray*)fhitscolmap.GetValue("CAST");
  if (!fhc) cout <<"!!!WARNING Hit Collection of CAL not found by PamVMCCaloDig!!!!"<<endl;

  Int_t GivenCaloCalib = 0; //@@@ should be given as input par @@@

  Int_t time=3;
  Int_t type=101;

  fdberr = fsql->Query_GL_PARAM(time,type); 

  if(fdberr<0){
    cout<<"No such record in DB for CAL: time="<<time<<" type="<<type<<endl;
    ThrowCalFileWarning("CAL");

  } else {

    fquery.str("");
    fquery << fsql->GetPAR()->PATH.Data() << "/";
    fquery << fsql->GetPAR()->NAME.Data();

    ThrowCalFileUsage("CAL",fquery.str().c_str());

    FILE *f;
    f =  fopen(fquery.str().c_str(),"rb");

    if(f==NULL) ThrowCalFileWarning("CAL");
    else{
      memset(fCalomip,0,4224*sizeof(fCalomip[0][0][0]));
      
      for (Int_t m = 0; m < 2 ; m++ ){
        for (Int_t k = 0; k < 22; k++ ){
          for (Int_t l = 0; l < 96; l++ ){
            fread(&fCalomip[m][k][l],sizeof(fCalomip[m][k][l]),1,f);
          };
        };
      };
      fclose(f);
      
      // determine which calibration has to be used 
      // and load it for each section
      

      UInt_t idcalib, calibno;
      UInt_t utime = 0;
      TString calname;
      for (UInt_t s=0; s<4; s++){
        
        // clear calo calib variables for section s
        
        ClearCaloCalib(s);
        
        if ( GivenCaloCalib ){
          
          // a time has been given as input on the command line 
          // so retrieve the calibration that preceed that time
          
          fsql->Query_GL_CALO_CALIB(GivenCaloCalib,utime,s);
            
          calibno = fsql->GetCaloCalib()->EV_ROOT;
          idcalib = fsql->GetCaloCalib()->ID_ROOT_L0;
          
          // determine path and name and entry of the calibration file
          
          printf("\n");
          printf(" ** SECTION %i **\n",s);
          
          fsql->Query_GL_ROOT(idcalib);
          
          fquery.str("");
          fquery << fsql->GetROOT()->PATH.Data() << "/";
          fquery << fsql->GetROOT()->NAME.Data();
          
          calname = (TString)fquery.str().c_str();
          
          printf("\n Section %i : using  file %s calibration at entry %i: \n",s,calname.Data(),calibno);
          
        } else {

          fdberr = 0;
          fdberr = fsql->Query_GL_PARAM(1,104); 

          fquery.str("");
          fquery << fsql->GetPAR()->PATH.Data() << "/";
          fquery << fsql->GetPAR()->NAME.Data(); 
          
          printf("\n Section %i : using default calorimeter calibration: \n %s \n",s,fquery.str().c_str());

          calname = (TString)fquery.str().c_str();
          calibno = s;
          
        };
        
        // load calibration variables in memory
        
        CaloLoadCalib(s,calname,calibno);
        
      };
      //
      // at this point we have in memory the calorimeter calibration 
      // and we can save it to disk in the correct format and use it to digitize the data
      //

      DigitizeCALOCALIB();
    }
  }        
}


void PamVMCCaloDig::ClearCaloCalib(Int_t s){
  
  fcstwerr[s] = 0;
  fcperror[s] = 0.;
  for ( Int_t d=0 ; d<11 ;d++  ){
    Int_t pre = -1;
    for ( Int_t j=0; j<96 ;j++){
      if ( j%16 == 0 ) pre++;
      fcalped[s][d][j] = 0.;
      fcstwerr[s] = 0.;
      fcperror[s] = 0.;
      fcalgood[s][d][j] = 0.;
      fcalthr[s][d][pre] = 0.;
      fcalrms[s][d][j] = 0.;
      fcalbase[s][d][pre] = 0.;
      fcalvar[s][d][pre] = 0.;
    };
  };
  return;
}

Int_t PamVMCCaloDig::CaloLoadCalib(Int_t s,TString calname, UInt_t calibno){
 
  UInt_t e = 0;
  switch(s){
  case 0: e = 0;
    break;
  case 1: e = 2;
    break;
  case 2: e = 3;
    break;
  case 3: e = 1;
    break;
  default:
    break;
  }
  fcfile.open(calname.Data());
  ThrowCalFileUsage("CAL",calname.Data());
  if ( !fcfile ){
    ThrowCalFileWarning("CAL");
    return(-107);
  };
  fcfile.close();
  
  fcrfile = new TFile(calname.Data());
  if ( !fcrfile ){
    ThrowCalFileWarning("CAL");
    return(-108);
  }
  TTree *tr = (TTree*)fcrfile->Get("CalibCalPed");
  if ( !tr ){
    cout<<"!!!WARNING Tree CalibCalPed in file "<<calname.Data()
        <<" was NOT found!!!"<<endl;
    return(-109);
  }
  
  TBranch *calo = tr->GetBranch("CalibCalPed");
  
  pamela::CalibCalPedEvent *ce = 0;
  tr->SetBranchAddress("CalibCalPed", &ce);
 
  Long64_t ncalibs = calo->GetEntries();
 
  if ( !ncalibs ){
    cout<<"!!!WARNING No entries in calibration file!!!"<<endl; 
    return(-110); 
  }
  
  calo->GetEntry(calibno);
  
  if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
    fcstwerr[s] = ce->cstwerr[s];
    fcperror[s] = ce->cperror[s];
    for ( Int_t d=0 ; d<11 ;d++  ){
      Int_t pre = -1;
      for ( Int_t j=0; j<96 ;j++){
        if ( j%16 == 0 ) pre++;
        fcalped[s][d][j] = ce->calped[e][d][j];
        fcalgood[s][d][j] = ce->calgood[e][d][j];
        fcalthr[s][d][pre] = ce->calthr[e][d][pre];
        fcalrms[s][d][j] = ce->calrms[e][d][j];
        fcalbase[s][d][pre] = ce->calbase[e][d][pre];
        fcalvar[s][d][pre] = ce->calvar[e][d][pre];
      };
    };
  } else {
    printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
    fcrfile->Close();
    return(-111);
  };
  fcrfile->Close();
  return(0);
}

void PamVMCCaloDig::DigitizeCALOCALIB() {

  
  // Header of the four sections
  
  fSecCalo[0] = 0xAA00; // XE
  fSecCalo[1] = 0xB100; // XO
  fSecCalo[2] = 0xB600; // YE
  fSecCalo[3] = 0xAD00; // YO
  
  // length of the data is 0x1215 
  
  fSecCaloLength[0] = 0x1215; // XE
  fSecCaloLength[1] = 0x1215; // XO
  fSecCaloLength[2] = 0x1215; // YE
  fSecCaloLength[3] = 0x1215; // YO 

  Short_t CRC;

  for(Int_t sec=0; sec < 4; sec++){
    //
    // sec =    0 -> XE      1 -> XO        2-> YE         3 -> YO
    //
    fData.clear();

    fData.push_back(fSecCalo[sec]);      
    fData.push_back(fSecCaloLength[sec]);
    
    if ( sec==1 ){
      FillCalPedBack(&fData,sec);
      FillCalThrBack(&fData,sec);
      FillCalRmsBack(&fData,sec);
      FillCalBaseVarBack(&fData,sec);
    } else {
      FillCalPedNorm(&fData,sec);
      FillCalThrNorm(&fData,sec);
      FillCalRmsNorm(&fData,sec);
      FillCalBaseVarNorm(&fData,sec);
    }
    //CRC
    
    fData.push_back(0x0000);
    CRC = 0;

  
    USBuffer::const_iterator p = fData.begin();
    while( p!= fData.end() ){
      CRC=crc(CRC,*p);
      p++;
    }
    fData.push_back(CRC);

    //Writing calibrations
    DigitizePSCU(fData.size()*2,0x18);
    fraw->WritePSCU(&fDataPSCU);
    fraw->CopyUShortToBuff(&fData);

    //No padding for now...
  }
}

void PamVMCCaloDig::FillCalPedNorm(USBuffer *buff, Int_t sec){

  Int_t chksum = 0;
  for (Int_t plane=0; plane < 11; plane++){
      for (Int_t strip=0; strip < 96; strip++){   
        chksum += (Int_t)fcalped[sec][plane][strip];
        buff->push_back((Int_t)fcalped[sec][plane][strip]);  
        buff->push_back((Int_t)fcalgood[sec][plane][strip]);   
      }
     
  }//3-2114
  
  buff->push_back(((UShort_t)chksum));
  buff->push_back(0x0000);
  buff->push_back(((UShort_t)((Int_t)(chksum >> 16))));
}

void PamVMCCaloDig::FillCalPedBack(USBuffer *buff, Int_t sec){

  Int_t chksum = 0;
  for (Int_t plane=0; plane < 11; plane++){
      for (Int_t strip=95; strip >= 0; strip--){   
        chksum += (Int_t)fcalped[sec][plane][strip];
        buff->push_back((Int_t)fcalped[sec][plane][strip]);  
        buff->push_back((Int_t)fcalgood[sec][plane][strip]);    
      }     
  }//3-2114

  buff->push_back(((UShort_t)chksum));
  buff->push_back(0x0000);
  buff->push_back(((UShort_t)((Int_t)(chksum >> 16))));
}


void PamVMCCaloDig::FillCalThrNorm(USBuffer *buff, Int_t sec){

  Int_t chksum = 0;
  for (Int_t plane=0; plane < 11; plane++){
      for (Int_t strip=0; strip < 6; strip++){   
        chksum += (Int_t)fcalthr[sec][plane][strip];
          buff->push_back(0x0000);  
          buff->push_back((Int_t)fcalthr[sec][plane][strip]);   
      }
     
  }//2118-2249

  buff->push_back(0x0000);
  buff->push_back(((UShort_t)chksum));
  buff->push_back(0x0000);
  buff->push_back(((UShort_t)((Int_t)(chksum >> 16))));
}

void PamVMCCaloDig::FillCalThrBack(USBuffer *buff, Int_t sec){

  Int_t chksum = 0;
  for (Int_t plane=0; plane < 11; plane++){
      for (Int_t strip=5; strip >= 0; strip--){   
        chksum += (Int_t)fcalthr[sec][plane][strip];
          buff->push_back(0x0000);  
          buff->push_back((Int_t)fcalthr[sec][plane][strip]);    
      }     
  }//2118-2249

  buff->push_back(0x0000);
  buff->push_back(((UShort_t)chksum));
  buff->push_back(0x0000);
  buff->push_back(((UShort_t)((Int_t)(chksum >> 16))));
}

void PamVMCCaloDig::FillCalRmsNorm(USBuffer *buff, Int_t sec){

  for (Int_t plane=0; plane < 11; plane++){
      for (Int_t strip=0; strip < 96; strip++){   
          buff->push_back(0x0000);  
          buff->push_back((Int_t)fcalrms[sec][plane][strip]);   
      }
  }//2254-4365
}

void PamVMCCaloDig::FillCalRmsBack(USBuffer *buff, Int_t sec){

  for (Int_t plane=0; plane < 11; plane++){
      for (Int_t strip=95; strip >= 0; strip--){   
          buff->push_back(0x0000);  
          buff->push_back((Int_t)fcalrms[sec][plane][strip]);    
      }     
  }//2254-4365
}

void PamVMCCaloDig::FillCalBaseVarNorm(USBuffer *buff, Int_t sec){

  for (Int_t plane=0; plane < 11; plane++){
      for (Int_t strip=0; strip < 6; strip++){   
          buff->push_back(0x0000);
          buff->push_back(((Int_t)fcalbase[sec][plane][strip]));   
          buff->push_back(0x0000);  
          buff->push_back(((Int_t)fcalvar[sec][plane][strip]));  
      }
  }//4366-4629
}

void PamVMCCaloDig::FillCalBaseVarBack(USBuffer *buff, Int_t sec){

  for (Int_t plane=0; plane < 11; plane++){
      for (Int_t strip=5; strip >= 0; strip--){   
        buff->push_back(0x0000);  
        buff->push_back((Int_t)fcalbase[sec][plane][strip]);   
        buff->push_back(0x0000);  
        buff->push_back((Int_t)fcalvar[sec][plane][strip]);     
      }     
  }//4366-4629
}

void PamVMCCaloDig::Digitize(){

  cout<<"Digitizing CALO..."<<endl;

  Int_t ModCalo = 0; // 0 is RAW, 1 is COMPRESS, 2 is FULL    
                     //####@@@@ should be given as input par @@@@####


  // call different routines depending on the acq mode you want to simulate
  
  switch ( ModCalo ){
  case 0:
    this->DigitizeCaloRaw();
    break;
  case 1:
    this->DigitizeCaloCompress();
    break;
  case 2:
    this->DigitizeCaloFull();
    break;
  };
}

void PamVMCCaloDig::DigitizeCaloRaw(){

  
  // some variables
  
  Float_t ens = 0.;
  UInt_t adcsig = 0;
  UInt_t adcbase = 0;
  UInt_t adc = 0;
  Int_t pre = 0;
  UInt_t l = 0;
  UInt_t lpl = 0;
  Double_t pedenoise;
  Float_t rms = 0.;
  Float_t pedestal = 0.;  
  UInt_t tmpadc[11];
  USBuffer tempdig;

  static const Float_t CALOGeV2MIPratio = 0.0001059994;

  // Header of the four sections
  
  fSecCalo[0] = 0xEA08; // XE
  fSecCalo[1] = 0xF108; // XO
  fSecCalo[2] = 0xF608; // YE 
  fSecCalo[3] = 0xED08; // YO
  
  // length of the data is 0x0428 in RAW mode
  
  fSecCaloLength[0] = 0x0428; // XE
  fSecCaloLength[1] = 0x0428; // XO
  fSecCaloLength[2] = 0x0428; // YE
  fSecCaloLength[3] = 0x0428; // YO

  fData.clear();
  
  for (Int_t sec=0; sec < 4; sec++){
    
    // sec =    0 -> XE      1 -> XO        2-> YE         3 -> YO
    
    tempdig.clear();
    
    // First of all we have section header and packet length
    if (sec !=1){
      tempdig.push_back(fSecCalo[sec]);//1
      tempdig.push_back(fSecCaloLength[sec]);//2
    }
    
    // selftrigger coincidences - in the future we should add here some code 
    // to simulate timing response of pre-amplifiers
    // For now from application point of view we have 1 hit = 1 strip (nHits <=4224)
    // It's possible to introduce 1 hit = 1 particle like in other detectors... if need
    
    for (Int_t autoplane=0; autoplane < 7; autoplane++) tempdig.push_back(0x0000); //3-8
    
    // here comes data
    // Section XO is read in the opposite direction respect to the others
   
    pre = -1;
    
    //l=(sec<2)?1:0;
    l = 0;                 // XE and XO are Y planes 
    if ( sec < 2 ) l = 1;  // while YE and  YO are X planes
    
    for (Int_t strip=0; strip < 96; strip++){   
      
      // which is the pre for this strip?
      
      if (strip%16 == 0) pre++;
      
      for (Int_t plane=0; plane < 11; plane++){
        
        switch (sec){
        case 0:
        case 3:
          lpl = plane * 2;
          break;
        case 1:
        case 2:
          lpl = (plane * 2) + 1;
          break;
        default:
          break;
        }
        
        // get the energy in GeV from the simulation for that strip
        
        ens = GetCaloErel(sec,plane,strip);     
        
        // convert it into ADC channels
        
        adcsig = Int_t(ens*fCalomip[l][lpl][strip]/CALOGeV2MIPratio);
        
        // sum baselines
        
        adcbase = (UInt_t)fcalbase[sec][plane][pre]; 
        
        // add noise and pedestals
        
        pedestal = fcalped[sec][plane][strip];
        rms = fcalrms[sec][plane][strip]/4.;
        
        // Add random gaussian noise of RMS rms and Centered in the pedestal
         
        pedenoise = gRandom->Gaus((Double_t)pedestal,(Double_t)rms); 
        
        // Sum all contribution
        
        adc = adcsig + adcbase + (Int_t)round(pedenoise);
        
        // Signal saturation
        
        if ( adc > 0x7FFF ) adc = 0x7FFF;
        
        // save value
        if (adcsig>0){
          cout <<"STRIP HIT"<<endl
               <<"SEC: "<<sec<<" PLANE: "<<plane<<" STRIP:"<<strip<<endl 
               <<"adcsig: "<< adcsig<<" adc: "<<adc<<endl;
        }
        if ( sec != 1 ) tempdig.push_back(adc);
        else tmpadc[plane]= adc;
        
      }

      if ( sec == 1 ) for (Int_t i=10; i>=0; i--) tempdig.push_back(tmpadc[i]);
    }
    
    if ( sec == 1 ){
      tempdig.push_back(fSecCaloLength[1]);//2
      tempdig.push_back(fSecCalo[1]);//1 
      reverse(tempdig.begin(), tempdig.end());
    }
    // here we calculate and save the CRC
    //++++++ COPY TO DIGITIZER'S BUFFER +++//
    Short_t CRC = 0;

    USBuffer::const_iterator p = tempdig.begin();
    while( p!= tempdig.end() ){
      CRC=crc(CRC,*p);
      fData.push_back(*p);
      p++;
    }
    fData.push_back(CRC);
  }
}


Float_t PamVMCCaloDig::GetCaloErel(Int_t sec, Int_t plane, Int_t strip){
  
  // determine plane and strip
  
  Int_t mplane = 0;

  switch (sec){
  case 0:
    mplane = plane * 4 + 1; // it must be 0, 4, 8, ... (+1)  from plane = 0, 11
    break;
  case 1:
    mplane = plane * 4 + 2 + 1; // it must be 2, 6, 10, ... (+1)  from plane = 0, 11
    break;
  case 2:
    mplane = plane * 4 + 3 + 1; // it must be 3, 7, 11, ... (+1)  from plane = 0, 11
    break;
  case 3:
    mplane = plane * 4 + 1 + 1; // it must be 1, 5, 9, ... (+1)  from plane = 0, 11
    break;
  }
  
 
  Float_t Erel = 0.;


  // search energy release in gpamela output
  if (fhc){ 
    PamVMCDetectorHit *hit = (PamVMCDetectorHit*)fhc->At((mplane-1)*96+strip);
    if (hit) Erel=hit->GetEREL();
  } else cout<<"CALO HIT Collection pointer not found"<<endl;

  //if (sec == 3) Erel = 0.01; else Erel = 0.;

  return Erel;
}


void PamVMCCaloDig::DigitizeCaloCompress(){

  cout<<"!!!WARNING PamVMCCaloDig COMPRESS MODE STILL NOT IMPLEMENTED!!!"<<endl;  
  
  this->DigitizeCaloRaw();
  return;

}

void PamVMCCaloDig::DigitizeCaloFull(){

  cout<<"!!!WARNING PamVMCCaloDig FULL MODE STILL NOT IMPLEMENTED!!!"<<endl;  
  
  this->DigitizeCaloRaw();
  return;

}
1	pamelats	1.1	#include "PamVMCCaloDig.h"
2			#include <TTree.h>
3			#include "CalibCalPedEvent.h"
4
5			ClassImp(PamVMCCaloDig)
6
7			extern "C"{
8			short crc(short, short);
9			};
10
11
12			void PamVMCCaloDig::LoadCalib(){
13
14			cout<<"Loading Tracker Calibrations..."<<endl;
15
16
17
18			fhc =(TClonesArray*)fhitscolmap.GetValue("CAST");
19			if (!fhc) cout <<"!!!WARNING Hit Collection of CAL not found by PamVMCCaloDig!!!!"<<endl;
20
21			Int_t GivenCaloCalib = 0; //@@@ should be given as input par @@@
22
23			Int_t time=3;
24			Int_t type=101;
25
26			fdberr = fsql->Query_GL_PARAM(time,type);
27
28			if(fdberr<0){
29			cout<<"No such record in DB for CAL: time="<<time<<" type="<<type<<endl;
30			ThrowCalFileWarning("CAL");
31
32			} else {
33
34			fquery.str("");
35			fquery << fsql->GetPAR()->PATH.Data() << "/";
36			fquery << fsql->GetPAR()->NAME.Data();
37
38			ThrowCalFileUsage("CAL",fquery.str().c_str());
39
40			FILE *f;
41			f = fopen(fquery.str().c_str(),"rb");
42
43			if(f==NULL) ThrowCalFileWarning("CAL");
44			else{
45			memset(fCalomip,0,4224*sizeof(fCalomip[0][0][0]));
46
47			for (Int_t m = 0; m < 2 ; m++ ){
48			for (Int_t k = 0; k < 22; k++ ){
49			for (Int_t l = 0; l < 96; l++ ){
50			fread(&fCalomip[m][k][l],sizeof(fCalomip[m][k][l]),1,f);
51			};
52			};
53			};
54			fclose(f);
55
56			// determine which calibration has to be used
57			// and load it for each section
58
59
60			UInt_t idcalib, calibno;
61			UInt_t utime = 0;
62			TString calname;
63			for (UInt_t s=0; s<4; s++){
64
65			// clear calo calib variables for section s
66
67			ClearCaloCalib(s);
68
69			if ( GivenCaloCalib ){
70
71			// a time has been given as input on the command line
72			// so retrieve the calibration that preceed that time
73
74			fsql->Query_GL_CALO_CALIB(GivenCaloCalib,utime,s);
75
76			calibno = fsql->GetCaloCalib()->EV_ROOT;
77			idcalib = fsql->GetCaloCalib()->ID_ROOT_L0;
78
79			// determine path and name and entry of the calibration file
80
81			printf("\n");
82			printf(" SECTION %i \n",s);
83
84			fsql->Query_GL_ROOT(idcalib);
85
86			fquery.str("");
87			fquery << fsql->GetROOT()->PATH.Data() << "/";
88			fquery << fsql->GetROOT()->NAME.Data();
89
90			calname = (TString)fquery.str().c_str();
91
92			printf("\n Section %i : using file %s calibration at entry %i: \n",s,calname.Data(),calibno);
93
94			} else {
95
96			fdberr = 0;
97			fdberr = fsql->Query_GL_PARAM(1,104);
98
99			fquery.str("");
100			fquery << fsql->GetPAR()->PATH.Data() << "/";
101			fquery << fsql->GetPAR()->NAME.Data();
102
103			printf("\n Section %i : using default calorimeter calibration: \n %s \n",s,fquery.str().c_str());
104
105			calname = (TString)fquery.str().c_str();
106			calibno = s;
107
108			};
109
110			// load calibration variables in memory
111
112			CaloLoadCalib(s,calname,calibno);
113
114			};
115			//
116			// at this point we have in memory the calorimeter calibration
117			// and we can save it to disk in the correct format and use it to digitize the data
118			//
119
120			DigitizeCALOCALIB();
121			}
122			}
123			}
124
125
126			void PamVMCCaloDig::ClearCaloCalib(Int_t s){
127
128			fcstwerr[s] = 0;
129			fcperror[s] = 0.;
130			for ( Int_t d=0 ; d<11 ;d++ ){
131			Int_t pre = -1;
132			for ( Int_t j=0; j<96 ;j++){
133			if ( j%16 == 0 ) pre++;
134			fcalped[s][d][j] = 0.;
135			fcstwerr[s] = 0.;
136			fcperror[s] = 0.;
137			fcalgood[s][d][j] = 0.;
138			fcalthr[s][d][pre] = 0.;
139			fcalrms[s][d][j] = 0.;
140			fcalbase[s][d][pre] = 0.;
141			fcalvar[s][d][pre] = 0.;
142			};
143			};
144			return;
145			}
146
147			Int_t PamVMCCaloDig::CaloLoadCalib(Int_t s,TString calname, UInt_t calibno){
148
149			UInt_t e = 0;
150			switch(s){
151			case 0: e = 0;
152			break;
153			case 1: e = 2;
154			break;
155			case 2: e = 3;
156			break;
157			case 3: e = 1;
158			break;
159			default:
160			break;
161			}
162			fcfile.open(calname.Data());
163			ThrowCalFileUsage("CAL",calname.Data());
164			if ( !fcfile ){
165			ThrowCalFileWarning("CAL");
166			return(-107);
167			};
168			fcfile.close();
169
170			fcrfile = new TFile(calname.Data());
171			if ( !fcrfile ){
172			ThrowCalFileWarning("CAL");
173			return(-108);
174			}
175			TTree tr = (TTree)fcrfile->Get("CalibCalPed");
176			if ( !tr ){
177			cout<<"!!!WARNING Tree CalibCalPed in file "<<calname.Data()
178			<<" was NOT found!!!"<<endl;
179			return(-109);
180			}
181
182			TBranch *calo = tr->GetBranch("CalibCalPed");
183
184			pamela::CalibCalPedEvent *ce = 0;
185			tr->SetBranchAddress("CalibCalPed", &ce);
186
187			Long64_t ncalibs = calo->GetEntries();
188
189			if ( !ncalibs ){
190			cout<<"!!!WARNING No entries in calibration file!!!"<<endl;
191			return(-110);
192			}
193
194			calo->GetEntry(calibno);
195
196			if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
197			fcstwerr[s] = ce->cstwerr[s];
198			fcperror[s] = ce->cperror[s];
199			for ( Int_t d=0 ; d<11 ;d++ ){
200			Int_t pre = -1;
201			for ( Int_t j=0; j<96 ;j++){
202			if ( j%16 == 0 ) pre++;
203			fcalped[s][d][j] = ce->calped[e][d][j];
204			fcalgood[s][d][j] = ce->calgood[e][d][j];
205			fcalthr[s][d][pre] = ce->calthr[e][d][pre];
206			fcalrms[s][d][j] = ce->calrms[e][d][j];
207			fcalbase[s][d][pre] = ce->calbase[e][d][pre];
208			fcalvar[s][d][pre] = ce->calvar[e][d][pre];
209			};
210			};
211			} else {
212			printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
213			fcrfile->Close();
214			return(-111);
215			};
216			fcrfile->Close();
217			return(0);
218			}
219
220			void PamVMCCaloDig::DigitizeCALOCALIB() {
221
222
223			// Header of the four sections
224
225			fSecCalo[0] = 0xAA00; // XE
226			fSecCalo[1] = 0xB100; // XO
227			fSecCalo[2] = 0xB600; // YE
228			fSecCalo[3] = 0xAD00; // YO
229
230			// length of the data is 0x1215
231
232			fSecCaloLength[0] = 0x1215; // XE
233			fSecCaloLength[1] = 0x1215; // XO
234			fSecCaloLength[2] = 0x1215; // YE
235			fSecCaloLength[3] = 0x1215; // YO
236
237			Short_t CRC;
238
239			for(Int_t sec=0; sec < 4; sec++){
240			//
241			// sec = 0 -> XE 1 -> XO 2-> YE 3 -> YO
242			//
243			fData.clear();
244
245			fData.push_back(fSecCalo[sec]);
246			fData.push_back(fSecCaloLength[sec]);
247
248			if ( sec==1 ){
249			FillCalPedBack(&fData,sec);
250			FillCalThrBack(&fData,sec);
251			FillCalRmsBack(&fData,sec);
252			FillCalBaseVarBack(&fData,sec);
253			} else {
254			FillCalPedNorm(&fData,sec);
255			FillCalThrNorm(&fData,sec);
256			FillCalRmsNorm(&fData,sec);
257			FillCalBaseVarNorm(&fData,sec);
258			}
259			//CRC
260
261			fData.push_back(0x0000);
262			CRC = 0;
263
264
265			USBuffer::const_iterator p = fData.begin();
266			while( p!= fData.end() ){
267			CRC=crc(CRC,*p);
268			p++;
269			}
270			fData.push_back(CRC);
271
272			//Writing calibrations
273			DigitizePSCU(fData.size()*2,0x18);
274			fraw->WritePSCU(&fDataPSCU);
275			fraw->CopyUShortToBuff(&fData);
276
277			//No padding for now...
278			}
279			}
280
281			void PamVMCCaloDig::FillCalPedNorm(USBuffer *buff, Int_t sec){
282
283			Int_t chksum = 0;
284			for (Int_t plane=0; plane < 11; plane++){
285			for (Int_t strip=0; strip < 96; strip++){
286			chksum += (Int_t)fcalped[sec][plane][strip];
287			buff->push_back((Int_t)fcalped[sec][plane][strip]);
288			buff->push_back((Int_t)fcalgood[sec][plane][strip]);
289			}
290
291			}//3-2114
292
293			buff->push_back(((UShort_t)chksum));
294			buff->push_back(0x0000);
295			buff->push_back(((UShort_t)((Int_t)(chksum >> 16))));
296			}
297
298			void PamVMCCaloDig::FillCalPedBack(USBuffer *buff, Int_t sec){
299
300			Int_t chksum = 0;
301			for (Int_t plane=0; plane < 11; plane++){
302			for (Int_t strip=95; strip >= 0; strip--){
303			chksum += (Int_t)fcalped[sec][plane][strip];
304			buff->push_back((Int_t)fcalped[sec][plane][strip]);
305			buff->push_back((Int_t)fcalgood[sec][plane][strip]);
306			}
307			}//3-2114
308
309			buff->push_back(((UShort_t)chksum));
310			buff->push_back(0x0000);
311			buff->push_back(((UShort_t)((Int_t)(chksum >> 16))));
312			}
313
314
315			void PamVMCCaloDig::FillCalThrNorm(USBuffer *buff, Int_t sec){
316
317			Int_t chksum = 0;
318			for (Int_t plane=0; plane < 11; plane++){
319			for (Int_t strip=0; strip < 6; strip++){
320			chksum += (Int_t)fcalthr[sec][plane][strip];
321			buff->push_back(0x0000);
322			buff->push_back((Int_t)fcalthr[sec][plane][strip]);
323			}
324
325			}//2118-2249
326
327			buff->push_back(0x0000);
328			buff->push_back(((UShort_t)chksum));
329			buff->push_back(0x0000);
330			buff->push_back(((UShort_t)((Int_t)(chksum >> 16))));
331			}
332
333			void PamVMCCaloDig::FillCalThrBack(USBuffer *buff, Int_t sec){
334
335			Int_t chksum = 0;
336			for (Int_t plane=0; plane < 11; plane++){
337			for (Int_t strip=5; strip >= 0; strip--){
338			chksum += (Int_t)fcalthr[sec][plane][strip];
339			buff->push_back(0x0000);
340			buff->push_back((Int_t)fcalthr[sec][plane][strip]);
341			}
342			}//2118-2249
343
344			buff->push_back(0x0000);
345			buff->push_back(((UShort_t)chksum));
346			buff->push_back(0x0000);
347			buff->push_back(((UShort_t)((Int_t)(chksum >> 16))));
348			}
349
350			void PamVMCCaloDig::FillCalRmsNorm(USBuffer *buff, Int_t sec){
351
352			for (Int_t plane=0; plane < 11; plane++){
353			for (Int_t strip=0; strip < 96; strip++){
354			buff->push_back(0x0000);
355			buff->push_back((Int_t)fcalrms[sec][plane][strip]);
356			}
357			}//2254-4365
358			}
359
360			void PamVMCCaloDig::FillCalRmsBack(USBuffer *buff, Int_t sec){
361
362			for (Int_t plane=0; plane < 11; plane++){
363			for (Int_t strip=95; strip >= 0; strip--){
364			buff->push_back(0x0000);
365			buff->push_back((Int_t)fcalrms[sec][plane][strip]);
366			}
367			}//2254-4365
368			}
369
370			void PamVMCCaloDig::FillCalBaseVarNorm(USBuffer *buff, Int_t sec){
371
372			for (Int_t plane=0; plane < 11; plane++){
373			for (Int_t strip=0; strip < 6; strip++){
374			buff->push_back(0x0000);
375			buff->push_back(((Int_t)fcalbase[sec][plane][strip]));
376			buff->push_back(0x0000);
377			buff->push_back(((Int_t)fcalvar[sec][plane][strip]));
378			}
379			}//4366-4629
380			}
381
382			void PamVMCCaloDig::FillCalBaseVarBack(USBuffer *buff, Int_t sec){
383
384			for (Int_t plane=0; plane < 11; plane++){
385			for (Int_t strip=5; strip >= 0; strip--){
386			buff->push_back(0x0000);
387			buff->push_back((Int_t)fcalbase[sec][plane][strip]);
388			buff->push_back(0x0000);
389			buff->push_back((Int_t)fcalvar[sec][plane][strip]);
390			}
391			}//4366-4629
392			}
393
394			void PamVMCCaloDig::Digitize(){
395
396			cout<<"Digitizing CALO..."<<endl;
397
398			Int_t ModCalo = 0; // 0 is RAW, 1 is COMPRESS, 2 is FULL
399			//####@@@@ should be given as input par @@@@####
400
401
402			// call different routines depending on the acq mode you want to simulate
403
404			switch ( ModCalo ){
405			case 0:
406			this->DigitizeCaloRaw();
407			break;
408			case 1:
409			this->DigitizeCaloCompress();
410			break;
411			case 2:
412			this->DigitizeCaloFull();
413			break;
414			};
415			}
416
417			void PamVMCCaloDig::DigitizeCaloRaw(){
418
419
420			// some variables
421
422			Float_t ens = 0.;
423			UInt_t adcsig = 0;
424			UInt_t adcbase = 0;
425			UInt_t adc = 0;
426			Int_t pre = 0;
427			UInt_t l = 0;
428			UInt_t lpl = 0;
429			Double_t pedenoise;
430			Float_t rms = 0.;
431			Float_t pedestal = 0.;
432			UInt_t tmpadc[11];
433			USBuffer tempdig;
434
435			static const Float_t CALOGeV2MIPratio = 0.0001059994;
436
437			// Header of the four sections
438
439			fSecCalo[0] = 0xEA08; // XE
440			fSecCalo[1] = 0xF108; // XO
441			fSecCalo[2] = 0xF608; // YE
442			fSecCalo[3] = 0xED08; // YO
443
444			// length of the data is 0x0428 in RAW mode
445
446			fSecCaloLength[0] = 0x0428; // XE
447			fSecCaloLength[1] = 0x0428; // XO
448			fSecCaloLength[2] = 0x0428; // YE
449			fSecCaloLength[3] = 0x0428; // YO
450
451			fData.clear();
452
453			for (Int_t sec=0; sec < 4; sec++){
454
455			// sec = 0 -> XE 1 -> XO 2-> YE 3 -> YO
456
457			tempdig.clear();
458
459			// First of all we have section header and packet length
460			if (sec !=1){
461			tempdig.push_back(fSecCalo[sec]);//1
462			tempdig.push_back(fSecCaloLength[sec]);//2
463			}
464
465			// selftrigger coincidences - in the future we should add here some code
466			// to simulate timing response of pre-amplifiers
467			// For now from application point of view we have 1 hit = 1 strip (nHits <=4224)
468			// It's possible to introduce 1 hit = 1 particle like in other detectors... if need
469
470			for (Int_t autoplane=0; autoplane < 7; autoplane++) tempdig.push_back(0x0000); //3-8
471
472			// here comes data
473			// Section XO is read in the opposite direction respect to the others
474
475			pre = -1;
476
477			//l=(sec<2)?1:0;
478			l = 0; // XE and XO are Y planes
479			if ( sec < 2 ) l = 1; // while YE and YO are X planes
480
481			for (Int_t strip=0; strip < 96; strip++){
482
483			// which is the pre for this strip?
484
485			if (strip%16 == 0) pre++;
486
487			for (Int_t plane=0; plane < 11; plane++){
488
489			switch (sec){
490			case 0:
491			case 3:
492			lpl = plane * 2;
493			break;
494			case 1:
495			case 2:
496			lpl = (plane * 2) + 1;
497			break;
498			default:
499			break;
500			}
501
502			// get the energy in GeV from the simulation for that strip
503
504			ens = GetCaloErel(sec,plane,strip);
505
506			// convert it into ADC channels
507
508			adcsig = Int_t(ens*fCalomip[l][lpl][strip]/CALOGeV2MIPratio);
509
510			// sum baselines
511
512			adcbase = (UInt_t)fcalbase[sec][plane][pre];
513
514			// add noise and pedestals
515
516			pedestal = fcalped[sec][plane][strip];
517			rms = fcalrms[sec][plane][strip]/4.;
518
519			// Add random gaussian noise of RMS rms and Centered in the pedestal
520
521			pedenoise = gRandom->Gaus((Double_t)pedestal,(Double_t)rms);
522
523			// Sum all contribution
524
525			adc = adcsig + adcbase + (Int_t)round(pedenoise);
526
527			// Signal saturation
528
529			if ( adc > 0x7FFF ) adc = 0x7FFF;
530
531			// save value
532			if (adcsig>0){
533			cout <<"STRIP HIT"<<endl
534			<<"SEC: "<<sec<<" PLANE: "<<plane<<" STRIP:"<<strip<<endl
535			<<"adcsig: "<< adcsig<<" adc: "<<adc<<endl;
536			}
537			if ( sec != 1 ) tempdig.push_back(adc);
538			else tmpadc[plane]= adc;
539
540			}
541
542			if ( sec == 1 ) for (Int_t i=10; i>=0; i--) tempdig.push_back(tmpadc[i]);
543			}
544
545			if ( sec == 1 ){
546			tempdig.push_back(fSecCaloLength[1]);//2
547			tempdig.push_back(fSecCalo[1]);//1
548			reverse(tempdig.begin(), tempdig.end());
549			}
550			// here we calculate and save the CRC
551			//++++++ COPY TO DIGITIZER'S BUFFER +++//
552			Short_t CRC = 0;
553
554			USBuffer::const_iterator p = tempdig.begin();
555			while( p!= tempdig.end() ){
556			CRC=crc(CRC,*p);
557			fData.push_back(*p);
558			p++;
559			}
560			fData.push_back(CRC);
561			}
562			}
563
564
565			Float_t PamVMCCaloDig::GetCaloErel(Int_t sec, Int_t plane, Int_t strip){
566
567			// determine plane and strip
568
569			Int_t mplane = 0;
570
571			switch (sec){
572			case 0:
573			mplane = plane * 4 + 1; // it must be 0, 4, 8, ... (+1) from plane = 0, 11
574			break;
575			case 1:
576			mplane = plane * 4 + 2 + 1; // it must be 2, 6, 10, ... (+1) from plane = 0, 11
577			break;
578			case 2:
579			mplane = plane * 4 + 3 + 1; // it must be 3, 7, 11, ... (+1) from plane = 0, 11
580			break;
581			case 3:
582			mplane = plane * 4 + 1 + 1; // it must be 1, 5, 9, ... (+1) from plane = 0, 11
583			break;
584			}
585
586
587			Float_t Erel = 0.;
588
589
590			// search energy release in gpamela output
591			if (fhc){
592			PamVMCDetectorHit hit = (PamVMCDetectorHit)fhc->At((mplane-1)*96+strip);
593			if (hit) Erel=hit->GetEREL();
594			} else cout<<"CALO HIT Collection pointer not found"<<endl;
595
596			//if (sec == 3) Erel = 0.01; else Erel = 0.;
597
598			return Erel;
599			}
600
601
602
603			void PamVMCCaloDig::DigitizeCaloCompress(){
604
605			cout<<"!!!WARNING PamVMCCaloDig COMPRESS MODE STILL NOT IMPLEMENTED!!!"<<endl;
606
607			this->DigitizeCaloRaw();
608			return;
609
610			}
611
612			void PamVMCCaloDig::DigitizeCaloFull(){
613
614			cout<<"!!!WARNING PamVMCCaloDig FULL MODE STILL NOT IMPLEMENTED!!!"<<endl;
615
616			this->DigitizeCaloRaw();
617			return;
618
619			}