WebSVN – f9daq – Diff – /belle2/masterclass/src/Blab2.cc


#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 <TBufferJSON.h>
#include <TMath.h>
#include <vector>
#include "BParticle.h"
#include "BEvent.h"
 
class Hdr{
public:
  int id;
  int len;
  int progress;
};
 
std::vector<int> histogram(int n, ...){
std::vector<int> result;
int val = 0;
   va_list ap;
   int i;
 
   va_start(ap, n);
   for(i = 0; i < n; i++) {
      result.push_back(  va_arg(ap, int) );
   }
   va_end(ap);
return result;
}
 
class Blab2 {
public:
const char *names[12]={"photon", "electron", "pion", "muon", "kaon", "proton", "J/Psi", "D", "D*", "B", "Phi","Lambda0"};
 
UInt_t fNeve;
UInt_t fNfirst;
UInt_t fPrint;
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  selector(int pin, int charge, SIMPLEPID particlename,  std::vector<int> hid, int pout );
int  combiner(int id0 ,int id1 , int same, SIMPLEPID particlename, double min, double max, int hid, int id );
int  combiner(int id0 ,int id1 , int same, SIMPLEPID particlename, double min, double max, std::vector<int>hid, int id );
int  combiner3(int id0 ,int id1 , int id2, int same, SIMPLEPID particlename, double min, double max, int hid, int id );
int  combiner3(int id0 ,int id1 , int id2, int same, SIMPLEPID particlename, double min, double max, std::vector<int>hid, int id );
int  fix_mass(int id);
int  Fill(std::vector<int> hid, std::vector<BParticle *> p);
void plist(int i);
 
 
ClassDef ( Blab2, 1 )
};
 
ClassImp(Blab2)
 
Blab2::Blab2():fNfirst(0), fNeve(0), fData(0), fPrint(0) {
 
 Process();
};
 
 
void Blab2::h1d(const char *varname, const char *name, int nbins, double min, double max, int id ){
   TString svar(varname);
   TString axis[]={"mass (GeV/c2)",
                "momentum (GeV/c)",
                "energy (GeV)","charge",
                "identity",
                "momentum (GeV/c)",
                "momentum (GeV/c)",
                "momentum (GeV/c)",
                "momentum (GeV/c)",
                "angle (deg.)",
                "cos(theta)", "#Deltamass_{1}(GeV/c2)", "#Deltamass_{2}(GeV/c2)", "#Delta mass_{3}(GeV/c2)"};
   fHtype[id] = 0;
   if (svar.Contains("GetMass"    )) fHtype[id]=0;
   if (svar.Contains("GetMomentum")) fHtype[id]=1;
   if (svar.Contains("GetEnergy"  )) fHtype[id]=2;
   if (svar.Contains("GetCharge"  )) fHtype[id]=3;
   if (svar.Contains("GetPid"     )) fHtype[id]=4;
   if (svar.Contains("GetXMomentum")) fHtype[id]=5;
   if (svar.Contains("GetYMomentum")) fHtype[id]=6;
   if (svar.Contains("GetZMomentum")) fHtype[id]=7;
   if (svar.Contains("GetTransverseMomentum")) fHtype[id]=8;
   if (svar.Contains("GetTheta"))             fHtype[id]=9;
   if (svar.Contains("GetCosTheta"))          fHtype[id]=10;
   if (svar.Contains("GetDeltaMass1"    )) fHtype[id]=11;
   if (svar.Contains("GetDeltaMass2"    )) fHtype[id]=12;
   if (svar.Contains("GetDeltaMass3"    )) fHtype[id]=13;
 
   //fH[id]= new TH1F(TString::Format("h%d",id), TString::Format("%s;%s;N",name,axis[fHtype[id]].Data()), nbins, min, max);
   if (fHtype[id]==4) {
     fH[id]= new TH1F(TString::Format("h%d",id), TString::Format("%s;%s;N",name,axis[fHtype[id]].Data()), 11, 0, 11);
     for (int i=0;i<11;i++) fH[id]->GetXaxis()->SetBinLabel(i+1,names[i]);
   } else {
     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(std::vector<int> id, std::vector<BParticle *> particle){
  for (int i=0; i< id.size(); i++){
  int hid = id[i];
  BParticle *p= particle[0];
  if (hid>=0 && fH[hid]) {
      double val;
      switch (fHtype[hid]){
      case 0 : val  = p->GetMass(); break;
      case 1 : val  = p->GetMomentum(); break;
      case 2 : val  = p->e(); break;
      case 3 : val  = p->charge(); break;
      case 4 : val  = p->pid(); break;
      case 5 : val  = p->px(); break;
      case 6 : val  = p->py(); break;
      case 7 : val  = p->pz(); break;
      case 8 : val  = p->GetTransverseMomentum(); break;
      case 9 : val  = (p->GetMomentum()!=0) ? p->pz()/p->GetMomentum() : 0; val = 180.0*TMath::ACos(val)/TMath::Pi(); break;
      case 10: val  = (p->GetMomentum()!=0) ? p->pz()/p->GetMomentum() : 0; break;
      case 11 :
      case 12 :
      case 13 : {
        int idx = fHtype[hid]-10;
        val = (particle.size()>idx) ? p->GetMass() - particle[idx]->GetMass() : 0; break;
      }
      default: val  = 0 ; break;
   }
   fH[hid]->Fill(val);
}  
   
   }
 
return 0;
}
int Blab2::combiner(int _p0, int _p1,int same, SIMPLEPID pid, double min, double max, int hid, int _p2 ){
std::vector<int> a;
return combiner(_p0,_p1,same,pid,min,max,a,_p2);
}
 
int Blab2::combiner3(int _p0, int _p1, int _p2, int same, SIMPLEPID pid, double min, double max, int hid, int _p3 ){
std::vector<int> a;
return combiner3(_p0,_p1,_p2, same,pid,min,max,a,_p3);
}
 
int Blab2::combiner(int _p0, int _p1,int same, SIMPLEPID pid, double min, double max, std::vector<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();
 std::vector<BParticle *> pall(3);
 int nprt=0;
 
 for(TIter next1(fList[_p0]);BParticle * p1 =(BParticle *) next1();) {
    // the second loop
   // in the case the second parti 
   for(TIter next2 = (_p0!=_p1 && same==0) ?  TIter(fList[_p1]): TIter(next1) ; BParticle * p2 =(BParticle *) next2();) {  
      if (p1==p2) continue;     // do not use the same particle in the combinations
      BParticle  p = *p1 + *p2; // Combine two particles into a new particle
      pall[0] = &p;
      pall[1] = p1;
      pall[2] = p2;
 
      if (p.InMassRange(min, max)){
            Fill(hid, pall);
            p.SetPid(pid); // set PID to particlename to fix the particle mass
            p.SetEnergyFromPid();
            TClonesArray& list = *fList[_p2];          
            new (list[nprt++]) BParticle ( p ); // create a new entry in kslist list of particles
           
      }
       
   }
               
 }
 return _p2;
}
 
 
int Blab2::combiner3(int _p0, int _p1,int _p2, int same, SIMPLEPID pid, double min, double max, std::vector<int> hid, int _p3 ){
   // Loop over all the particles in both lists.
 if (_p0 < 0 ) _p0 =0;
 if (_p1 < 0 ) _p1 =0;
 if (_p2 < 0 ) _p2 =0;
 
 
 fList[_p3]->Clear();
 std::vector<BParticle *> pall(4);  
 
 int nprt=0;
 
 for(TIter next1(fList[_p0]);BParticle * p1 =(BParticle *) next1();) {
    // the second loop
   // in the case the second parti 
   for(TIter next2 = (_p0!=_p1 && same==0) ?  TIter(fList[_p1]): TIter(next1) ; BParticle * p2 =(BParticle *) next2();) {  
      if (p1==p2) continue;     // do not use the same particle in the combinations
      for(TIter next3 = (_p1!=_p2 && same==0) ?  TIter(fList[_p2]): TIter(next2) ; BParticle * p3 =(BParticle *) next3();) {  
        if (p2==p3) continue;     // do not use the same particle in the combinations
        BParticle  p =  *p1 + *p2 + *p3 ; // Combine two particles into a new particle  
        pall[0] = &p;
        pall[1] = p1;
        pall[2] = p2;
        pall[3] = p3;
 
           
        if (p.InMassRange(min, max)){
            Fill(hid, pall);
            p.SetPid(pid); // set PID to particlename to fix the particle mass
            p.SetEnergyFromPid();
            TClonesArray& list = *fList[_p3];          
            new (list[nprt++]) BParticle ( p ); // create a new entry in kslist list of particles       
        }
      }
       
   }
               
 }
 return _p3;
}
 
 
int Blab2::selector(int pin, int charge, SIMPLEPID type ,  int  hid, int pout ){
std::vector<int> a;
return selector(pin,charge,type,a,pout);
}
int Blab2::selector(int pin, int charge, SIMPLEPID type ,  std::vector<int> hid, int pout ){
 if (pin < 0 ) pin =0;
 
  fList[pout]->Clear();
  int nprt=0;
  std::vector<BParticle *> pall(1);
 
  for(TIter next(fList[pin]); BParticle * p =(BParticle *) next();) {
        if (p->charge()== charge || charge > 1){
          if ( p->pid()== type || type == ALL ) {
            TClonesArray& list = *fList[pout];
            new (list[nprt++]) BParticle ( *p );
            pall[0] = p ;  
            Fill(hid, pall);
          }
        }
  }      
   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(){
 
char sList[0xFFFF];
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);
 
 
 
send_message(0, TString::Format("<br>Number of Events in the file %lld<br>\n", t->GetEntries() ),0);
TStopwatch timer;
timer.Start();
int nev  = 0;
int i    =TMath::Min(fNfirst, (UInt_t) t-> GetEntries());
int cNeve=TMath::Min(fNfirst+fNeve, (UInt_t) t-> GetEntries());
int fPart = fPrint;
int totaltracks = 0;
while (i<cNeve){
 t-> GetEntry(i); // Read the content of the event from the file
 fList[0]= mevent->GetParticleList();
 
 event();
 
 int progress = (100*i)/cNeve;
 if (i%10000==0) send_message(2, TString::Format("Event %d\n",i), progress);
 
 int nprt=0;
 if (nev>100) fPrint = 0; // disable particle prints for huge numer of events
 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);
 for(TIter next(fList[0]); BParticle * p =(BParticle *) next();) {
   nprt++;
   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()] );
 }
 if (fPart) fHnprt->Fill(nprt);
 totaltracks += nprt;
 if (fPrint) {
   sprintf(sList,"%s</table>",sList);
   send_message(0, TString(sList),progress);
   nev++;
 }
 mevent-> Clear();  // Clear the memory.
 for (int k=0;k<100;k++) if (fList[k]!=0) fList[k]->Clear();
 i++;
}
double avgtracks=(i)?float(totaltracks)/i:0;
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);
 
if (fPart) send_message(1,TBufferJSON::ConvertToJSON(fHnprt),100 );
 
 
for (int i=0;i<100;i++) if (fH[i]!=0) send_message(1,TBufferJSON::ConvertToJSON(fH[i]),100 );
 
TDatime d;
timer.Stop();
send_message(3, TString::Format("'%s', %d, %f, %f", d.AsSQLString(),i, timer.RealTime(), timer.CpuTime() ),100);
 
}
 
 
 

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

Subversion Repositories f9daq

(root)/belle2/masterclass/src/Blab2.cc @ 344 – Rev 331 → 333