/[PAMELA software]/calo/flight/CaloElectron/src/CaloHough.cpp
ViewVC logotype

Annotation of /calo/flight/CaloElectron/src/CaloHough.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (hide annotations) (download)
Tue Aug 4 13:59:12 2009 UTC (15 years, 4 months ago) by mocchiut
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +3 -3 lines
Changed to work with GCC 4.x (gfortran) + ROOT >= 5.24

1 pam-fi 1.1 #include <CaloHough.h>
2     ////////////////////////////////////////////////////////////
3     /// INITIALIZATIONS & GLOBAL PARAMETERS
4     ////////////////////////////////////////////////////////////
5    
6     CaloHough_parameters * CaloHough_parameters::_parameters = 0;
7    
8     /**
9     * Set external parameters to default values
10     */
11     void CaloHough_parameters::SetDefault(){
12    
13     cout << "CaloHough --> SET DEFAULT PARAMETERS" << endl;
14    
15     isimu = 0;
16    
17     // Q_UP = 4.;// mip
18     // Q_DN = 0.;
19     // M_UP = 3;
20     // M_DN = 0;
21     Q_UP = 10000.;
22     Q_DN = 100.;
23     M_UP = 3000;
24     M_DN = 0;
25    
26     h_qtot[0] = new TH2F("houghqtot_0","Hit pattern",96,-0.5,95.5,22,-0.5,21.5);
27     h_qtot[1] = new TH2F("houghqtot_1","Hit pattern",96,-0.5,95.5,22,-0.5,21.5);
28    
29     int n1 = 150;
30     int n2 = 150;
31    
32     h_par[0] = new TH2F("houghpar_0","par1 vs par2",n1,-40.,40.,n2,-10.,10.);
33     h_par[1] = new TH2F("houghpar_1","par1 vs par2",n1,-40.,40.,n2,-10.,10.);
34    
35     h_par_fine[0] = new TH2F("houghparfine_0","par1 vs par2",10*n1,-40.,40.,10*n2,-10.,10.);
36     h_par_fine[1] = new TH2F("houghparfine_1","par1 vs par2",10*n1,-40.,40.,10*n2,-10.,10.);
37    
38     h_par1[0] = new TH1F("houghpar1_0","par1",n1,-10.,10.);
39     h_par1[1] = new TH1F("houghpar1_1","par1",n1,-10.,10.);
40    
41     h_par2[0] = new TH1F("houghpar2_0","par2",n2,-40.,40.);
42     h_par2[1] = new TH1F("houghpar2_1","par2",n2,-40.,40.);
43     };
44     ////////////////////////////////////////////////////////////
45     /// CLASS IMPLEMENTATION
46     ////////////////////////////////////////////////////////////
47     /**
48     * Reset event
49     * it has to be called for each event
50     */
51     void CaloHough::Reset(){
52     return;
53     }
54    
55     /**
56     * \brief Set event
57     */
58     bool CaloHough::Set(CaloLevel1 *cl1){
59    
60     // bool debug = CaloHough_parameters::Get()->debug;
61     float Q_UP = CaloHough_parameters::Get()->Q_UP;
62     int M_UP = CaloHough_parameters::Get()->M_UP;
63     float Q_DN = CaloHough_parameters::Get()->Q_DN;
64     int M_DN = CaloHough_parameters::Get()->M_DN;
65    
66     cout << "bool CaloHough::Set(CaloLevel1*)"<<endl;
67    
68     if(!cl1)return false;
69    
70     TH2F* h_qtot[2];
71     TH2F* h_par[2];
72     TH2F* h_par_fine[2];
73     TH1F* h_par1[2];
74     TH1F* h_par2[2];
75     for(int iv=0; iv<2; iv++){
76     h_qtot[iv] = CaloHough_parameters::Get()->h_qtot[iv];
77     h_par[iv] = CaloHough_parameters::Get()->h_par[iv];
78     h_par_fine[iv] = CaloHough_parameters::Get()->h_par_fine[iv];
79     h_par1[iv] = CaloHough_parameters::Get()->h_par1[iv];
80     h_par2[iv] = CaloHough_parameters::Get()->h_par2[iv];
81     if(!h_qtot[iv] || !h_par[iv] || !h_par1[iv] || !h_par2[iv] || !h_par_fine[iv])return false;
82     h_qtot[iv]->Reset();
83     h_par[iv]->Reset();
84     h_par_fine[iv]->Reset();
85     h_par1[iv]->Reset();
86     h_par2[iv]->Reset();
87     }
88    
89     // -----------------------------------
90     // set charge collected in each plane
91     // -----------------------------------
92     //
93     CaloStrip cstrip;
94     int isimu = CaloHough_parameters::Get()->isimu;
95     if(isimu==0)cstrip = CaloStrip(cl1,false);
96     if(isimu==1)cstrip = CaloStrip(cl1,true);
97     //
98     vector<int> ids[22][2]; //id strip 0-95
99     vector<float> xy[22][2]; //cluster coordinate
100     int islast = -1;
101     int iplast = -1;
102     int mult = 0;
103     float qsum = 0.;
104     float xyqsum = 0.;
105     float issum = 0.;
106     int ncluster[] = {0,0};
107     for(int ih=0;ih<cl1->istrip;ih++){
108     //
109     int iv=-1;
110     int ip=-1;
111     int is=-1;
112     cl1->DecodeEstrip(ih,iv,ip,is);
113     float xyz[3];
114     cstrip.Set(iv,ip,is);
115     // cout << " * "<<setw(3)<<iv<<setw(3)<<ip<<setw(3)<<is<<endl;
116     xyz[0] = cstrip.GetX();
117     xyz[1] = cstrip.GetY();
118     xyz[2] = cstrip.GetZ();
119     //
120     // ...................................... metodo 2
121     if( is != islast+1 || ip != iplast ){ // new cluster
122     if(
123     qsum < Q_UP &&
124     qsum > Q_DN &&
125     mult < M_UP &&
126     mult > M_DN &&
127     true){ // store the previous cluster
128     xy[ip][iv].push_back( xyqsum/qsum );
129     ids[ip][iv].push_back( (int)(0.5+issum/qsum) );
130     h_qtot[iv]->Fill( issum/qsum, (float)(21-ip), qsum);
131     ncluster[iv]++;
132     cout << " >> "<<qsum<<endl;
133     }// ... and reset counters
134     mult = 0;
135     qsum = 0.;
136     xyqsum = 0.;
137     issum = 0.;
138     }
139    
140     mult++;
141     qsum += cstrip.GetE();
142     xyqsum += xyz[iv] * cstrip.GetE();
143     issum += (float)is * cstrip.GetE();
144     islast = is;
145     iplast = ip;
146    
147     // ...................................... metodo 1
148     // if( cstrip.GetE() < Q_UP ){
149     // h_qtot[iv]->Fill( (float)is, (float)(21-ip), cstrip.GetE());
150     // ids[ip][iv].push_back(is);
151     // }
152    
153     }
154     cout << "n.clusters X "<<ncluster[0]<<endl;
155     cout << "n.clusters Y "<<ncluster[0]<<endl;
156     //
157     // cstrip.Set(0,0,0);
158     // cstrip.Set(0,21,0);
159     cstrip.Set(0,11,0);
160     float zref = cstrip.GetZ();//reference plane
161     //
162     for(int iv=0; iv<2; iv++){
163     for(int ip1=0; ip1<22-1; ip1++){
164     // for(int ip2=ip1+1; ip2<22; ip2++){
165     for(int ip2=ip1+1; ip2<min(22,ip1+3); ip2++){
166     // cout << "ip1 "<<ip1<<" "<<ids[ip1][iv].size()<<endl;
167     // cout << "ip2 "<<ip2<<" "<<ids[ip2][iv].size()<<endl;
168 mocchiut 1.2 for(int is1=0; is1<(int)ids[ip1][iv].size(); is1++){
169     for(int is2=0; is2<(int)ids[ip2][iv].size(); is2++){
170 pam-fi 1.1
171     // ...................................... metodo 1
172     // int ids1 = ids[ip1][iv][is1];
173     // int ids2 = ids[ip2][iv][is2];
174     // float xyz1[3];
175     // float xyz2[3];
176     // cstrip.Set(iv,ip1,ids1);
177     // xyz1[0] = cstrip.GetX();
178     // xyz1[1] = cstrip.GetY();
179     // xyz1[2] = cstrip.GetZ();
180     // cstrip.Set(iv,ip2,ids2);
181     // xyz2[0] = cstrip.GetX();
182     // xyz2[1] = cstrip.GetY();
183     // xyz2[2] = cstrip.GetZ();
184    
185     // ...................................... metodo 2
186 mocchiut 1.2 // int ids1 = ids[ip1][iv][is1];
187 pam-fi 1.1 int ids2 = ids[ip2][iv][is2];
188     float xyz1[3];
189     float xyz2[3];
190     cstrip.Set(iv,ip1,is1);
191     xyz1[iv] = xy[ip1][iv][is1];
192     xyz1[2] = cstrip.GetZ();
193     cstrip.Set(iv,ip2,ids2);
194     xyz2[iv] = xy[ip2][iv][is2];
195     xyz2[2] = cstrip.GetZ();
196     //
197     float par1 = 0.;
198     if(fabs(xyz1[2]-xyz2[2])>0.)par1 = (xyz1[iv]-xyz2[iv])/(xyz1[2]-xyz2[2]);
199     float par2 = xyz1[iv] - par1*xyz1[2];
200     par2 = par1*zref + par2;
201     h_par[iv]->Fill(par2,par1); // parameter space
202     h_par_fine[iv]->Fill(par2,par1); // parameter space
203     h_par1[iv]->Fill(par1);
204     h_par2[iv]->Fill(par2);
205     // cout << setw(15)<<par1<< setw(15)<<par2<<endl;
206     }
207     }
208     }
209     }
210     //
211     for( int ibx=1; ibx <= h_par[iv]->GetXaxis()->GetNbins(); ibx++){
212     for( int iby=1; iby <= h_par[iv]->GetYaxis()->GetNbins(); iby++){
213     // get bin content
214     float content = h_par[iv]->GetBinContent(ibx,iby);
215     if( content > 6){
216     double xmi = h_par[iv]->GetXaxis()->GetBinLowEdge(ibx);
217     double xma = h_par[iv]->GetXaxis()->GetBinUpEdge(ibx);
218     double ymi = h_par[iv]->GetYaxis()->GetBinLowEdge(iby);
219     double yma = h_par[iv]->GetYaxis()->GetBinUpEdge(iby);
220     h_par_fine[iv]->GetXaxis()->SetRangeUser(xmi,xma);
221     h_par_fine[iv]->GetYaxis()->SetRangeUser(ymi,yma);
222     double par1 = h_par_fine[iv]->GetMean(1);
223     double par2 = h_par_fine[iv]->GetMean(2);
224     cout << iv << " <><><><><><><>< "<<setw(15)<<par1<<setw(15)<<par2<<endl;
225     }
226     }
227     }
228    
229     //
230     for(int ip=0; ip<22; ip++){
231     ids[ip][iv].clear();
232     xy[ip][iv].clear();
233     }
234     }//end loop on views
235    
236     return true;
237     }
238     ClassImp(CaloHough);
239     ClassImp(CaloHough_parameters);
240    

  ViewVC Help
Powered by ViewVC 1.1.23