source: trunk/MagicSoft/Simulation/Detector/Starfield/convertcorsika.cxx@ 10100

//////////////////////////////////////////////////////////////////////////////////////
//     Write photon data in binary corsika-like format. 
//     Output will be used as an input for the reflector as if it was a shower event.
//
//////////////////////////////////////////////////////////////////////////////////////

#include "convertcorsika.h"

COREventHeader cerevth;
CORParticle cerphot;
CORStatfile cerstat;
float outbuf[273];

int convertcorsika(int id, int photnum, photon phot[], float inttime_s, int verbose,char output_name[]){

  int i,filenum;
  int lastblock;
  float x_cm, y_cm, u, v, lambda_nm, t_ns;
  char cor_dir[60];
  char stat_dir[60];  
  char cor_file[20]="cer";
  char stat_file[20]="sta";
  
  filenum=id; 

  // Event corsika file
  
  //File labeling.
  
  strcpy (cor_dir, "./");
  
  strcat(cor_file, output_name);
    
  strcat(cor_dir, cor_file);

  
  // Fill the header of the corsika-like event. 
  //The fields in evt.fill are: event number, primary identifier, total energy (GeV),
  // first target identifier, first z interaction(cm) , momentum x, momentum y, momentum z, 
  // zenith angle, azimuth angle, coreposition x (cm), coreposition y(cm)
  
  ofstream cerfile;
  //FILE cerfilec;
  
  if(!cerfile) {
     cerr<<"Cannot create cerenkov file.\n";
    exit(1);
    } 
  
  cerfile.open(cor_dir, ios::out|ios::binary);

  /*if((cerfilec=fopen(cor_dir,"wb"))==NULL){
    printf("C-style::Cannot create cerenkov file.\n");
    exit(1);
    }*/
  
  //____________________________________________________________________________________
  
  cout<<"Writing binary Cherenkov file"<<" "<<cor_dir<<" "<<"..."<<endl;
  
  cerevth.fill (1.0, 
                99.0, // primary id
                inttime_s * 1e9, // instead of the total energy: integration time (ns)
                1.0,
                1.e7,
                0.0,
                0.0,
                0.0,
                0.0, // theta is always 0. !
                0.0, // phi is always 0. !
                0.0,
                0.0);
  strcpy(cerevth.EVTH, "RUNH");
                
  cerevth.write(cerfile);

  strcpy(cerevth.EVTH, "EVTH");
  cerevth.write(cerfile);
                
  // Here we fill the information for each photon as a corsika particle. 
  // The fields in cerphot.fill are: wavelentgh (nm), x core position (cm), y core position (cm),
  // u cosine director, v cosine director, time since first interaction (ns), 
  // height (cm) of first interaction.          
  //photnum=0;
  for(i=0;i<photnum;i++){  
    
    lambda_nm = phot[i].lambda_nm;    
    x_cm = phot[i].x_m * 100.0;
    y_cm = phot[i].y_m * 100.0;
 
    u = phot[i].u;
    v = phot[i].v;  
    t_ns = phot[i].arrtime_sec * 1e9;

    cerphot.fill(lambda_nm,
                 x_cm,
                 y_cm,
                 u,
                 v,
                 t_ns,
                 1.e7); 
    
     
    cerphot.write(cerfile);
  }
  lastblock=39-photnum%39;

  for(i=0;i<lastblock;i++){
  cerphot.fill(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
  cerphot.write(cerfile);  
  }
  strcpy(cerevth.EVTH, "EVTE");
  cerevth.write(cerfile);
  strcpy(cerevth.EVTH, "RUNE");
  cerevth.write(cerfile);

  cerfile.close();  
  
  cout<<"Done."<<endl;
  
  // Statistics corsika file
  
  strcpy (stat_dir, "./");
    
  //File labeling. 
  
  strcat(stat_file, output_name);
    
  strcat(stat_dir, stat_file);
  
  cerfile.open (stat_dir, ios::out|ios::binary);
  
  cout<<"Writing binary statistics file "<<" "<<stat_dir<<" "<<"..."<<endl;
  
  cerstat.write(cerfile); 
  
  cerfile.close(); 
  
  cout<<"Done."<<endl;  
  
  return 0;
      
}