FairRoot/PandaRoot
Functions
analysis_digi_cluster_7gammaAnalysis_fwendcap.C File Reference
#include <stdio.h>
#include <stdlib.h>
#include <vector>
#include "TF1.h"

Go to the source code of this file.

Functions

int analysis_digi_cluster_fwendcap ()
 

Function Documentation

int analysis_digi_cluster_fwendcap ( )

highestenergy_sumOf7clusters += highestenergy[k]; if (k==6) //eventually write the summation of the 7 highest cluster enegies

cout <<"highestenergy_sumOf7clusters = "<< highestenergy_sumOf7clusters <<'
'<<endl;

Definition at line 7 of file analysis_digi_cluster_7gammaAnalysis_fwendcap.C.

References acos(), bump_array, c10, cluster_array, cluster_energy, crystal_did, ctime, digi_array, digi_ene, digi_phi, digi_theta, Double_t, ene_mc, fsim, PndEmcCluster::GetEnergyCorrected(), PndEmcCluster::GetModule(), PndMCTrack::GetMomentum(), h1, h10, h10d, h2, h20, h20d, hene_mc, i, id_d, mctrack, mctrack_array, module_d, ncounts, ndigi, PndEmcCluster::NumberOfDigis(), outfile, p4mom, phi_mc, photon_momentum, Pi, printf(), row_did, rtime, sqrt(), theta_mc, timer, tsim, TString, PndEmcCluster::where(), PndEmcCluster::x(), PndEmcCluster::y(), and PndEmcCluster::z().

8 {
9  TStopwatch timer;
10  timer.Start();
11 
12  gROOT->LoadMacro("$VMCWORKDIR/macro/mvd/Tools.C");
13  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
14 
15  TFile* fsim = new TFile("/home/moini/myPandaRoot/myCodes/ROOTfiles/h_c_channel/emc_complete_onlyFwEndCap.root");
16 
17  TTree *tsim=(TTree *) fsim->Get("pndsim") ;
18  TClonesArray* mctrack_array=new TClonesArray("PndMCTrack");
19  tsim->SetBranchAddress("MCTrack",&mctrack_array);
20 
21  TClonesArray* digi_array=new TClonesArray("PndEmcDigi");
22  tsim->SetBranchAddress("EmcDigi",&digi_array);
23 
24  TClonesArray* cluster_array=new TClonesArray("PndEmcCluster");
25  tsim->SetBranchAddress("EmcCluster",&cluster_array);
26 
27  TClonesArray* bump_array=new TClonesArray("PndEmcBump");
28  tsim->SetBranchAddress("EmcBump",&bump_array);
29 
30  Double_t theta_mc, phi_mc, ene_mc,
31  digi_ene, digi_theta, digi_phi;
32  TVector3 photon_momentum;
33  TLorentzVector p4mom;
34  Double_t clusterThresholdEnergy=0;
35 
36  TH1F *hene_mc= new TH1F("hene_mc","MC energy (GeV)",100,0.0,1.05);
37  TH1F *h10= new TH1F("h10"," Theta (MC)",180,0.,180.);
38  TH1F *h20= new TH1F("h20"," Phi (MC)",360,0.,360.);
39  TH1F *h10d= new TH1F("h10d"," Theta (data)",180,0.,180.);
40  TH1F *h20d= new TH1F("h20d"," Phi (data)",360,0.,360.);
41  TH1F *h1= new TH1F("h1"," Theta difference",100,-2.,2.);
42  TH1F *h2= new TH1F("h2"," Phi difference",100,-2.,2.);
43 
44  Int_t ncounts = tsim->GetEntriesFast();
45  cout << "Number of events = " << ncounts<< endl;
46  //ncounts=100; //in case you want to analyse specific number of events
47 
48  Int_t module_d, id_d, crystal_did, row_did;
49  TH1F *HISTdigi_ene= new TH1F("HISTdigi_ene","FwEndCap energy (digi)",2000,0.0,2.05); // deposited energy per event in the whole FwEndCap
50  TH1F *HISTene_d= new TH1F("HISTene_d","energy (digi)",1000,0.0001,1.01); // deposited energy per digi in the whole FwEndCap
51  TH2F *HISTxpad_x_d = new TH2F("HISTxpad_x_d","xpad vs x (digi)",200,200,300,2*38.5,-38.5,38.5); // here x means: map of the x-position of the hit crystal
52  TH2F *HISTypad_y_d = new TH2F("HISTypad_y_d","ypad vs y (digi)",200,200,300,2*38.5,-38.5,38.5); // here y means: map of the y-position of the hit crystal
53  TH2F *HISTrow_x_d = new TH2F("HISTrow_x_d","row vs x (digi)",200,200,300,2*38.5,-38.5,38.5); // here x means: map of the x-position of the hit crystal
54  TH2F *HISTrow_crystal_d = new TH2F("HISTrow_crystal_d","row vs column (digi)",2*38.5,-38.5,38.5,2*38.5,-38.5,38.5);
55  TH2F *HISTy_x_d = new TH2F("HISTy_x_d","y vs x (digi map)",200,200,300,200,200,300); // here x and y mean: map of the y- and x-position of the hit crystal
56  TH1F *HISTOdigi_ene11_12= new TH1F("HISTOdigi_ene11_12","crystal [11,12] dep. energy (GeV)",100,0.0001,1.01);
57  TH1F *HISTOdigi_ene12_12= new TH1F("HISTOdigi_ene12_12","crystal [12,12] dep. energy (GeV)",100,0.0001,1.01);
58  TH1F *HISTOdigi_ene13_12= new TH1F("HISTOdigi_ene13_12","crystal [13,12] dep. energy (GeV)",100,0.0001,1.01);
59  TH1F *HISTOdigi_ene14_12= new TH1F("HISTOdigi_ene14_12","crystal [14,12] dep. energy (GeV)",100,0.0001,1.01);
60  TH1F *HISTtotdigis_fromClusterAnalysis = new TH1F("HISTtotdigis_fromClusterAnalysis","total number of digis per event",100,0,25);
61  TH1I *HISTmultiplicity_5MeVtrigThresh_d = new TH1I("HISTmultiplicity_5MeVtrigThresh_d","Crystal digi multiplicity per event (3 MeV detection threshold and considering 5 MeV trigger threshold)",40,0.5,40.5);
62  TH1I *HISTmultiplicity_10MeVtrigThresh_d = new TH1I("HISTmultiplicity_10MeVtrigThresh_d","Crystal digi multiplicity per event (3 MeV detection threshold and considering 10 MeV trigger threshold)",40,0.5,40.5);
63 
64  TH1I *HISTnumberOfClusters = new TH1I("HISTnumberOfClusters","Total number of clusters per event, having more than or equal to one digi per cluster",100,0,20);
65  TH1I *HISTclustNo = new TH1I("HISTclustNo","Total number of clusters per event, having more than one digi per cluster and satisfying the cluster threshold energy",1000,0,20);
66  TH1I *HISTndigisOfHighestEnergyCluster = new TH1I("HISTndigisOfHighestEnergyCluster","Number of digis per event for clusters with highest energy deposition",1000,0,30);
67  TH1F *HISTmax_cluster_energy= new TH1F("HISTOmax_cluster_energy","Maximum cluster energy deposition per event",500,0.,1.05);
68  TH1F *HISTmax_cluster_theta = new TH1F("HISTmax_cluster_theta","#theta per event of the cluster with maximum deposition energy",100,0,30);
69  TH1F *HISTmax_cluster_phi = new TH1F("HISTmax_cluster_phi","#phi per event of the cluster with maximum deposition energy",370,-185,185);
70  TH1F *HISTmax_cluster_x = new TH1F("HISTmax_cluster_x","x per event of the cluster with maximum deposition energy",220,-110,110);
71  TH1F *HISTmax_cluster_y = new TH1F("HISTmax_cluster_y","y per event of the cluster with maximum deposition energy",220,-110,110);
72  TH1F *HISTmax_cluster_z = new TH1F("HISTmax_cluster_z","z per event of the cluster with maximum deposition energy",80,208,216);
73  TH2F *HISTmax_cluster_y_x = new TH2F("HISTmax_cluster_y_x","y vs x per event of the cluster with maximum deposition energy",220,-110,110,220,-110,110);
74  TH1F *HISTtheta_MonteCarlo = new TH1F("HISTtheta_MonteCarlo","Monte Carlo #theta per event of the cluster with maximum deposition energy",100,0,30);
75  TH1F *HISTphi_MonteCarlo = new TH1F("HISTphi_MonteCarlo","Monte Carlo #phi per event of the cluster with maximum deposition energy",370,-185,185);
76  TH1F *HISTmax_cluster_theta_diff = new TH1F("HISTmax_cluster_theta_diff","Difference between Monte Carlo and reconstructed #theta per event of the cluster with maximum deposition energy",100,-5,5);
77  TH1F *HISTmax_cluster_phi_diff = new TH1F("HISTmax_cluster_phi_diff","Difference between Monte Carlo and reconstructed #phi per event of the cluster with maximum deposition energy",360,-5,5);
78  TH2F *HISTtheta_MonteCarlo_cluster = new TH2F("HISTtheta_MonteCarlo_cluster","Monte Carlo versus cluster #theta of the cluster with maximum deposition energy",100,0,30,100,0,30);
79  TH2F *HISTphi_MonteCarlo_cluster = new TH2F("HISTphi_MonteCarlo_cluster","Monte Carlo versus cluster #phi of the cluster with maximum deposition energy",370,-185,185,370,-185,185);
80  TH1F *HISTmax_cluster_x_diff = new TH1F("HISTmax_cluster_x_diff","Difference between Monte Carlo and reconstructed x per event of the cluster with maximum deposition energy",100,-5,5);
81  TH1F *HISTmax_cluster_y_diff = new TH1F("HISTmax_cluster_y_diff","Difference between Monte Carlo and reconstructed y per event of the cluster with maximum deposition energy",100,-5,5);
82  TH2F *HISTclusterSize_clusterPhi = new TH2F("HISTclusterSize_clusterPhi","number of digis vs cluster #phi for clusters with highest energy deposition",370,-185,185,100,0,25);
83  TH2F *HISTclusterEnergy_clusterPhi = new TH2F("HISTclusterEnergy_clusterPhi","cluster energy deposition vs cluster #phi for clusters with highest energy deposition",370,-185,185,500,0.,1.05);
84  TH2F *HISTnumberOfclusters_clusterPhi = new TH2F("HISTnumberOfclusters_clusterPhi","number of clusters vs cluster #phi for clusters with highest energy deposition",370,-185,185,100,0,20);
85  TH2F *HISTtotdigis_clusterPhi = new TH2F("HISTtotdigis_clusterPhi","total number of digis vs cluster #phi (#phi of clusters with highest energy deposition)",370,-185,185,100,0,25);
86 
87  TH1F *HISTinvariantMass = new TH1F("HISTinvariantMass","invariant mass of every two clusters (7 highest energy clusters in total)",3200,0,3.2);
88  TH2F *HISTinvariantMass_openingAngle = new TH2F("HISTinvariantMass_openingAngle","opening angle of the two clusters versus the invariant mass of the assumed two gammas that make the two clusters",3200,0,3.2,140,0,50);
89  TH1F *HIST_Egamma1Egamma2 = new TH1F("HIST_Egamma1Egamma2","E_{#gamma1} #times E_{#gamma2}",8000,0,8.);
90  TH2F *HIST_Egamma1vsEgamma2 = new TH2F("HIST_Egamma1vsEgamma2","E_{#gamma2} versus E_{#gamma1}",8000,0,8.,3200,0,3.2);
91  TH2F *HIST_Egamma1Egamma2_openingAngle = new TH2F("HIST_Egamma1Egamma2_openingAngle","opening angle of the two clusters versus E_{#gamma1} #times E_{#gamma2}",8000,0,8.,140,0,50);
92  TH2F *HISTinvariantMass_Egamma1Egamma2 = new TH2F("HISTinvariantMass_Egamma1Egamma2","E_{#gamma1} #times E_{#gamma2} versus the invariant mass of the two gammas",3200,0,3.2,8000,0,8.);
93  TH2F *HIST_Egamma1_Egamma1Egamma2 = new TH2F("HIST_Egamma1_Egamma1Egamma2","E_{#gamma1} #times E_{#gamma2} versus E_{#gamma1}",6400,0,6.4,8000,0,8.);
94  TH2F *HIST_Egamma1_openingAngle = new TH2F("HIST_Egamma1_openingAngle","opening angle of the two clusters versus E_{#gamma1}",6400,0,6.4,140,0,50);
95  TH3F *HISTinvariantMass_openingAngle_Egamma1Egamma2 = new TH3F("HISTinvariantMass_openingAngle_Egamma1Egamma2","",150,0,1.5,140,0,50,400,0,4.);
96 
97 
98  Int_t totaldigis=0; //counter on the total number of digis (fired detectors) in all events
99  Int_t efficiencyCounter_d=0;
100  Int_t atLeastOneHitPerEventCounter_d=0;
101 
102  for (int j = 0; j < ncounts; j++) {
103 
104  Int_t efficiencyFlag10MeV_d=0;
105  Int_t efficiencyFlag5MeV_d=0;
106  Int_t atLeastOneHitPerEventFlag_d=0;
107  Double_t sum_digi_ene=0;
108 
109  tsim->GetEntry(j);
110 
111  PndMCTrack *mctrack=(PndMCTrack *) mctrack_array->At(0);
112  photon_momentum=mctrack->GetMomentum();
113  theta_mc=photon_momentum.Theta()*(180./TMath::Pi());
114  phi_mc=photon_momentum.Phi()*(180./TMath::Pi());
115  //p4mom=mctrack->Get4Momentum();
116 /*
118  Int_t ndigis = digi_array->GetEntries();
119  totaldigis+=ndigis;
120  if (ndigis!=0)
121  atLeastOneHitPerEventFlag_d = 1;
122 
123  for (Int_t i=0; i<ndigis; i++)
124  {
125  PndEmcDigi *digi=(PndEmcDigi*)digi_array->At(i);
126  TVector3 digi_pos=digi->where();
127  module_d = digi->GetModule();
128  if (module_d==3){
129  //dene_mc->Fill(p4mom.E());
130 
131  digi_ene=digi->GetEnergy();
132  HISTene_d->Fill(digi_ene);
133  if (digi_ene >= 0.010) // 10 MeV trigger threshold used (here only) to calculate the efficiency of the FwEndCap geometry
134  efficiencyFlag10MeV_d = 1;
135  if (digi_ene >= 0.005) // 5 MeV trigger threshold used (here only) to calculate the efficiency of the FwEndCap geometry
136  efficiencyFlag5MeV_d = 1;
137 
138  //filling some histograms for a few crystals
139  if (digi->GetXPad()==11 && digi->GetYPad()==12)
140  HISTOdigi_ene11_12->Fill(digi_ene);
141  else if (digi->GetXPad()==12 && digi->GetYPad()==12)
142  HISTOdigi_ene12_12->Fill(digi_ene);
143  else if (digi->GetXPad()==13 && digi->GetYPad()==12)
144  HISTOdigi_ene13_12->Fill(digi_ene);
145  else if (digi->GetXPad()==14 && digi->GetYPad()==12)
146  HISTOdigi_ene14_12->Fill(digi_ene);
147 
148  if (digi_ene >= 0.010)
149  sum_digi_ene+=digi_ene;
150 
151  //Double_t digi_z=digi_pos.Z();
152  //hz_d->Fill(digi_z);
153 
154  digi_theta=(digi->GetTheta())*(180./TMath::Pi());
155  digi_phi=(digi->GetPhi())*(180./TMath::Pi());
156  //h10digi->Fill(theta_mc);
157  //h20digi->Fill(phi_mc);
158  //h10ddigi->Fill(digi_theta);
159  //h20ddigi->Fill(digi_phi);
160  //h1d->Fill(digi_theta-theta_mc);
161  //h2d->Fill(digi_phi-phi_mc);
162 
163  id_d = digi->GetDetectorId();
164  crystal_did = -(id_d%10000 - 36); //crystal column number (the same as GetXPad())
165  row_did = (id_d/1000000)%100 - 37; //crystal row number (the same as GetYPad())
166 
167  HISTxpad_x_d->Fill(digi->GetThetaInt(),digi->GetXPad());
168  HISTypad_y_d->Fill(digi->GetPhiInt(),digi->GetYPad());
169  HISTrow_x_d->Fill(digi->GetThetaInt(),digi->GetYPad());
170 
171  //hid_d->Fill(id_d);
172  //hmodule_d->Fill(module_d);
173  HISTrow_crystal_d->Fill(digi->GetXPad(),digi->GetYPad());
174  HISTy_x_d->Fill(digi->GetThetaInt(),digi->GetPhiInt());
175  }
176  }
177 
178  if (efficiencyFlag10MeV_d==1) //if the trigger threshold is satisfied
179  HISTdigi_ene->Fill(sum_digi_ene); // deposited energy per event in all digis (whole FwEndCap)
180 
181  if (efficiencyFlag10MeV_d==1)
182  {
183  efficiencyCounter_d+=1;
184  HISTmultiplicity_10MeVtrigThresh_d->Fill(ndigis); //10 MEV trigger threshold and 3 MeV detection threshold (in the level of digi this 3 MeV is already implemented)
185  }
186  if (efficiencyFlag5MeV_d==1)
187  HISTmultiplicity_5MeVtrigThresh_d->Fill(ndigis); //5 MEV trigger threshold and 3 MeV detection threshold (in the level of digi this 3 MeV is already implemented)
188 
189  if (atLeastOneHitPerEventFlag_d==1)
190  atLeastOneHitPerEventCounter_d+=1;
191 */
192 
194  Int_t NumberOfClusters, ClustNo, totdigis_fromClusterAnalysis, ndigisOfHighestEnergyCluster;
195  NumberOfClusters=ClustNo=totdigis_fromClusterAnalysis=ndigisOfHighestEnergyCluster=0; //set per event
196  Int_t idClusterSatisfiedThresholdEnergy=-1;
197  Double_t highestEnergy=0; // (set per event)
198  Double_t highestenergy_sumOf7clusters=0; // (set per event)
199 
200  Double_t highestenergy[7], cluster_x_7[7], cluster_y_7[7],cluster_z_7[7],cluster_dist_7[7];
201  for (int k=0; k<7;k++){
202  highestenergy[k]=0; // (set per event)
203  cluster_x_7[k]=cluster_y_7[k]=cluster_z_7[k]=cluster_dist_7[k]=-1000; // (set per event)
204  }
205 
206  Int_t originalClusterIndex[7];
207  for (int k=0; k<7;k++)
208  originalClusterIndex[k]=-1; //stores the original indices of the 7 clusters with biggest energies (set per event)
209 
210  if (cluster_array->GetEntriesFast() > 0)
211  {
212  NumberOfClusters=cluster_array->GetEntriesFast();
213  std::vector <double> clusterEnergies;
214  std::vector <double> clusterEnergies_original;
215 
216  for (Int_t i=0; i<NumberOfClusters; i++) //loop over all clusters of the present event under analysis (ndigi>=1)
217  {
218  PndEmcCluster *cluster=(PndEmcCluster*)cluster_array->At(i);
219  if (cluster->GetModule()==3)
220  {
221  Double_t cluster_energy=cluster->GetEnergyCorrected();//energy();
222  clusterEnergies.push_back(cluster_energy); //first push all the cluster energies into the vector; no sorting yet
223  clusterEnergies_original.push_back(cluster_energy); //make another one, in order to keep track of indices
224 
225  Int_t ndigi=cluster->NumberOfDigis(); // number of fired detectors inside the cluster
226  totdigis_fromClusterAnalysis += ndigi;
227  if (ndigi >= 1 && cluster_energy >= clusterThresholdEnergy) //by this definition for a cluster, ndigi should be more than 1
228  ClustNo++;
229 
230  /*if (cluster_energy > highestEnergy) //finding the cluster with highest deposited energy
231  {
232  idClusterSatisfiedThresholdEnergy = i;
233  highestEnergy = cluster_energy;
234  ndigisOfHighestEnergyCluster = ndigi;
235  }*/
236  }
237  }
238 
239  /*
240  // Per-event analysis of the cluster with highest deposited energy
241  if (idClusterSatisfiedThresholdEnergy != -1) {
242  PndEmcCluster *cluster=(PndEmcCluster*)cluster_array->At(idClusterSatisfiedThresholdEnergy);
243  TVector3 cluster_pos=cluster->where();
244  Double_t cluster_theta=cluster->theta()*(180./TMath::Pi()); //cluster_theta=cluster_pos.Theta()*(180./TMath::Pi());
245  Double_t cluster_phi=cluster->phi()*(180./TMath::Pi()); //cluster_phi=cluster_pos.Phi()*(180./TMath::Pi());
246  Double_t cluster_x=cluster->x();
247  Double_t cluster_y=cluster->y();
248  Double_t cluster_z=cluster->z();
249  Double_t cluster_energy=cluster->GetEnergyCorrected();//energy();
250  //sum_clus_ene+=cluster_energy; //sum up the highest deposited cluster energy of all events
251  HISTmax_cluster_energy->Fill(cluster_energy); //fill the highest deposited cluster energy per event
252  HISTmax_cluster_theta->Fill(cluster_theta); // reconstructed
253  HISTmax_cluster_theta_diff->Fill(cluster_theta-theta_mc);
254  //if (cluster_theta>10)
255  // {
256  HISTmax_cluster_phi->Fill(cluster_phi); // reconstructed
257  HISTmax_cluster_phi_diff->Fill(cluster_phi-phi_mc);
258  HISTmax_cluster_x->Fill(cluster_x); // reconstructed
259  HISTmax_cluster_y->Fill(cluster_y); // reconstructed
260  HISTmax_cluster_y_x->Fill(cluster_x,cluster_y);
261  HISTmax_cluster_z->Fill(cluster_z); // reconstructed
262  Double_t x_mc = cluster_z*TMath::Tan(theta_mc/180.*TMath::Pi())*TMath::Cos(phi_mc/180.*TMath::Pi());
263  Double_t y_mc = cluster_z*TMath::Tan(theta_mc/180.*TMath::Pi())*TMath::Sin(phi_mc/180.*TMath::Pi());
264  HISTmax_cluster_x_diff->Fill(cluster_x-x_mc);
265  HISTmax_cluster_y_diff->Fill(cluster_y-y_mc);
266  //}
267  HISTtheta_MonteCarlo->Fill(theta_mc);
268  HISTphi_MonteCarlo->Fill(phi_mc);
269  HISTtheta_MonteCarlo_cluster->Fill(cluster_theta,theta_mc);
270  HISTphi_MonteCarlo_cluster->Fill(cluster_phi,phi_mc);
271  HISTclusterSize_clusterPhi->Fill(cluster_phi,ndigisOfHighestEnergyCluster);
272  HISTclusterEnergy_clusterPhi->Fill(cluster_phi,cluster_energy);
273  HISTnumberOfclusters_clusterPhi->Fill(cluster_phi,NumberOfClusters);
274  HISTtotdigis_clusterPhi->Fill(cluster_phi,totdigis_fromClusterAnalysis);
275  }*/
276 
277  // Per-event analysis of the 7 clusters with highest deposited energy
278  std::sort(clusterEnergies.begin(), clusterEnergies.end()); //now sort it
279 
280  if (NumberOfClusters>=7) //finding at least 7 clusters
281  {
282  for (int k=0; k<7; k++) //finding the first 7 biggest cluster energies (k=0 corresponds to the cluster with biggest energy)
283  {
284  highestenergy[k] = clusterEnergies[(NumberOfClusters-1)-k];
285  //cout << "decsendingly sorted energy of cluster # " << k+1 << " = " << highestenergy[k] << '\n';
289 
290  for (int i=0; i<NumberOfClusters; i++)//finding the original cluster indices of the 7 biggest energy clusters
291  if (!(clusterEnergies_original[i]>highestenergy[k]) && !(clusterEnergies_original[i]<highestenergy[k]))//practically the equality sign
292  originalClusterIndex[k]=i;
293 
294  PndEmcCluster *cluster=(PndEmcCluster*)cluster_array->At(originalClusterIndex[k]);
295  TVector3 cluster_pos=cluster->where();
296  cluster_x_7[k]=cluster->x();
297  cluster_y_7[k]=cluster->y();
298  cluster_z_7[k]=cluster->z();
299  cluster_dist_7[k]=cluster_pos.Mag();//distance of the center of the cluster to the target point
300 
301  }
302 
303  for (int k0=1; k0<7; k0++)
304  {
305  Double_t cosTheta = (1./cluster_dist_7[0]/cluster_dist_7[k0])*(cluster_x_7[0]*cluster_x_7[k0]+cluster_y_7[0]*cluster_y_7[k0]+cluster_z_7[0]*cluster_z_7[k0]);
306  Double_t invMass = sqrt(2*highestenergy[0]*highestenergy[k0]*(1-cosTheta));
307  HISTinvariantMass->Fill(invMass);
308  HISTinvariantMass_openingAngle->Fill(invMass,acos(cosTheta)*180/3.1415926535);
309  HIST_Egamma1Egamma2->Fill(highestenergy[0]*highestenergy[k0]);// (E_gamma1)(E_gamma2)
310  HIST_Egamma1vsEgamma2->Fill(highestenergy[0],highestenergy[k0]);// E_gamma1 versus E_gamma2
311  HIST_Egamma1Egamma2_openingAngle->Fill(highestenergy[0]*highestenergy[k0],acos(cosTheta)*180/3.1415926535);
312  HISTinvariantMass_Egamma1Egamma2->Fill(invMass,highestenergy[0]*highestenergy[k0]);
313  HIST_Egamma1_Egamma1Egamma2->Fill(highestenergy[0],highestenergy[0]*highestenergy[k0]);
314  HIST_Egamma1_openingAngle->Fill(highestenergy[0],acos(cosTheta)*180/3.1415926535);
315  HISTinvariantMass_openingAngle_Egamma1Egamma2->Fill(invMass,acos(cosTheta)*180/3.1415926535,highestenergy[0]*highestenergy[k0]);
316  }
317  for (int k1=2; k1<7; k1++)
318  {
319  Double_t cosTheta = (1./cluster_dist_7[1]/cluster_dist_7[k1])*(cluster_x_7[1]*cluster_x_7[k1]+cluster_y_7[1]*cluster_y_7[k1]+cluster_z_7[1]*cluster_z_7[k1]);
320  Double_t invMass = sqrt(2*highestenergy[1]*highestenergy[k1]*(1-cosTheta));
321  HISTinvariantMass->Fill(invMass);
322  HISTinvariantMass_openingAngle->Fill(invMass,acos(cosTheta)*180/3.1415926535);
323  HIST_Egamma1Egamma2->Fill(highestenergy[1]*highestenergy[k1]);// (E_gamma1)(E_gamma2)
324  HIST_Egamma1vsEgamma2->Fill(highestenergy[1],highestenergy[k1]);// E_gamma1 versus E_gamma2
325  HIST_Egamma1Egamma2_openingAngle->Fill(highestenergy[1]*highestenergy[k1],acos(cosTheta)*180/3.1415926535);
326  HISTinvariantMass_Egamma1Egamma2->Fill(invMass,highestenergy[1]*highestenergy[k1]);
327  HIST_Egamma1_Egamma1Egamma2->Fill(highestenergy[1],highestenergy[1]*highestenergy[k1]);
328  HIST_Egamma1_openingAngle->Fill(highestenergy[1],acos(cosTheta)*180/3.1415926535);
329  HISTinvariantMass_openingAngle_Egamma1Egamma2->Fill(invMass,acos(cosTheta)*180/3.1415926535,highestenergy[1]*highestenergy[k1]);
330  }
331  for (int k2=3; k2<7; k2++)
332  {
333  Double_t cosTheta = (1./cluster_dist_7[2]/cluster_dist_7[k2])*(cluster_x_7[2]*cluster_x_7[k2]+cluster_y_7[2]*cluster_y_7[k2]+cluster_z_7[2]*cluster_z_7[k2]);
334  Double_t invMass = sqrt(2*highestenergy[2]*highestenergy[k2]*(1-cosTheta));
335  HISTinvariantMass->Fill(invMass);
336  HISTinvariantMass_openingAngle->Fill(invMass,acos(cosTheta)*180/3.1415926535);
337  HIST_Egamma1Egamma2->Fill(highestenergy[2]*highestenergy[k2]);// (E_gamma1)(E_gamma2)
338  HIST_Egamma1vsEgamma2->Fill(highestenergy[2],highestenergy[k2]);// E_gamma1 versus E_gamma2
339  HIST_Egamma1Egamma2_openingAngle->Fill(highestenergy[2]*highestenergy[k2],acos(cosTheta)*180/3.1415926535);
340  HISTinvariantMass_Egamma1Egamma2->Fill(invMass,highestenergy[2]*highestenergy[k2]);
341  HIST_Egamma1_Egamma1Egamma2->Fill(highestenergy[2],highestenergy[2]*highestenergy[k2]);
342  HIST_Egamma1_openingAngle->Fill(highestenergy[2],acos(cosTheta)*180/3.1415926535);
343  HISTinvariantMass_openingAngle_Egamma1Egamma2->Fill(invMass,acos(cosTheta)*180/3.1415926535,highestenergy[2]*highestenergy[k2]);
344  }
345  for (int k3=4; k3<7; k3++)
346  {
347  Double_t cosTheta = (1./cluster_dist_7[3]/cluster_dist_7[k3])*(cluster_x_7[3]*cluster_x_7[k3]+cluster_y_7[3]*cluster_y_7[k3]+cluster_z_7[3]*cluster_z_7[k3]);
348  Double_t invMass = sqrt(2*highestenergy[3]*highestenergy[k3]*(1-cosTheta));
349  HISTinvariantMass->Fill(invMass);
350  HISTinvariantMass_openingAngle->Fill(invMass,acos(cosTheta)*180/3.1415926535);
351  HIST_Egamma1Egamma2->Fill(highestenergy[3]*highestenergy[k3]);// (E_gamma1)(E_gamma2)
352  HIST_Egamma1vsEgamma2->Fill(highestenergy[3],highestenergy[k3]);// E_gamma1 versus E_gamma2
353  HIST_Egamma1Egamma2_openingAngle->Fill(highestenergy[3]*highestenergy[k3],acos(cosTheta)*180/3.1415926535);
354  HISTinvariantMass_Egamma1Egamma2->Fill(invMass,highestenergy[3]*highestenergy[k3]);
355  HIST_Egamma1_Egamma1Egamma2->Fill(highestenergy[3],highestenergy[3]*highestenergy[k3]);
356  HIST_Egamma1_openingAngle->Fill(highestenergy[3],acos(cosTheta)*180/3.1415926535);
357  HISTinvariantMass_openingAngle_Egamma1Egamma2->Fill(invMass,acos(cosTheta)*180/3.1415926535,highestenergy[3]*highestenergy[k3]);
358  }
359  for (int k4=5; k4<7; k4++)
360  {
361  Double_t cosTheta = (1./cluster_dist_7[4]/cluster_dist_7[k4])*(cluster_x_7[4]*cluster_x_7[k4]+cluster_y_7[4]*cluster_y_7[k4]+cluster_z_7[4]*cluster_z_7[k4]);
362  Double_t invMass = sqrt(2*highestenergy[4]*highestenergy[k4]*(1-cosTheta));
363  HISTinvariantMass->Fill(invMass);
364  HISTinvariantMass_openingAngle->Fill(invMass,acos(cosTheta)*180/3.1415926535);
365  HIST_Egamma1Egamma2->Fill(highestenergy[4]*highestenergy[k4]);// (E_gamma1)(E_gamma2)
366  HIST_Egamma1vsEgamma2->Fill(highestenergy[4],highestenergy[k4]);// E_gamma1 versus E_gamma2
367  HIST_Egamma1Egamma2_openingAngle->Fill(highestenergy[4]*highestenergy[k4],acos(cosTheta)*180/3.1415926535);
368  HISTinvariantMass_Egamma1Egamma2->Fill(invMass,highestenergy[4]*highestenergy[k4]);
369  HIST_Egamma1_Egamma1Egamma2->Fill(highestenergy[4],highestenergy[4]*highestenergy[k4]);
370  HIST_Egamma1_openingAngle->Fill(highestenergy[4],acos(cosTheta)*180/3.1415926535);
371  HISTinvariantMass_openingAngle_Egamma1Egamma2->Fill(invMass,acos(cosTheta)*180/3.1415926535,highestenergy[4]*highestenergy[k4]);
372  }
373  for (int k5=6; k5<7; k5++)
374  {
375  Double_t cosTheta = (1./cluster_dist_7[5]/cluster_dist_7[k5])*(cluster_x_7[5]*cluster_x_7[k5]+cluster_y_7[5]*cluster_y_7[k5]+cluster_z_7[5]*cluster_z_7[k5]);
376  Double_t invMass = sqrt(2*highestenergy[5]*highestenergy[k5]*(1-cosTheta));
377  HISTinvariantMass->Fill(invMass);
378  HISTinvariantMass_openingAngle->Fill(invMass,acos(cosTheta)*180/3.1415926535);
379  HIST_Egamma1Egamma2->Fill(highestenergy[5]*highestenergy[k5]);// (E_gamma1)(E_gamma2)
380  HIST_Egamma1vsEgamma2->Fill(highestenergy[5],highestenergy[k5]);// E_gamma1 versus E_gamma2
381  HIST_Egamma1Egamma2_openingAngle->Fill(highestenergy[5]*highestenergy[k5],acos(cosTheta)*180/3.1415926535);
382  HISTinvariantMass_Egamma1Egamma2->Fill(invMass,highestenergy[5]*highestenergy[k5]);
383  HIST_Egamma1_Egamma1Egamma2->Fill(highestenergy[5],highestenergy[5]*highestenergy[k5]);
384  HIST_Egamma1_openingAngle->Fill(highestenergy[5],acos(cosTheta)*180/3.1415926535);
385  HISTinvariantMass_openingAngle_Egamma1Egamma2->Fill(invMass,acos(cosTheta)*180/3.1415926535,highestenergy[5]*highestenergy[k5]);
386  }
387  }
388  }
389 
390  //hz_clust->Fill(cluster->z());
391 
392  if (NumberOfClusters>0)
393  HISTnumberOfClusters->Fill(NumberOfClusters); //number of clusters per event, having more than or equal to one digi per cluster
394  if (ClustNo>0)
395  HISTclustNo->Fill(ClustNo); //number of clusters per event, having more than one digi per cluster and satisfying the cluster threshold energy
396  if (ndigisOfHighestEnergyCluster>0)
397  HISTndigisOfHighestEnergyCluster->Fill(ndigisOfHighestEnergyCluster); //number of digis per event for clusters with highest energy deposition
398 
399  HISTtotdigis_fromClusterAnalysis->Fill(totdigis_fromClusterAnalysis);
400 
401  if (j%1000 ==0)
402  cout<<"Evt. No. "<< j<<endl;
403  }
404  //cout<<"totaldigis = "<<totaldigis<<"\n";
405  //cout<<"efficiencyCounter_d = "<<efficiencyCounter_d<<" , atLeastOneHitPerEventCounter_d = "<<atLeastOneHitPerEventCounter_d<<"\n";
406  //cout<<"FwEndCap efficiency (more than 10 MeV deposition in at least one crystal) = "<< (double)efficiencyCounter_d/atLeastOneHitPerEventCounter_d<<"\n";
407 
408 
410  gStyle->SetPalette(1);
411 
412 
413 
414 /*
415  TCanvas* c3 = new TCanvas("c3","", 100, 100, 800, 800);
416  c3->Divide(3,2);
417  c3->cd(1);
418  HISTOdigi_ene11_12->Draw();
419  c3->cd(2);
420  HISTOdigi_ene12_12->Draw();
421  c3->cd(3);
422  HISTOdigi_ene13_12->Draw();
423  c3->cd(4);
424  HISTOdigi_ene14_12->Draw();
425  c3->cd(5);
426  HISTene_d->Draw();
427  HISTene_d->GetXaxis()->SetTitle("FwEndCap Deposited Energy per digi (digi analysis)");
428  c3->cd(6);
429  HISTdigi_ene->SetLineColor(1);
430  HISTdigi_ene->Draw();
431  HISTdigi_ene->GetXaxis()->SetTitle("FwEndCap Deposited Energy per event (digi analysis)");
432  */
433 
434  /*
435  TCanvas* c4 = new TCanvas("c4","", 100, 100, 800, 800);
436  c4->Divide(3,2);
437  c4->cd(1);
438  HISTy_x_d->Draw("colz");
439  c4->cd(2);
440  HISTrow_crystal_d->Draw("colz");
441  c4->cd(3);
442  HISTxpad_x_d->Draw("colz");
443  c4->cd(4);
444  HISTypad_y_d->Draw("colz");
445  c4->cd(5);
446  HISTrow_x_d->Draw("colz");
447  */
448 
449  /*
450  TCanvas* c5 = new TCanvas("c5","", 100, 100, 800, 800);
451  c5->Divide(4,4);
452  c5->cd(1);
453  HISTnumberOfClusters->Draw("");
454  c5->cd(2);
455  HISTclustNo->Draw("");
456  c5->cd(3);
457  HISTndigisOfHighestEnergyCluster->Draw("");
458  c5->cd(4);
459  HISTmax_cluster_energy->Draw("");
460  c5->cd(5);
461  HISTtheta_MonteCarlo->Draw("");
462  c5->cd(6);
463  HISTphi_MonteCarlo->Draw("");
464  c5->cd(7);
465  HISTtheta_MonteCarlo_cluster->Draw("colz");
466  c5->cd(8);
467  HISTphi_MonteCarlo_cluster->Draw("colz");
468  c5->cd(9);
469  HISTmax_cluster_theta->Draw("");
470  c5->cd(10);
471  HISTmax_cluster_phi->Draw("");
472  c5->cd(11);
473  HISTmax_cluster_theta_diff->Draw("");
474  c5->cd(12);
475  HISTmax_cluster_phi_diff->Draw("");
476  c5->cd(13);
477  HISTmax_cluster_x->Draw("");
478  c5->cd(14);
479  HISTmax_cluster_y->Draw("");
480  c5->cd(15);
481  //HISTmax_cluster_z->Draw("");
482  HISTmax_cluster_x_diff->Draw("");
483  c5->cd(16);
484  HISTmax_cluster_y_diff->Draw("");
485 
486  TCanvas* c55 = new TCanvas("c55","", 100, 100, 800, 800);
487  HISTmax_cluster_y_x->Draw("colz");
488  */
489 
490  /* TCanvas* c6 = new TCanvas("c6","", 100, 100, 800, 800);
491  c6->Divide(2,2);
492  c6->cd(1);
493  HISTclusterSize_clusterPhi->Draw("colz");
494  c6->cd(2);
495  HISTclusterEnergy_clusterPhi->Draw("colz");
496  c6->cd(3);
497  HISTnumberOfclusters_clusterPhi->Draw("colz");
498  c6->cd(4);
499  HISTtotdigis_clusterPhi->Draw("colz");
500  //HISTtotdigis_fromClusterAnalysis->Draw();
501  */
502 
503  /*
504  TCanvas* c8 = new TCanvas("c8","", 100, 100, 1200, 800);
505  c8->Divide(2,2);
506  c8->cd(1);
507  HISTmultiplicity_5MeVtrigThresh_d->Draw();
508  c8->cd(2);
509  HISTmultiplicity_10MeVtrigThresh_d->Draw();
510 */
511 
512 
513  TCanvas* c10 = new TCanvas("c10","", 100, 100, 1200, 800);
514  c10->Divide(2,4);
515  c10->cd(1)->SetLogy();
516  HISTinvariantMass->Draw();
517  c10->cd(2)->SetLogy();
518  HIST_Egamma1Egamma2->Draw();
519  c10->cd(3)->SetLogz();
520  HIST_Egamma1vsEgamma2->Draw("colz");
521  c10->cd(4)->SetLogz();
522  HISTinvariantMass_openingAngle->Draw("colz");
523  c10->cd(5)->SetLogz();
524  HIST_Egamma1Egamma2_openingAngle->Draw("colz");
525  c10->cd(6)->SetLogz();
526  HIST_Egamma1_openingAngle->Draw("colz");
527  c10->cd(7)->SetLogz();
528  HIST_Egamma1_Egamma1Egamma2->Draw("colz");
529  c10->cd(8)->SetLogz();
530  HISTinvariantMass_Egamma1Egamma2->Draw("colz");
531 
532  TCanvas* c100 = new TCanvas("c100","", 100, 100, 1200, 800);
533  HISTinvariantMass_openingAngle_Egamma1Egamma2->Draw("Box");
534  HISTinvariantMass_openingAngle_Egamma1Egamma2->GetXaxis()->SetTitle("m_{#gamma#gamma}");
535  HISTinvariantMass_openingAngle_Egamma1Egamma2->GetXaxis()->SetTitleOffset(1.5);
536  HISTinvariantMass_openingAngle_Egamma1Egamma2->GetYaxis()->SetTitle("opening angle");
537  HISTinvariantMass_openingAngle_Egamma1Egamma2->GetYaxis()->SetTitleOffset(1.5);
538  HISTinvariantMass_openingAngle_Egamma1Egamma2->GetZaxis()->SetTitle("E_{#gamma1} #times E_{#gamma2}");
539  HISTinvariantMass_openingAngle_Egamma1Egamma2->GetZaxis()->SetTitleOffset(1.1);
540 
542  TString outfile= "./invMass_histo_onlyFwEndCapIncluded_GetEnergyCorrected.root";
543  TFile* file_out = new TFile(outfile,"RECREATE");
544 
545  HISTinvariantMass->Write();
546  HISTinvariantMass_openingAngle->Write();
547  HIST_Egamma1Egamma2->Write();
548  HIST_Egamma1vsEgamma2->Write();
549  HIST_Egamma1Egamma2_openingAngle->Write();
550  HISTinvariantMass_Egamma1Egamma2->Write();
551  HIST_Egamma1_Egamma1Egamma2->Write();
552  HIST_Egamma1_openingAngle->Write();
553  HISTinvariantMass_openingAngle_Egamma1Egamma2->Write();
554 
555  file_out->Close();
556 
557  timer.Stop();
558  Double_t rtime = timer.RealTime();
559  Double_t ctime = timer.CpuTime();
560  printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
561  return 0;
562 }
mctrack
PndMCTrack * mctrack
Definition: analysis_point_hit_fwendcap.C:250
phi_mc
Double_t phi_mc
Definition: anal_hit_digi_cluster_fwendcap.C:31
tsim
TTree * tsim
Definition: anal_hit_digi_cluster_fwendcap.C:10
acos
friend F32vec4 acos(const F32vec4 &a)
Definition: P4_F32vec4.h:113
printf
printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime)
theta_mc
Double_t theta_mc
Definition: anal_hit_digi_cluster_fwendcap.C:31
crystal_did
Double_t crystal_did
Definition: anal_hit_digi_cluster_fwendcap.C:39
id_d
Int_t id_d
Definition: anal_hit_digi_cluster_fwendcap.C:43
PndEmcCluster::where
TVector3 where() const
Definition: PndEmcCluster.cxx:111
i
Int_t i
Definition: run_full.C:25
photon_momentum
TVector3 photon_momentum
Definition: anal_hit_digi_cluster_fwendcap.C:45
PndEmcCluster::GetModule
Short_t GetModule() const
Definition: PndEmcCluster.cxx:282
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
digi_array
TClonesArray * digi_array
Definition: anal_hit_digi_cluster_fwendcap.C:24
c10
TCanvas * c10
Definition: plot_invariantmass.C:42
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
ncounts
Int_t ncounts
Definition: anal_hit_digi_cluster_fwendcap.C:155
h1
TH1F * h1
Definition: anal_hit_digi_cluster_fwendcap.C:121
h10d
TH1F * h10d
Definition: anal_hit_digi_cluster_fwendcap.C:119
PndEmcCluster::z
Double_t z() const
Definition: PndEmcCluster.cxx:135
ene_mc
Double_t ene_mc
Definition: anal_hit_digi_cluster_fwendcap.C:31
digi_ene
Double_t digi_ene
Definition: anal_hit_digi_cluster_fwendcap.C:39
cluster_array
TClonesArray * cluster_array
Definition: anal_hit_digi_cluster_fwendcap.C:27
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
timer
TStopwatch timer
Definition: hit_dirc.C:51
cluster_energy
Double_t cluster_energy
Definition: anal_hit_digi_cluster_fwendcap.C:31
ndigi
int ndigi
Definition: anal_hit_digi_cluster_fwendcap.C:48
h20d
TH1F * h20d
Definition: anal_hit_digi_cluster_fwendcap.C:120
hene_mc
TH1F * hene_mc
Definition: anal_hit_digi_cluster_fwendcap.C:92
PndEmcCluster::x
Double_t x() const
Definition: PndEmcCluster.cxx:123
row_did
Double_t row_did
Definition: anal_hit_digi_cluster_fwendcap.C:39
PndEmcCluster
a cluster (group of neighboring crystals) of hit emc crystals
Definition: PndEmcCluster.h:29
h10
TH1F * h10
Definition: anal_hit_digi_cluster_fwendcap.C:117
ctime
Double_t ctime
Definition: hit_dirc.C:114
fsim
TFile * fsim
Definition: anal_hit_digi_cluster_fwendcap.C:8
digi_theta
Double_t digi_theta
Definition: anal_hit_digi_cluster_fwendcap.C:39
PndEmcCluster::NumberOfDigis
Int_t NumberOfDigis() const
Definition: PndEmcCluster.cxx:341
h2
TH1F * h2
Definition: anal_hit_digi_cluster_fwendcap.C:122
PndEmcCluster::GetEnergyCorrected
Double_t GetEnergyCorrected() const
Definition: PndEmcCluster.cxx:379
h20
TH1F * h20
Definition: anal_hit_digi_cluster_fwendcap.C:118
bump_array
TClonesArray * bump_array
Definition: bump_analys.C:14
mctrack_array
TClonesArray * mctrack_array
Definition: anal_hit_digi_cluster_fwendcap.C:11
p4mom
TLorentzVector p4mom
Definition: anal_hit_digi_cluster_fwendcap.C:46
rtime
Double_t rtime
Definition: hit_dirc.C:113
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
module_d
Double_t module_d
Definition: anal_hit_digi_cluster_fwendcap.C:39
digi_phi
Double_t digi_phi
Definition: anal_hit_digi_cluster_fwendcap.C:39
PndEmcCluster::y
Double_t y() const
Definition: PndEmcCluster.cxx:129
outfile
TString outfile
Definition: createRootFscGeometryFile.C:57
Generated on Sun Aug 18 2019 03:02:21 for FairRoot/PandaRoot by   doxygen 1.8.8