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

Subversion Repositories f9daq

(root)/belle2/masterclass/src/Blab2.cc – Rev 277 → 319