Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

CorsikaFileShowerGenerator - source file

// $Id: CorsikaFileShowerGenerator.cc,v 1.10 2005/10/21 13:31:30 moreggia Exp $
// Author: M.Dattoli   2004/08/31 

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  Id: CorsikaFileShowerGenerator                                           *
 *  Package: ShowerGenerators                                                *
 *  Coordinator: Sergio.Bottai, Dmitry.Naumov                                *
 *                                                                           *
 *****************************************************************************/

//_____________________________________________________________________________
//   Corsika File Shower Generator
//   =============================
//
//   <extensive description>
//
//   Config file parameters
//   ======================
//   
//   FlagThinning [string]:  flag for thinning option; 
//			     can be t (thinning enabled)
//			     or f (thinning disabled)
//
//   HowManyShowers [int] : number of showers simulated in corsika run
//
//   ParticleFileName [string]: name of corsika normal particle output file
//
//   LongFileName [string]: name of corsika longitudinal output file
//

#include "CorsikaFileShowerGenerator.hh"
#include "MCTruth.hh"
#include <ctime>
#include <cstdlib>
#include <cmath>
 
using namespace sou;
using namespace TMath;
using namespace EConst;

ClassImp(CorsikaFileShowerGenerator)

//______________________________________________________________________________
 CorsikaFileShowerGenerator::CorsikaFileShowerGenerator(): 
           EventGenerator("Corsikafile"),EsafMsgSource(), fTrack(0), 
           fParticleFile(NULL), fLongFile(NULL), fTruth(0) 
    // 
    // Constructor
    // 
{
    fParticleFileName = "";
    fLongFileName = "";
}


//______________________________________________________________________________
 CorsikaFileShowerGenerator::~CorsikaFileShowerGenerator() 
    // 
    // Destructor
    // 
{
    if ( fTrack )  delete fTrack;
    if ( fTruth )  delete fTruth;
        fTrack=NULL;
        fTruth=NULL;
}



//______________________________________________________________________________
 void CorsikaFileShowerGenerator::ReplaceInputFile(const char* pfn,const char* lfn) 
    //
    // Replace input files
    // 
{

    // Particle file
    if ( fParticleFile )
	 Msg(EsafMsg::Panic)<<"Can not change name when particle file is already open!"<<MsgDispatch;
    if (!pfn) 
    {
	Msg(EsafMsg::Panic)<<" Particle file name ignored!" <<MsgDispatch;
        return;
    }
    Conf()->ReplaceStr("CorsikaFileShowerGenerator.ParticleFileName",pfn);
    Msg(EsafMsg::Info)<<"Input file name changed to" << pfn <<MsgDispatch;

    //Longitudinal file: info about energy loss, n electrons, ncharged, ncherenkov, depth
    if ( fLongFile )
	Msg(EsafMsg::Panic)<<"Can not change name when particle file is already open!"<<MsgDispatch;
    if (!lfn) 
    { 
	Msg(EsafMsg::Panic)<<" Particle file name ignored!" <<MsgDispatch;
        return;
    }
    Conf()->ReplaceStr("CorsikaFileShowerGenerator.LongFileName",lfn);
    Msg(EsafMsg::Info)<<"Input file name changed to" << lfn <<MsgDispatch;

}



//______________________________________________________________________________
 void CorsikaFileShowerGenerator::Open() 
    // 
    // Open file
    // 
{

  //particle file
    if ( Conf()->IsDefined( "CorsikaFileShowerGenerator.ParticleFileName" ) )
        fParticleFileName = Conf()->GetStr("CorsikaFileShowerGenerator.ParticleFileName");

    else
	Msg(EsafMsg::Panic)<<"Unknown input particle file name in CorsikaFileShowerGenerator"<<MsgDispatch;

    fParticleFile = fopen( fParticleFileName.c_str(), "rb" );

    if( fParticleFile == NULL ) 
	Msg(EsafMsg::Panic)<<"Can't open file "+fParticleFileName << MsgDispatch;
    else
        Msg(EsafMsg::Info) << "File: " << fParticleFileName << " successfully opened " << MsgDispatch;

  //longitudinal file
    if ( Conf()->IsDefined( "CorsikaFileShowerGenerator.LongFileName" ) )
        fLongFileName = Conf()->GetStr("CorsikaFileShowerGenerator.LongFileName");

    else
	Msg(EsafMsg::Panic)<<"Unknown input longitudinal ditribution file name in CorsikaFileShowerGenerator"<<MsgDispatch;

    fLongFile = fopen( fLongFileName.c_str(), "rb" );

    if( fLongFile == NULL ) 
	Msg(EsafMsg::Panic)<<"Can't open file "+fLongFileName << MsgDispatch;
    else
        Msg(EsafMsg::Info) << "File: " << fLongFileName << " successfully opened " << MsgDispatch;


}


//______________________________________________________________________________
 void CorsikaFileShowerGenerator::Close() 
    // 
    // Close file
    // 
{

    if( fParticleFile ) fclose(fParticleFile ) ;
    if( fLongFile  ) fclose(fLongFile ) ;

    
}


//______________________________________________________________________________
 PhysicsData* CorsikaFileShowerGenerator::Get() 
{
    // 
    // Get a new event
    // 

    // build and return a track
    if ( !fParticleFile || !fLongFile) Open();
 
    if ( Conf()->IsDefined( "CorsikaFileShowerGenerator.FlagThinning" ) )
        thinning = Conf()->GetStr("CorsikaFileShowerGenerator.FlagThinning");
    else Msg(EsafMsg::Panic) <<"thinning flag missing in config file" << MsgDispatch;
    if (thinning=="f") {
	n1 = 272;
	nmax = 7;
	Msg(EsafMsg::Info)<<"Thinning disabled"<<MsgDispatch;
	}

	else if (thinning=="t") {
	n1 = 311;
	nmax = 8;
	Msg(EsafMsg::Info)<<"Thinning enabled"<<MsgDispatch;
	}

	n2 = n1 + 1;
 
        ConfigFileParser *pConfig = Config::Get()->GetCF("Generators","CorsikaFileShowerGenerator");
        howmany = pConfig->GetNum("CorsikaFileShowerGenerator.HowManyShowers");

        //random choice of event to be processed among simulated events
        srand(time(0));
        event =(rand()%(Int_t)(howmany))+1;
        Msg(EsafMsg::Debug) << "Processing event n : " << event << MsgDispatch;
 

       if ( LoadTrack(event) ) {
        Msg(EsafMsg::Debug) << "Load Track Event" << MsgDispatch;    
        return fTrack;
        }
       else {
        Msg(EsafMsg::Panic) << "End of file " << fParticleFileName << " reached? " << MsgDispatch;
        return (ShowerTrack*)0;
     }
 
       return fTrack;   
       }

//______________________________________________________________________________
 void CorsikaFileShowerGenerator::ReadParticleHeader() 
{
    //
    // Read particle file header.
    //

  Msg(EsafMsg::Debug) << "Reading particle header.." << MsgDispatch;
    while(!feof(fParticleFile))
    {
        //Reading BlockHeader (273,21) or (312,21)

        nin = fread(&r,4,1,fParticleFile);
        for (k=1; k<=21;k++)
        {
            pos = ftell(fParticleFile); 

            //Reading SubBlockHeader   
            nin = fread(&header,sizeof(header),1,fParticleFile);

            //Reading data in SubBlock; Storing them into map
            if (header[0]=='R'&&header[3]=='H')  
            {

                float data_runh[n1];
                for (i=0; i<n1; i++) {data_runh[i]=0;}
                //Reading RUNH Block (once for run)
                nin = fread (&data_runh, sizeof(data_runh), 1, fParticleFile);
		for (i=0; i<=4; i++) {
                aatm[i]=data_runh[253+i];
                batm[i]=data_runh[258+i];
                catm[i]=data_runh[263+i];
		}
            }  // end of if "runh"

            else if (header[0]=='R'&&header[3]=='E') 
            {
                float data_rune[n1];

                for (i=0; i<n1; i++) {data_rune[i]=0;}
                //reading RUNE Block (once per run)
                nin = fread (&data_rune, sizeof(data_rune), 1,fParticleFile );
            } // end of if "rune"

            else if (header[0]=='E'&&header[3]=='H') 
            {
                float data_evth[n1];

                for (i=0; i<n1; i++) {data_evth[i]=0;}

                //Reading EVTH Block (1 for event)
                nin = fread (&data_evth, sizeof(data_evth), 1,fParticleFile );
                ev_num= (Int_t)(data_evth[0]);

                if (ev_num==event) 
                {
		Msg(EsafMsg::Debug) << "Reading event n " << ev_num << MsgDispatch;
                    z_firstint = data_evth[5];  // cm
		    Msg(EsafMsg::Debug)<< "z_firstint: " << z_firstint <<MsgDispatch; 
                    thick_firstint = GetThickness(abs(z_firstint)); //g/cm^2
		    Msg(EsafMsg::Debug)<< "thick_firstint: " << thick_firstint <<MsgDispatch; 
                    epolangd = data_evth[9]*180/Pi();
                    eaziangd = data_evth[10]*180/Pi();
                    eldircos = Sin(data_evth[9])*Cos(data_evth[10]);
                    emdircos = Sin(data_evth[9])*Sin(data_evth[10]);
                    endircos = Cos(data_evth[9]);
		    energy = data_evth[2];
                    Msg (EsafMsg::Debug) << "epolangd: " << epolangd
					<< "n eaziangd: " << eaziangd 
					<< "n eldircos: " << eldircos 
					<< "n emdircos: " << emdircos 
					<< "n endircos: " << endircos 
					<< "n energy: " << energy << MsgDispatch;
                    
                }  //end of if (ev_num==event)
            } //end of if "evth"

            else if (header[0]=='E'&&header[3]=='E') 
            {
                float data_evte[n1];

                for (i=0; i<n1; i++) {data_evte[i]=0;}
                //Reading EVTE Block (1 for event)
                nin = fread (&data_evte, sizeof(data_evte), 1, fParticleFile);

            } // end of if "evte"

            else if (header[0]=='L'&&header[3]=='G') 
            {
                float data_long[n1];
                for (i=0; i<n1; i++) {data_long[i]=0;}

                //Reading LONG Block
                nin = fread (&data_long, sizeof(data_long), 1, fParticleFile);
            } // end of if "long"

            else 
            {
                fseek(fParticleFile,pos,SEEK_SET);
                float data_block[n2];
                for (i=0; i<n2; i++) {data_block[i]=0;}

                //Reading DATA BLOCK
                nin = fread (&data_block, sizeof(data_block), 1,fParticleFile );

            } // end of else... DATA BLOCK

        } // end of for (k=1; k<=21;k++)

        k=1;

        //Reading BlockEnd  (273,21) o (312,21) 
        nin = fread(&r,4,1,fParticleFile);
    }

}

//______________________________________________________________________________
 void CorsikaFileShowerGenerator::ReadLong()
    //
    // Read long file
    //
{

Msg(EsafMsg::Debug) << "Reading Long File.." << MsgDispatch;
                for ( i=0; i<10; i++) 
                { en[i]=0;
                  part[i]=0;
		}
    nshow=0;
    pos = ftell(fLongFile); 
    while (!feof(fLongFile) && nshow<=event)
    {
       fscanf(fLongFile,"%*s%*s%*s%f%*s%*s%*s%f%*s%*s%*s%i",&nstep,&step_long,&nshow);
       fgets(longitudinal,sizeof(longitudinal),fLongFile);
       pos=ftell(fLongFile);
       Msg(EsafMsg::Debug)<< "nstep, step_long:  " << nstep << "; t"<< step_long << MsgDispatch;
       fgets(longitudinal,sizeof(longitudinal),fLongFile);
                                                                                                                            
       if (nshow!=event)
         {
           for (j=1; j<= 2*nstep+12; j++ ){
             fgets(longitudinal,sizeof(longitudinal),fLongFile);
             pos = ftell(fLongFile);
             fseek(fLongFile,pos,SEEK_SET);
           }
           fseek(fLongFile,pos,SEEK_SET);
           nshow = nshow+1;
         }
                                                                                                                            
       else
         {
           for ( j=1; j <= nstep; j++)
          {
             for ( i=0; i<10; i++)
              {
               fscanf(fLongFile,"%f",&part[i]);
               pos=ftell(fLongFile);
               fseek(fLongFile,pos,SEEK_SET);
              } // reading 1 raw
             particle.insert(pair<int,float>(j,part[0]));
             particle.insert(pair<int,float>(j,part[0] - step_long));
             particle.insert(pair<int,float>(j,part[3] + part[2]));
             particle.insert(pair<int,float>(j,part[7]));
             particle.insert(pair<int,float>(j,part[9]));
         } // reading block
                                                                                                                            
             fgets(longitudinal,sizeof(longitudinal),fLongFile);
             fgets(longitudinal,sizeof(longitudinal),fLongFile);
             fgets(longitudinal,sizeof(longitudinal),fLongFile);
                                                                                                                            
           for ( j=1; j <=nstep; j++)
            {
              for ( i=0; i<10; i++)
               { fscanf(fLongFile,"%f",&en[i]); 
                } // reading 1 raw
              particle.insert(pair<int,float>(j,en[2]));
       	    } // reading block
           
	      fgets(longitudinal,sizeof(longitudinal),fLongFile);
	      fgets(longitudinal,sizeof(longitudinal),fLongFile);
	      fgets(longitudinal,sizeof(longitudinal),fLongFile);

	      pos=ftell(fLongFile);
              fseek(fLongFile,pos,SEEK_SET);
              fscanf(fLongFile,"%*s%*s%f%f%f%*f%*f%*f",&Nemax,&to,&tmax);
              for (int i = 1; i<=8; i++) {
	      fgets(longitudinal,sizeof(longitudinal),fLongFile);
		}

             nshow=nshow+1;
         } // end of else 
                                                                                                                            
       pos=ftell(fLongFile);
     }
                                                                                                                            

}


//______________________________________________________________________________
 bool CorsikaFileShowerGenerator::LoadHeader()
    //
    // Fill header object
    //
{
    fHeader.clear();
    ReadParticleHeader();

    x0comes = xi = 0;
    y0comes = yi = 0;
    if (z_firstint>=0) {
	    z0comes = z_firstint;}
    else {
 	z0comes = - z_firstint;}
    zmax = GetAltitude(tmax);

    vertical_distance = (abs(z_firstint) - zmax);
    coeff_f = GetCoeff(vertical_distance);
    
    xmax = GetXstep(coeff_f);
    ymax = GetYstep(coeff_f);
    float coeff = GetCoeff(z0comes);

    fHeader.insert(pair<string,float>("EPRMENRG",energy));
    fHeader.insert(pair<string,float>("EPRMDPTH",thick_firstint));
    fHeader.insert(pair<string,float>("EX0COMES",x0comes));  
    fHeader.insert(pair<string,float>("EY0COMES",y0comes)); 
    fHeader.insert(pair<string,float>("EZ0COMES",z0comes)); 
    fHeader.insert(pair<string,float>("EPOLANGD",epolangd)); 
    fHeader.insert(pair<string,float>("EAZIANGD",eaziangd));
    fHeader.insert(pair<string,float>("ELDIRCOS",eldircos)); 
    fHeader.insert(pair<string,float>("EMDIRCOS",emdircos));
    fHeader.insert(pair<string,float>("ENDIRCOS",endircos));
    fHeader.insert(pair<string,float>("EXDMPPOS",GetXstep(coeff)));
    fHeader.insert(pair<string,float>("EYDMPPOS",GetYstep(coeff)));
    fHeader.insert(pair<string,float>("EZDMPPOS",0.));
    fHeader.insert(pair<string,float>("EXMAXPOS",xmax));  
    fHeader.insert(pair<string,float>("EYMAXPOS",ymax)); 
    fHeader.insert(pair<string,float>("EZMAXPOS",zmax)); 
    return true;

}

//______________________________________________________________________________
 bool CorsikaFileShowerGenerator::LoadTrack(Int_t event)
    // 
    // Load one event from file
    //
{

    ReadLong();
    if (!LoadHeader()) return false;


    fTrack = new ShowerTrack();

    fTrack->fEthrEl = 0.;
    fTrack->fEnergy = fHeader["EPRMENRG"]*GeV;   //in eV
    fTrack->fTheta  = fHeader["EPOLANGD"]*deg;
    fTrack->fPhi    = fHeader["EAZIANGD"]*deg;
    fTrack->fX1     = fHeader["EPRMDPTH"]*g/cm2;
    
    fTrack->fDirVers[0] = fHeader["ELDIRCOS"];
    fTrack->fDirVers[1] = fHeader["EMDIRCOS"];
    fTrack->fDirVers[2] = fHeader["ENDIRCOS"];
    
    fTrack->fInitPos[0] = fHeader["EX0COMES"]*cm;
    fTrack->fInitPos[1] = fHeader["EY0COMES"]*cm;
    fTrack->fInitPos[2] = fHeader["EZ0COMES"]*cm;


    if(fHeader["HITGROUN"]==0) 
    {  
        fTrack->fHitGround = false;
    }
    else  
    {
        fTrack->fHitGround = true;
    }



    coeff_i = 0;
    xi=0;
    yi=0;
    zi = abs(z_firstint);
    timei = 0;
    for(j=1; j<=nstep; j++)
    {


        int count=0;
        for (pnum=particle.lower_bound(j); pnum!=particle.upper_bound(j); ++pnum)
        {
            step[count] = pnum->second;
            count++;
        }


// step_bound (g/cm^2)
        fStep.fXi = (step[1]+ thick_firstint)*gram/cm2;
        fStep.fXf = (step[1]+ thick_firstint+step_long)*gram/cm2;

// step_bound (cartesian coordinates) 
        
        thickness = step_long*endircos;
        zf = GetAltitude(step[1]+ thick_firstint+thickness);
                                                                               
        vertical_distance = abs(zi-zf);
        coeff_f = GetCoeff(vertical_distance);
        xf = GetXstep(coeff_f);
        yf = GetYstep(coeff_f);

       if (zi>0) {
        fStep.fXYZi[2] = zi*cm;}
        else {
 	fStep.fXYZi[2] = 0;}

        if (zf>0) {
        fStep.fXYZf[2] = zf*cm;}
        else {
 	fStep.fXYZf[2] = 0;}
 if (fStep.fXYZi[2]!=0 || fStep.fXYZf[2]!=0 ) 
{
//        fStep.fXYZf[2] = zf;
        fStep.fXYZi[0] = xi*cm;
        fStep.fXYZf[0] = xf*cm;
        fStep.fXYZi[1] = yi*cm;
        fStep.fXYZf[1] = yf*cm;
// time
        distance = sqrt(abs(pow(xi,2) + pow(yi,2) + pow(zi,2) - (pow(xf,2) + pow(yf,2) + pow(zf,2))));

        fStep.fTimei = timei*nanosecond; 
        timef = GetStepTime(distance) + timei;
        fStep.fTimef = timef*nanosecond;
        fStep.fNelectrons = step[2];
        fStep.fNcharged = step[3];
        fStep.fNcherenkov = step[4];
        fStep.fEloss = step[5];

// age 

        fStep.fAgei = GetAge(step[1] + thick_firstint);
        fStep.fAgef = GetAge(step[1] + thick_firstint+step_long);

        fTrack->Add(fStep);
        coeff_i = coeff_f;
        xi = xf;
        yi = yf;
        zi = zf;
        timei = timef;
} 


   }

    return true;

}

//______________________________________________________________________________
 MCTruth* CorsikaFileShowerGenerator::GetTruth() 
    // 
    // Get pointer to MonteCarlo Truth
    // 
{


    if (!fTrack->Size()) return NULL;
    if (!fTruth) {
         fTruth = new MCTruth();}

    fTruth->SetEnergy(fHeader["EPRMENRG"]*GeV);   //in eV
    fTruth->SetThetaPhi(fHeader["EPOLANGD"]*deg,fHeader["EAZIANGD"]*deg);
    EVector v;
    v[0] = fHeader["EX0COMES"]*cm;
    v[1] = fHeader["EY0COMES"]*cm;
    v[2] = fHeader["EZ0COMES"]*cm;
    fTruth->SetFirstInt(v, fHeader["EPRMDPTH"]);
    v[0] = fHeader["EXDMPPOS"]*cm;
    v[1] = fHeader["EYDMPPOS"]*cm;
    v[2] = fHeader["EZDMPPOS"]*cm;
                                                                                                                             
                                                                                                                             
    // in case no value for ShowerMax is present, put arbitrary number 0
                                                                                                                             
    v[0] = 0;
    v[1] = 0;
    v[2] = 0;
    map<string,double>::iterator showmax;
    showmax = fHeader.find(string("EXMAXPOS"));
    if (showmax != fHeader.end()) v[0] = fHeader["EXMAXPOS"]*cm;
    showmax = fHeader.find(string("EYMAXPOS"));
    if (showmax != fHeader.end()) v[1] = fHeader["EYMAXPOS"]*cm;
    showmax = fHeader.find(string("EZMAXPOS"));
    if (showmax != fHeader.end()) v[2] = fHeader["EZMAXPOS"]*cm;
                                                                                                                             
    fTruth->SetShowerMax(v,fHeader["ESHOWMAX"]);

  return fTruth;
                                                                                                                           
}

//______________________________________________________________________________
 float CorsikaFileShowerGenerator::GetAltitude(float thick)
{
float h;
 if (thick >= GetThickness(4.e+05)) {
	h = -catm[0]*log((thick - aatm[0])/batm[0]);}
 if (thick >= GetThickness(10.e+05) && thick < GetThickness(4.e+05)) {
	h = -catm[1]*log((thick - aatm[1])/batm[1]);}
 if (thick >= GetThickness(40.e+05) && thick< GetThickness(10.e+05)) {
	h = -catm[2]*log((thick - aatm[2])/batm[2]);}
 if (thick >= GetThickness(100.e+05) && thick< GetThickness(40.e+05)) {
	h = -catm[3]*log((thick - aatm[3])/batm[3]);}
 if (thick <= GetThickness(100.e+05) ) {
	h = (aatm[4]-thick)*catm[4]/batm[4];}

 return h;
}

//______________________________________________________________________________
 float CorsikaFileShowerGenerator::GetThickness(float z1)
{
 float thick_f;
 float z = abs(z1);
 if (z<= 4.e+05) { 

	thick_f = aatm[0]+batm[0]*exp(-z/catm[0]);}
 else if (z> 4.e+05 && z<= 10.e+05) {
        thick_f = aatm[1]+batm[1]*exp(-z/catm[1]);}
 else if (z> 10.e+05 && z<= 40.e+05) {
        thick_f = aatm[2]+batm[2]*exp(-z/catm[2]);}
 else if (z> 40.e+05 && z<= 100.e+05) {
        thick_f = aatm[3]+batm[3]*exp(-z/catm[3]);}
 else if (z> 100.e+05) {
	thick_f = aatm[4]-batm[4]*z/catm[4];}
 return thick_f;
}

//______________________________________________________________________________
 float CorsikaFileShowerGenerator::GetCoeff(float h)
{
 float coeff = h/endircos;
 return coeff;
}



//______________________________________________________________________________
 float CorsikaFileShowerGenerator::GetXstep(float t)
{
 float x = xi + t*eldircos; 
 return x;
}


//______________________________________________________________________________
 float CorsikaFileShowerGenerator::GetYstep(float t)
{
 float y = yi + t*emdircos; 
 return y;
}

//______________________________________________________________________________
 float CorsikaFileShowerGenerator::GetStepTime(float h)
{
 float time = 0.1*h/Clight();
 return time;
}


//______________________________________________________________________________
 float CorsikaFileShowerGenerator::GetAge(float h)
{
 float age;
    float t = h-thick_firstint;
    age = 3*t/(t + 2*tmax);
 return age;
}
                                                                                                                            

About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005