Subversion Repositories f9daq

#include <TClonesArray.h>
#include <TH1F.h>
#include <TStopwatch.h>
#include <TDatime.h>
#include <TString.h>
#include <TFile.h>
#include <TTree.h>
#include <TBranch.h>
#include <TClonesArray.h>

#include "BParticle.h"
#include "BEvent.h"

class Hdr{
public:
  int id;
  int len;
  int progress;
};


class Blab2 {
public:
UInt_t fNeve;
UInt_t fNfirst;
int fData;
TH1F *fH[100];
UInt_t fHtype[100];
TClonesArray *fList[100];

Blab2();
~Blab2();
void Init();
void event();
void Process();
void h1d(const char *varname, const char *name, int nbins, double min, double max, int id );
int  selector(int pin, int charge, SIMPLEPID particlename,  int hid, int pout );
int  combiner(int id0 ,int id1 , SIMPLEPID particlename, double min, double max, int hid, int id );
int  fix_mass(int id);
int  Fill(int hid, BParticle *p);
void plist(int i);


ClassDef ( Blab2, 1 )
};

ClassImp(Blab2)

Blab2::Blab2():fNfirst(0), fNeve(0), fData(0){

 Process();
};


void Blab2::h1d(const char *varname, const char *name, int nbins, double min, double max, int id ){
   TString svar(varname);
   TString axis[3]={"mass (GeV/c2)","momentum (GeV/c)","energy (GeV)"};
   fHtype[id] = 0;
   if (svar.Contains("GetMass"    )) fHtype[id]=0;
   if (svar.Contains("GetMomentum")) fHtype[id]=1;
   if (svar.Contains("GetEnergy"  )) fHtype[id]=2;

   fH[id]= new TH1F(TString::Format("h%d",id), TString::Format("%s;%s;N",name,axis[fHtype[id]].Data()), nbins, min, max);
}

int Blab2::Fill(int hid, BParticle *p){
if (hid>=0 && fH[hid]) {
      float val;
      switch (fHtype[hid]){
      case 0: val = p->GetMass(); break;
      case 1: val = p->GetMomentum(); break;
      case 2: val = p->e(); break;
      default : val = 0 ; break;
      }
   fH[hid]->Fill(val);
}  

return 0;
}

int Blab2::combiner(int _p0, int _p1,SIMPLEPID pid, double min, double max, int hid, int _p2 ){
   // Loop over all the particles in both lists.
 if (_p0 < 0 ) _p0 =0;
 if (_p1 < 0 ) _p1 =0;
 

 fList[_p2]->Clear();
 int nprt=0;
 
 for(TIter next1(fList[_p0]);BParticle * p1 =(BParticle *) next1();) {
    // the second loop
   
   for(TIter next2 = (_p0!=_p1) ?  TIter(fList[_p1]): TIter(next1) ; BParticle * p2 =(BParticle *) next2();) {  
      if (p1==p2) continue;
      BParticle  p = *p1 + *p2; // Combine two particles into a new particle  
      if (p.InMassRange(min, max)){
            //SIMPLEPID pid = PION; // should be set to particlename
            p.SetPid(pid);
            TClonesArray& list = *fList[_p2];
            new (list[nprt++]) BParticle ( p ); // create a new entry in kslist list of particles
            Fill(hid, &p);
      }
       
   }
               
 }
 return _p2;
}


int Blab2::selector(int pin, int charge, SIMPLEPID type ,  int hid, int pout ){
 if (pin < 0 ) pin =0;
 
  fList[pout]->Clear();
  int nprt=0;
 
  for(TIter next(fList[pin]); BParticle * p =(BParticle *) next();) {
        if (p->charge()== charge && ( p->pid()== type || type == ALL )) {
          TClonesArray& list = *fList[pout];
          new (list[nprt++]) BParticle ( *p );
          Fill(hid, p);
        }
  }      
   return pout;
}


int Blab2::fix_mass(int pin){
   if (pin < 0 ) pin =0;
   for(TIter next(fList[pin]); BParticle * p =(BParticle *) next();)  p->SetEnergyFromPid();
   return pin;
}

void Blab2::plist(int i){
  fList[i]= new TClonesArray( "BParticle", 500 );
}
Blab2::~Blab2(){};


void send_message(int id, TString msg, int progress){
static Hdr hdr;

   hdr.id = id;
   hdr.len= msg.Length();
   hdr.progress= progress;
   fwrite(&hdr,3*sizeof(int),1,stdout);
   fwrite(msg.Data(),hdr.len,1,stdout);
   fflush(stdout);

}


void Blab2::Process(){

for (int i=0;i<100;i++) fH[i]=NULL;
for (int i=0;i<100;i++) fHtype[i]=0;
for (int i=0;i<100;i++) fList[i]=NULL;

Init();

TFile * f = new TFile(TString::Format("../../data/hadron-%d.root",fData)); // Open a data file
if(f->IsZombie()) {  send_message(0,TString::Format("File %d not found\n",fData), 0 );  return; }  
TTree * t =(TTree *) f-> Get( "T"); // Obtain a pointer to a tree of "event" data in the file
BEvent * mevent = new BEvent(); // Create a  "BEvent" object where data from the file will be loaded
TBranch * branch = t-> GetBranch( "BEvent"); // Obtain a branch for "BEvent" in the tree
branch-> SetAddress(&mevent); // Register created "BEvent" object to the branch
TH1F *fHnprt= new TH1F("h100", "Number of particles in the event;N particles;N events", 50, -0.5, 49.5);
TH1F *fHid= new TH1F("h101", "Particle type;ID;N particles", 6, -0, 6);

const char *names[6]={"PHOTON", "ELECTRON", "PION", "MUON", "KAON", "PROTON"};
for (int i=1;i<=6;i++) fHid->GetXaxis()->SetBinLabel(i,names[i-1]);

send_message(0, TString::Format("Number of Entries in the file %lld\n", t->GetEntries() ),0);
TStopwatch timer;
timer.Start();
int i    =TMath::Min(fNfirst, (UInt_t) t-> GetEntries());
int cNeve=TMath::Min(fNfirst+fNeve, (UInt_t) t-> GetEntries());
while (i<cNeve){
 t-> GetEntry(i); // Read the content of the event from the file
 fList[0]= mevent->GetParticleList();
 
 event();
 
 if (i%10000==0) send_message(2, TString::Format("Event %d\n",i), (100*i)/cNeve);

 int nprt=0;
 for(TIter next(fList[0]); BParticle * p =(BParticle *) next();) {
   nprt++;
   fHid->Fill(p->pid());
 }
 fHnprt->Fill(nprt);

 mevent-> Clear();  // Clear the memory.
 for (int k=0;k<100;k++) if (fList[k]!=0) fList[k]->Clear();
 i++;
}

send_message(2, TString::Format("Number of events processed %d\n", i ),100);
for (int i=0;i<100;i++) if (fH[i]!=0) send_message(1,TBufferJSON::ConvertToJSON(fH[i]),100 );
send_message(1,TBufferJSON::ConvertToJSON(fHnprt),100 );
send_message(1,TBufferJSON::ConvertToJSON(fHid),100 );
TDatime d;
timer.Stop();
send_message(3, TString::Format("'%s', %d, %f, %f", d.AsSQLString(),i, timer.RealTime(), timer.CpuTime() ),100);

}