Subversion Repositories f9daq

Rev

Rev 268 | Rev 273 | Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | RSS feed

  1. #include <TClonesArray.h>
  2. #include <TH1F.h>
  3. #include <TStopwatch.h>
  4. #include <TDatime.h>
  5. #include <TString.h>
  6. #include <TFile.h>
  7. #include <TTree.h>
  8. #include <TBranch.h>
  9. #include <TBufferJSON.h>
  10.  
  11. #include "BParticle.h"
  12. #include "BEvent.h"
  13.  
  14. class Hdr{
  15. public:
  16.   int id;
  17.   int len;
  18.   int progress;
  19. };
  20.  
  21.  
  22. class Blab2 {
  23. public:
  24. const char *names[11]={"photon", "electron", "pion", "muon", "kaon", "proton", "J/Psi", "D", "D*", "B", "Phi"};
  25.  
  26. UInt_t fNeve;
  27. UInt_t fNfirst;
  28. UInt_t fPrint;
  29. int fData;
  30. TH1F *fH[100];
  31. UInt_t fHtype[100];
  32. TClonesArray *fList[100];
  33.  
  34. Blab2();
  35. ~Blab2();
  36. void Init();
  37. void event();
  38. void Process();
  39. void h1d(const char *varname, const char *name, int nbins, double min, double max, int id );
  40. int  selector(int pin, int charge, SIMPLEPID particlename,  int hid, int pout );
  41. int  combiner(int id0 ,int id1 , int same, SIMPLEPID particlename, double min, double max, int hid, int id );
  42. int  fix_mass(int id);
  43. int  Fill(int hid, BParticle *p);
  44. void plist(int i);
  45.  
  46.  
  47. ClassDef ( Blab2, 1 )
  48. };
  49.  
  50. ClassImp(Blab2)
  51.  
  52. Blab2::Blab2():fNfirst(0), fNeve(0), fData(0), fPrint(0) {
  53.  
  54.  Process();
  55. };
  56.  
  57.  
  58. void Blab2::h1d(const char *varname, const char *name, int nbins, double min, double max, int id ){
  59.    TString svar(varname);
  60.    TString axis[]={"mass (GeV/c2)","momentum (GeV/c)","energy (GeV)","charge","identity","momentum (GeV/c)","momentum (GeV/c)","momentum (GeV/c)","momentum (GeV/c)"};
  61.    fHtype[id] = 0;
  62.    if (svar.Contains("GetMass"    )) fHtype[id]=0;
  63.    if (svar.Contains("GetMomentum")) fHtype[id]=1;
  64.    if (svar.Contains("GetEnergy"  )) fHtype[id]=2;
  65.    if (svar.Contains("GetCharge"  )) fHtype[id]=3;
  66.    if (svar.Contains("GetPid"     )) fHtype[id]=4;
  67.    if (svar.Contains("GetXMomentum")) fHtype[id]=5;
  68.    if (svar.Contains("GetYMomentum")) fHtype[id]=6;
  69.    if (svar.Contains("GetZMomentum")) fHtype[id]=7;
  70.    if (svar.Contains("GetTranverseMomentum")) fHtype[id]=8;
  71.  
  72.    fH[id]= new TH1F(TString::Format("h%d",id), TString::Format("%s;%s;N",name,axis[fHtype[id]].Data()), nbins, min, max);
  73.    
  74.    if (fHtype[id]==4) {
  75.      fH[id]= new TH1F(TString::Format("h%d",id), TString::Format("%s;%s;N",name,axis[fHtype[id]].Data()), 11, 0, 11);
  76.      for (int i=0;i<11;i++) fH[id]->GetXaxis()->SetBinLabel(i+1,names[i]);
  77.    } else {
  78.      fH[id]= new TH1F(TString::Format("h%d",id), TString::Format("%s;%s;N",name,axis[fHtype[id]].Data()), nbins, min, max);  
  79.    }
  80.    
  81.  
  82. }
  83.  
  84. int Blab2::Fill(int hid, BParticle *p){
  85. if (hid>=0 && fH[hid]) {
  86.       float val;
  87.       switch (fHtype[hid]){
  88.       case 0: val = p->GetMass(); break;
  89.       case 1: val = p->GetMomentum(); break;
  90.       case 2: val = p->e(); break;
  91.       case 3: val = p->charge(); break;
  92.       case 4: val = p->pid(); break;
  93.       case 5: val = p->px(); break;
  94.       case 6: val = p->py(); break;
  95.       case 7: val = p->pz(); break;
  96.       case 8: val = sqrt(p->px()*p->px()+p->py()*p->py()); break;
  97.       default : val = 0 ; break;
  98.       }
  99.    fH[hid]->Fill(val);
  100. }  
  101.  
  102. return 0;
  103. }
  104.  
  105. int Blab2::combiner(int _p0, int _p1,int same, SIMPLEPID pid, double min, double max, int hid, int _p2 ){
  106.    // Loop over all the particles in both lists.
  107.  if (_p0 < 0 ) _p0 =0;
  108.  if (_p1 < 0 ) _p1 =0;
  109.  
  110.  
  111.  fList[_p2]->Clear();
  112.  int nprt=0;
  113.  
  114.  for(TIter next1(fList[_p0]);BParticle * p1 =(BParticle *) next1();) {
  115.     // the second loop
  116.    // in the case the second parti
  117.    for(TIter next2 = (_p0!=_p1 && same==0) ?  TIter(fList[_p1]): TIter(next1) ; BParticle * p2 =(BParticle *) next2();) {  
  118.       if (p1==p2) continue;     // do not use the same particle in the combinations
  119.       BParticle  p = *p1 + *p2; // Combine two particles into a new particle  
  120.       if (p.InMassRange(min, max)){
  121.             Fill(hid, &p);
  122.             p.SetPid(pid); // set PID to particlename to fix the particle mass
  123.             p.SetEnergyFromPid();
  124.             TClonesArray& list = *fList[_p2];          
  125.             new (list[nprt++]) BParticle ( p ); // create a new entry in kslist list of particles
  126.            
  127.       }
  128.        
  129.    }
  130.                
  131.  }
  132.  return _p2;
  133. }
  134.  
  135.  
  136. int Blab2::selector(int pin, int charge, SIMPLEPID type ,  int hid, int pout ){
  137.  if (pin < 0 ) pin =0;
  138.  
  139.   fList[pout]->Clear();
  140.   int nprt=0;
  141.  
  142.   for(TIter next(fList[pin]); BParticle * p =(BParticle *) next();) {
  143.         if (p->charge()== charge || charge > 1){
  144.           if ( p->pid()== type || type == ALL ) {
  145.             TClonesArray& list = *fList[pout];
  146.             new (list[nprt++]) BParticle ( *p );
  147.             Fill(hid, p);
  148.           }
  149.         }
  150.   }      
  151.    return pout;
  152. }
  153.  
  154.  
  155. int Blab2::fix_mass(int pin){
  156.    if (pin < 0 ) pin =0;
  157.    for(TIter next(fList[pin]); BParticle * p =(BParticle *) next();)  p->SetEnergyFromPid();
  158.    return pin;
  159. }
  160.  
  161. void Blab2::plist(int i){
  162.   fList[i]= new TClonesArray( "BParticle", 500 );
  163. }
  164. Blab2::~Blab2(){};
  165.  
  166.  
  167. void send_message(int id, TString msg, int progress){
  168. static Hdr hdr;
  169.  
  170.    hdr.id = id;
  171.    hdr.len= msg.Length();
  172.    hdr.progress= progress;
  173.    fwrite(&hdr,3*sizeof(int),1,stdout);
  174.    fwrite(msg.Data(),hdr.len,1,stdout);
  175.    fflush(stdout);
  176.  
  177. }
  178.  
  179.  
  180. void Blab2::Process(){
  181.  
  182. char sList[0xFFFF];
  183. for (int i=0;i<100;i++) fH[i]=NULL;
  184. for (int i=0;i<100;i++) fHtype[i]=0;
  185. for (int i=0;i<100;i++) fList[i]=NULL;
  186.  
  187. Init();
  188.  
  189. TFile * f = new TFile(TString::Format("../../data/hadron-%d.root",fData)); // Open a data file
  190. if(f->IsZombie()) {  send_message(0,TString::Format("File %d not found\n",fData), 0 );  return; }  
  191. TTree * t =(TTree *) f-> Get( "T"); // Obtain a pointer to a tree of "event" data in the file
  192. BEvent * mevent = new BEvent(); // Create a  "BEvent" object where data from the file will be loaded
  193. TBranch * branch = t-> GetBranch( "BEvent"); // Obtain a branch for "BEvent" in the tree
  194. branch-> SetAddress(&mevent); // Register created "BEvent" object to the branch
  195. TH1F *fHnprt= new TH1F("h100", "Number of particles in the event;N particles;N events", 50, -0.5, 49.5);
  196.  
  197.  
  198.  
  199. send_message(0, TString::Format("<br>Number of Events in the file %lld<br>\n", t->GetEntries() ),0);
  200. TStopwatch timer;
  201. timer.Start();
  202. int nev  = 0;
  203. int i    =TMath::Min(fNfirst, (UInt_t) t-> GetEntries());
  204. int cNeve=TMath::Min(fNfirst+fNeve, (UInt_t) t-> GetEntries());
  205. while (i<cNeve){
  206.  t-> GetEntry(i); // Read the content of the event from the file
  207.  fList[0]= mevent->GetParticleList();
  208.  
  209.  event();
  210.  
  211.  int progress = (100*i)/cNeve;
  212.  if (i%10000==0) send_message(2, TString::Format("Event %d\n",i), progress);
  213.  
  214.  int nprt=0;
  215.  if (nev>100) fPrint = 0; // disable particle prints for huge numer of events
  216.  if (fPrint) sprintf(sList,"Primary particle list for Event %d<br/><table class='plist' ><tr><th>N<th>px(GeV/c)<th>py(GeV/c)<th>pz(GeV/c)<th>p(GeV/c)<th>Energy(GeV)<th>Charge<th>ID<th></tr>", i);
  217.  for(TIter next(fList[0]); BParticle * p =(BParticle *) next();) {
  218.    nprt++;
  219.    if (fPrint) sprintf(sList,"%s<tr><td>%d<td>%g<td>%g<td>%g<td>%g<td>%g<td>%1.0f<td>%s</tr>",sList,nprt, p->px(),p->py(),p->pz(),p->GetMomentum(),p->e(), p->charge(),names[p->pid()] );
  220.  }
  221.  fHnprt->Fill(nprt);
  222.  if (fPrint) {
  223.    sprintf(sList,"%s</table>",sList);
  224.    send_message(0, TString(sList),progress);
  225.    nev++;
  226.  }
  227.  mevent-> Clear();  // Clear the memory.
  228.  for (int k=0;k<100;k++) if (fList[k]!=0) fList[k]->Clear();
  229.  i++;
  230. }
  231.  
  232. send_message(2, TString::Format("Number of events processed %d\n", i ),100);
  233.  
  234. send_message(1,TBufferJSON::ConvertToJSON(fHnprt),100 );
  235.  
  236.  
  237. for (int i=0;i<100;i++) if (fH[i]!=0) send_message(1,TBufferJSON::ConvertToJSON(fH[i]),100 );
  238.  
  239. TDatime d;
  240. timer.Stop();
  241. send_message(3, TString::Format("'%s', %d, %f, %f", d.AsSQLString(),i, timer.RealTime(), timer.CpuTime() ),100);
  242.  
  243. }
  244.  
  245.  
  246.