Subversion Repositories f9daq

Rev

Rev 331 | 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. #include <TMath.h>
  11. #include <vector>
  12. #include "BParticle.h"
  13. #include "BEvent.h"
  14.  
  15. class Hdr{
  16. public:
  17.   int id;
  18.   int len;
  19.   int progress;
  20. };
  21.  
  22. std::vector<int> histogram(int n, ...){
  23. std::vector<int> result;
  24. int val = 0;
  25.    va_list ap;
  26.    int i;
  27.  
  28.    va_start(ap, n);
  29.    for(i = 0; i < n; i++) {
  30.       result.push_back(  va_arg(ap, int) );
  31.    }
  32.    va_end(ap);
  33. return result;
  34. }
  35.  
  36. class Blab2 {
  37. public:
  38. const char *names[12]={"photon", "electron", "pion", "muon", "kaon", "proton", "J/Psi", "D", "D*", "B", "Phi","Lambda0"};
  39.  
  40. UInt_t fNeve;
  41. UInt_t fNfirst;
  42. UInt_t fPrint;
  43. int fData;
  44. TH1F *fH[100];
  45. UInt_t fHtype[100];
  46. TClonesArray *fList[100];
  47.  
  48. Blab2();
  49. ~Blab2();
  50. void Init();
  51. void event();
  52. void Process();
  53. void h1d(const char *varname, const char *name, int nbins, double min, double max, int id );
  54. int  selector(int pin, int charge, SIMPLEPID particlename,  int hid, int pout );
  55. int  selector(int pin, int charge, SIMPLEPID particlename,  std::vector<int> hid, int pout );
  56. int  combiner(int id0 ,int id1 , int same, SIMPLEPID particlename, double min, double max, int hid, int id );
  57. int  combiner(int id0 ,int id1 , int same, SIMPLEPID particlename, double min, double max, std::vector<int>hid, int id );
  58. int  combiner3(int id0 ,int id1 , int id2, int same, SIMPLEPID particlename, double min, double max, int hid, int id );
  59. int  combiner3(int id0 ,int id1 , int id2, int same, SIMPLEPID particlename, double min, double max, std::vector<int>hid, int id );
  60. int  fix_mass(int id);
  61. int  Fill(std::vector<int> hid, std::vector<BParticle *> p);
  62. void plist(int i);
  63.  
  64.  
  65. ClassDef ( Blab2, 1 )
  66. };
  67.  
  68. ClassImp(Blab2)
  69.  
  70. Blab2::Blab2():fNfirst(0), fNeve(0), fData(0), fPrint(0) {
  71.  
  72.  Process();
  73. };
  74.  
  75.  
  76. void Blab2::h1d(const char *varname, const char *name, int nbins, double min, double max, int id ){
  77.    TString svar(varname);
  78.    TString axis[]={"mass (GeV/c2)",
  79.                 "momentum (GeV/c)",
  80.                 "energy (GeV)","charge",
  81.                 "identity",
  82.                 "momentum (GeV/c)",
  83.                 "momentum (GeV/c)",
  84.                 "momentum (GeV/c)",
  85.                 "momentum (GeV/c)",
  86.                 "angle (deg.)",
  87.                 "cos(theta)", "#Deltamass_{1}(GeV/c2)", "#Deltamass_{2}(GeV/c2)", "#Delta mass_{3}(GeV/c2)"};
  88.    fHtype[id] = 0;
  89.    if (svar.Contains("GetMass"    )) fHtype[id]=0;
  90.    if (svar.Contains("GetMomentum")) fHtype[id]=1;
  91.    if (svar.Contains("GetEnergy"  )) fHtype[id]=2;
  92.    if (svar.Contains("GetCharge"  )) fHtype[id]=3;
  93.    if (svar.Contains("GetPid"     )) fHtype[id]=4;
  94.    if (svar.Contains("GetXMomentum")) fHtype[id]=5;
  95.    if (svar.Contains("GetYMomentum")) fHtype[id]=6;
  96.    if (svar.Contains("GetZMomentum")) fHtype[id]=7;
  97.    if (svar.Contains("GetTransverseMomentum")) fHtype[id]=8;
  98.    if (svar.Contains("GetTheta"))             fHtype[id]=9;
  99.    if (svar.Contains("GetCosTheta"))          fHtype[id]=10;
  100.    if (svar.Contains("GetDeltaMass1"    )) fHtype[id]=11;
  101.    if (svar.Contains("GetDeltaMass2"    )) fHtype[id]=12;
  102.    if (svar.Contains("GetDeltaMass3"    )) fHtype[id]=13;
  103.  
  104.    //fH[id]= new TH1F(TString::Format("h%d",id), TString::Format("%s;%s;N",name,axis[fHtype[id]].Data()), nbins, min, max);
  105.    if (fHtype[id]==4) {
  106.      fH[id]= new TH1F(TString::Format("h%d",id), TString::Format("%s;%s;N",name,axis[fHtype[id]].Data()), 11, 0, 11);
  107.      for (int i=0;i<11;i++) fH[id]->GetXaxis()->SetBinLabel(i+1,names[i]);
  108.    } else {
  109.      fH[id]= new TH1F(TString::Format("h%d",id), TString::Format("%s;%s;N",name,axis[fHtype[id]].Data()), nbins, min, max);  
  110.    }
  111.    
  112.  
  113. }
  114.  
  115.  
  116.  
  117. int Blab2::Fill(std::vector<int> id, std::vector<BParticle *> particle){
  118.   for (int i=0; i< id.size(); i++){
  119.   int hid = id[i];
  120.   BParticle *p= particle[0];
  121.   if (hid>=0 && fH[hid]) {
  122.       double val;
  123.       switch (fHtype[hid]){
  124.       case 0 : val  = p->GetMass(); break;
  125.       case 1 : val  = p->GetMomentum(); break;
  126.       case 2 : val  = p->e(); break;
  127.       case 3 : val  = p->charge(); break;
  128.       case 4 : val  = p->pid(); break;
  129.       case 5 : val  = p->px(); break;
  130.       case 6 : val  = p->py(); break;
  131.       case 7 : val  = p->pz(); break;
  132.       case 8 : val  = p->GetTransverseMomentum(); break;
  133.       case 9 : val  = (p->GetMomentum()!=0) ? p->pz()/p->GetMomentum() : 0; val = 180.0*TMath::ACos(val)/TMath::Pi(); break;
  134.       case 10: val  = (p->GetMomentum()!=0) ? p->pz()/p->GetMomentum() : 0; break;
  135.       case 11 :
  136.       case 12 :
  137.       case 13 : {
  138.         int idx = fHtype[hid]-10;
  139.         val = (particle.size()>idx) ? p->GetMass() - particle[idx]->GetMass() : 0; break;
  140.       }
  141.       default: val  = 0 ; break;
  142.    }
  143.    fH[hid]->Fill(val);
  144. }  
  145.    
  146.    }
  147.  
  148. return 0;
  149. }
  150. int Blab2::combiner(int _p0, int _p1,int same, SIMPLEPID pid, double min, double max, int hid, int _p2 ){
  151. std::vector<int> a;
  152. return combiner(_p0,_p1,same,pid,min,max,a,_p2);
  153. }
  154.  
  155. int Blab2::combiner3(int _p0, int _p1, int _p2, int same, SIMPLEPID pid, double min, double max, int hid, int _p3 ){
  156. std::vector<int> a;
  157. return combiner3(_p0,_p1,_p2, same,pid,min,max,a,_p3);
  158. }
  159.  
  160. int Blab2::combiner(int _p0, int _p1,int same, SIMPLEPID pid, double min, double max, std::vector<int> hid, int _p2 ){
  161.    // Loop over all the particles in both lists.
  162.  if (_p0 < 0 ) _p0 =0;
  163.  if (_p1 < 0 ) _p1 =0;
  164.  
  165.  
  166.  fList[_p2]->Clear();
  167.  std::vector<BParticle *> pall(3);
  168.  int nprt=0;
  169.  
  170.  for(TIter next1(fList[_p0]);BParticle * p1 =(BParticle *) next1();) {
  171.     // the second loop
  172.    // in the case the second parti
  173.    for(TIter next2 = (_p0!=_p1 && same==0) ?  TIter(fList[_p1]): TIter(next1) ; BParticle * p2 =(BParticle *) next2();) {  
  174.       if (p1==p2) continue;     // do not use the same particle in the combinations
  175.       BParticle  p = *p1 + *p2; // Combine two particles into a new particle
  176.       pall[0] = &p;
  177.       pall[1] = p1;
  178.       pall[2] = p2;
  179.  
  180.       if (p.InMassRange(min, max)){
  181.             Fill(hid, pall);
  182.             p.SetPid(pid); // set PID to particlename to fix the particle mass
  183.             p.SetEnergyFromPid();
  184.             TClonesArray& list = *fList[_p2];          
  185.             new (list[nprt++]) BParticle ( p ); // create a new entry in kslist list of particles
  186.            
  187.       }
  188.        
  189.    }
  190.                
  191.  }
  192.  return _p2;
  193. }
  194.  
  195.  
  196. int Blab2::combiner3(int _p0, int _p1,int _p2, int same, SIMPLEPID pid, double min, double max, std::vector<int> hid, int _p3 ){
  197.    // Loop over all the particles in both lists.
  198.  if (_p0 < 0 ) _p0 =0;
  199.  if (_p1 < 0 ) _p1 =0;
  200.  if (_p2 < 0 ) _p2 =0;
  201.  
  202.  
  203.  fList[_p3]->Clear();
  204.  std::vector<BParticle *> pall(4);  
  205.  
  206.  int nprt=0;
  207.  
  208.  for(TIter next1(fList[_p0]);BParticle * p1 =(BParticle *) next1();) {
  209.     // the second loop
  210.    // in the case the second parti
  211.    for(TIter next2 = (_p0!=_p1 && same==0) ?  TIter(fList[_p1]): TIter(next1) ; BParticle * p2 =(BParticle *) next2();) {  
  212.       if (p1==p2) continue;     // do not use the same particle in the combinations
  213.       for(TIter next3 = (_p1!=_p2 && same==0) ?  TIter(fList[_p2]): TIter(next2) ; BParticle * p3 =(BParticle *) next3();) {  
  214.         if (p2==p3) continue;     // do not use the same particle in the combinations
  215.         BParticle  p =  *p1 + *p2 + *p3 ; // Combine two particles into a new particle  
  216.         pall[0] = &p;
  217.         pall[1] = p1;
  218.         pall[2] = p2;
  219.         pall[3] = p3;
  220.  
  221.            
  222.         if (p.InMassRange(min, max)){
  223.             Fill(hid, pall);
  224.             p.SetPid(pid); // set PID to particlename to fix the particle mass
  225.             p.SetEnergyFromPid();
  226.             TClonesArray& list = *fList[_p3];          
  227.             new (list[nprt++]) BParticle ( p ); // create a new entry in kslist list of particles      
  228.         }
  229.       }
  230.        
  231.    }
  232.                
  233.  }
  234.  return _p3;
  235. }
  236.  
  237.  
  238. int Blab2::selector(int pin, int charge, SIMPLEPID type ,  int  hid, int pout ){
  239. std::vector<int> a;
  240. return selector(pin,charge,type,a,pout);
  241. }
  242. int Blab2::selector(int pin, int charge, SIMPLEPID type ,  std::vector<int> hid, int pout ){
  243.  if (pin < 0 ) pin =0;
  244.  
  245.   fList[pout]->Clear();
  246.   int nprt=0;
  247.   std::vector<BParticle *> pall(1);
  248.  
  249.   for(TIter next(fList[pin]); BParticle * p =(BParticle *) next();) {
  250.         if (p->charge()== charge || charge > 1){
  251.           if ( p->pid()== type || type == ALL ) {
  252.             TClonesArray& list = *fList[pout];
  253.             new (list[nprt++]) BParticle ( *p );
  254.             pall[0] = p ;  
  255.             Fill(hid, pall);
  256.           }
  257.         }
  258.   }      
  259.    return pout;
  260. }
  261.  
  262.  
  263. int Blab2::fix_mass(int pin){
  264.    if (pin < 0 ) pin =0;
  265.    for(TIter next(fList[pin]); BParticle * p =(BParticle *) next();)  p->SetEnergyFromPid();
  266.    return pin;
  267. }
  268.  
  269. void Blab2::plist(int i){
  270.   fList[i]= new TClonesArray( "BParticle", 500 );
  271. }
  272. Blab2::~Blab2(){};
  273.  
  274.  
  275. void send_message(int id, TString msg, int progress){
  276. static Hdr hdr;
  277.  
  278.    hdr.id = id;
  279.    hdr.len= msg.Length();
  280.    hdr.progress= progress;
  281.    fwrite(&hdr,3*sizeof(int),1,stdout);
  282.    fwrite(msg.Data(),hdr.len,1,stdout);
  283.    fflush(stdout);
  284.  
  285. }
  286.  
  287.  
  288. void Blab2::Process(){
  289.  
  290. char sList[0xFFFF];
  291. for (int i=0;i<100;i++) fH[i]=NULL;
  292. for (int i=0;i<100;i++) fHtype[i]=0;
  293. for (int i=0;i<100;i++) fList[i]=NULL;
  294.  
  295. Init();
  296.  
  297. TFile * f = new TFile(TString::Format("../../data/hadron-%d.root",fData)); // Open a data file
  298. if(f->IsZombie()) {  send_message(0,TString::Format("File %d not found\n",fData), 0 );  return; }  
  299. TTree * t =(TTree *) f-> Get( "T"); // Obtain a pointer to a tree of "event" data in the file
  300. BEvent * mevent = new BEvent(); // Create a  "BEvent" object where data from the file will be loaded
  301. TBranch * branch = t-> GetBranch( "BEvent"); // Obtain a branch for "BEvent" in the tree
  302. branch-> SetAddress(&mevent); // Register created "BEvent" object to the branch
  303. TH1F *fHnprt= new TH1F("h100", "Number of particles in the event;N particles;N events", 50, -0.5, 49.5);
  304.  
  305.  
  306.  
  307. send_message(0, TString::Format("<br>Number of Events in the file %lld<br>\n", t->GetEntries() ),0);
  308. TStopwatch timer;
  309. timer.Start();
  310. int nev  = 0;
  311. int i    =TMath::Min(fNfirst, (UInt_t) t-> GetEntries());
  312. int cNeve=TMath::Min(fNfirst+fNeve, (UInt_t) t-> GetEntries());
  313. int fPart = fPrint;
  314. int totaltracks = 0;
  315. while (i<cNeve){
  316.  t-> GetEntry(i); // Read the content of the event from the file
  317.  fList[0]= mevent->GetParticleList();
  318.  
  319.  event();
  320.  
  321.  int progress = (100*i)/cNeve;
  322.  if (i%10000==0) send_message(2, TString::Format("Event %d\n",i), progress);
  323.  
  324.  int nprt=0;
  325.  if (nev>100) fPrint = 0; // disable particle prints for huge numer of events
  326.  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);
  327.  for(TIter next(fList[0]); BParticle * p =(BParticle *) next();) {
  328.    nprt++;
  329.    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()] );
  330.  }
  331.  if (fPart) fHnprt->Fill(nprt);
  332.  totaltracks += nprt;
  333.  if (fPrint) {
  334.    sprintf(sList,"%s</table>",sList);
  335.    send_message(0, TString(sList),progress);
  336.    nev++;
  337.  }
  338.  mevent-> Clear();  // Clear the memory.
  339.  for (int k=0;k<100;k++) if (fList[k]!=0) fList[k]->Clear();
  340.  i++;
  341. }
  342. double avgtracks=(i)?float(totaltracks)/i:0;
  343. send_message(0, TString::Format("Number of events processed: %d<br/>\nNumber of particles: %d<br/>\nAverage number of particles per event %f<br/>\n", i, totaltracks, avgtracks ),100);
  344.  
  345. if (fPart) send_message(1,TBufferJSON::ConvertToJSON(fHnprt),100 );
  346.  
  347.  
  348. for (int i=0;i<100;i++) if (fH[i]!=0) send_message(1,TBufferJSON::ConvertToJSON(fH[i]),100 );
  349.  
  350. TDatime d;
  351. timer.Stop();
  352. send_message(3, TString::Format("'%s', %d, %f, %f", d.AsSQLString(),i, timer.RealTime(), timer.CpuTime() ),100);
  353.  
  354. }
  355.  
  356.  
  357.