All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ICARUSPurityDQM_module.cc
Go to the documentation of this file.
1 
2 ////////////////////////////////////////////////////////////////////////
3 //
4 // ICARUSPurityDQM class
5 //
6 // Christian Farnese
7 //
8 ////////////////////////////////////////////////////////////////////////
9 
10 //#ifndef ICARUSPURITYDQM_H
11 //#define ICARUSPURITYDQM_H
12 
13 #include <iomanip>
14 #include <TH1F.h>
15 #include <TProfile.h>
16 #include <vector>
17 #include <string>
18 #include <array>
19 #include <fstream>
20 
21 //Framework includes
22 #include "art/Framework/Core/ModuleMacros.h"
23 #include "canvas/Persistency/Common/FindManyP.h"
24 #include "art/Framework/Principal/Event.h"
25 #include "fhiclcpp/ParameterSet.h"
26 #include "art/Framework/Principal/Handle.h"
27 #include "canvas/Persistency/Common/Ptr.h"
28 #include "canvas/Persistency/Common/PtrVector.h"
29 //#include "art/Framework/Services/Registry/ServiceHandle.h"
30 //#include "art/Framework/Services/Optional/TFileService.h"
31 #include "art_root_io/TFileService.h"
32 //#include "art/Framework/Services/Optional/TFileDirectory.h"
33 //#include "messagefacility/MessageLogger/MessageLogger.h"
34 
35 //LArSoft includes
36 #include "larcore/Geometry/Geometry.h"
37 #include "nusimdata/SimulationBase/MCTruth.h"
38 #include "nug4/ParticleNavigation/ParticleList.h"
39 #include "nug4/ParticleNavigation/EmEveIdCalculator.h"
40 #include "lardataobj/RecoBase/Hit.h"
41 #include "lardataobj/RecoBase/Cluster.h"
42 #include "lardataobj/RecoBase/Wire.h"
43 #include "lardataobj/RecoBase/Vertex.h"
44 #include "larsim/MCCheater/BackTracker.h"
45 #include "lardata/Utilities/AssociationUtil.h"
46 
47 #include "lardataobj/RawData/RawDigit.h"
48 #include "lardataobj/RawData/raw.h"
49 
50 //Database Connection Files
51 //#include "../../MetricManagerShim/MetricManager.hh"
52 //#include "../../MetricConfig/ConfigureRedis.hh"
53 
54 //purity info class
55 #include "sbnobj/Common/Analysis/TPCPurityInfo.hh"
56 
57 #include "art/Framework/Core/EDProducer.h"
58 #include <TMath.h>
59 #include <TH1F.h>
60 #include "TH2D.h"
61 #include "TProfile2D.h"
62 #include <TGraph.h>
63 #include <TGraphErrors.h>
64 #include <TGraphAsymmErrors.h>
65 #include "TF1.h"
66 #include "TCanvas.h"
67 #include "TNtuple.h"
68 
69 //Redis connection
70 //#include "sbndaq-redis-plugin/Utilities.h"
71 
72 class TH1F;
73 class TH2F;
74 
75 
76 ///Cluster finding and building
77 namespace icarus {
78 
79  class ICARUSPurityDQM : public art::EDProducer {
80 
81  public:
82 
83  explicit ICARUSPurityDQM(fhicl::ParameterSet const& pset);
84  virtual ~ICARUSPurityDQM();
85 
86  /// read access to event
87  void produce(art::Event& evt);
88  void beginJob();
89  void endJob();
90  int Nothere(std::vector<int>* a, int b);
91  int Notheref(std::vector<float>* a, float b);
92  Double_t FoundMeanLog(std::vector<float>* a,float b);
93 
94  private:
95 
96  TH1F* puritytpc0;
97  TH1F* puritytpc1;
98  TH1F* puritytpc2;
99  TH1F* puritytpc3;
100 
101 
102 
103 
104  TH1F* purityvalues;
105  TH1F* purityvalues2;
106  TH1F* h_hit_height;
107  TH1F* h_hit_area;
108  TH2D* h_hit_height_area;
109 
110  TH1F* h_basediff;
111  TH1F* h_base1;
112  TH1F* h_base2;
113  TH1F* h_basebase;
114  TH1F* h_ratio;
115  TH1F* h_ratio_after;
116  TH1F* h_ratio_after_2;
117  TH2D* h_ratio_after_3;
118  TH2D* h_ratio_3;
119 
120 
121  TH1F* h_rms;
122  TH1F* h_errors;
123  TH1F* h_hittime;
124  TH1F* h_hittime_2;
125  TH1F* h_hittime_3;
126 
127  //TH1F* fRun;
128  //TH2D* fRunSub;
129 
130  std::vector<art::InputTag> fDigitModuleLabel;
131  //short fPrintLevel;
132 
133  float fValoretaufcl;
134  float fmaxfcl;
135  float fminfcl;
136  float fdisfcl;
137  int fcryofcl;
138  int fquantevoltefcl;
139 
140  int fplanefcl;
141  int fthresholdfcl;
142  int fdumphitsfcl;
143  int fgruppifcl;
144 int fdwclusfcl;
145 int fdsclusfcl;
146  bool fPersistPurityInfo;
147 
148  TNtuple* purityTuple;
149  bool fFillAnaTuple;
150 
151  TTree* fpurTree;
152  int run_tree;
153  int subrun_tree;
154  int event_tree;
155  int tpc_tree;
156  int qhits_tree;
157  float dwire_tree;
158  float dtime_tree;
159  float earea_tree;
160  float marea_tree;
161  float slope_tree;
162  float errslope_tree;
163  float chi_tree;
164 
165 
166  }; // class ICARUSPurityDQM
167 
168 }
169 
170 //#endif
171 
172 namespace icarus{
173 
174  //--------------------------------------------------------------------
175  ICARUSPurityDQM::ICARUSPurityDQM(fhicl::ParameterSet const& pset)
176  : EDProducer(pset)
177  , fDigitModuleLabel (pset.get< std::vector<art::InputTag> > ("RawModuleLabel"))
178  , fValoretaufcl (pset.get< float > ("ValoreTauFCL"))
179  , fmaxfcl (pset.get< float > ("FracMaxFCL",0.9))
180  , fminfcl (pset.get< float > ("FraxMinFCL",0.05))
181  , fdisfcl (pset.get< float > ("DisFCL",3.0))
182  , fcryofcl (pset.get< int > ("CryostatFCL"))
183  , fquantevoltefcl (pset.get< int > ("QuanteFCL",5))
184  , fplanefcl (pset.get< int > ("PlaneFCL"))
185  , fthresholdfcl (pset.get< int > ("ThresholdFCL"))
186  , fdumphitsfcl (pset.get< int > ("DumpHitFCL",0))
187  , fgruppifcl (pset.get< int > ("GruppiFCL",8))
188  , fdwclusfcl (pset.get< int > ("DeltaWClusFCL",3))
189  , fdsclusfcl (pset.get< int > ("DeltaSClusFCL",100))
190  , fPersistPurityInfo (pset.get< bool > ("PersistPurityInfo",false))
191  , fFillAnaTuple (pset.get< bool > ("FillAnaTuple",false))
192  {
193 
194  //declare what we produce .. allow it to not be persistable to the event
195  if(fPersistPurityInfo)
196  produces< std::vector<anab::TPCPurityInfo> >("",art::Persistable::Yes);
197  else
198  produces< std::vector<anab::TPCPurityInfo> >("",art::Persistable::No);
199 
200  }
201 
202  //------------------------------------------------------------------
203  ICARUSPurityDQM::~ICARUSPurityDQM()
204  {
205 
206  }
207 
208  void ICARUSPurityDQM::beginJob()
209  {
210 
211  // get access to the TFile service
212  art::ServiceHandle<art::TFileService> tfs;
213 
214  purityvalues = tfs->make<TH1F>("purityvalues","purityvalues",20000,-10,10);
215  h_basediff = tfs->make<TH1F>("h_basediff","h_basediff",10000,-20,20);
216  h_basebase = tfs->make<TH1F>("h_basebase","h_basebase",10000,-20,20);
217  h_base1 = tfs->make<TH1F>("h_base1","h_base1",10000,-20,20);
218  h_base2 = tfs->make<TH1F>("h_base2","h_base2",10000,-20,20);
219 
220  h_rms = tfs->make<TH1F>("h_rms","h_rms",2000,0,20);
221  //fRun=tfs->make<TH1F>("fRun","Events per run", 4000,0.5 ,4000.5);
222  //fRunSub=tfs->make<TH2D>("fRunSub","Events per run", 4000,0.5 ,4000.5,50,0.5,50.5);
223  h_errors = tfs->make<TH1F>("h_errors","",1000,0,1000);
224 
225 
226  h_hittime = tfs->make<TH1F>("h_hittime","",2500,-0.5,2499.5);
227  h_hittime_2 = tfs->make<TH1F>("h_hittime_2","",2500,-0.5,2499.5);
228  h_hittime_3 = tfs->make<TH1F>("h_hittime_3","",2500,-0.5,2499.5);
229 
230  purityvalues2 = tfs->make<TH1F>("purityvalues2","purityvalues2",20000,-10,10);
231  puritytpc0 = tfs->make<TH1F>("puritytpc0","puritytpc0",20000,-10,10);
232  puritytpc1 = tfs->make<TH1F>("puritytpc1","puritytpc1",20000,-10,10);
233  puritytpc2 = tfs->make<TH1F>("puritytpc2","puritytpc2",20000,-10,10);
234  puritytpc3 = tfs->make<TH1F>("puritytpc3","puritytpc3",20000,-10,10);
235  h_hit_height = tfs->make<TH1F>("h_hit_height","h_hit_height",2000,0.5,2000.5);
236  h_hit_area = tfs->make<TH1F>("h_hit_area","h_hit_area",10000,0,10000);
237  h_hit_height_area=tfs->make<TH2D>("h_hit_height_area","", 300,0.5,300.5,2000,0.,2000.);
238 
239  h_ratio = tfs->make<TH1F>("h_ratio","h_ratio",20000,0,5);
240  h_ratio_after = tfs->make<TH1F>("h_ratio_after","h_ratio_after",20000,0,5);
241  h_ratio_after_2 = tfs->make<TH1F>("h_ratio_after_2","",2000,0,2);
242  h_ratio_after_3 = tfs->make<TH2D>("h_ratio_after_3","",1300,0,2600,2000,0,2);
243  h_ratio_3 = tfs->make<TH2D>("h_ratio_3","",1300,0,2600,2000,0,2);
244 
245 
246 
247  fpurTree = tfs->make<TTree>("puritytree","tree for the purity studies");
248  fpurTree->Branch("run",&run_tree,"run/I");
249  fpurTree->Branch("subrun",&subrun_tree,"subrun/I");
250  fpurTree->Branch("event",&event_tree,"event/I");
251  fpurTree->Branch("tpc",&tpc_tree,"tpc/I");
252  fpurTree->Branch("qhits",&qhits_tree,"qhits/I");
253  fpurTree->Branch("deltawire",&dwire_tree,"deltawire/F");
254  fpurTree->Branch("deltatime",&dtime_tree,"deltatime/F");
255  fpurTree->Branch("errarea",&earea_tree,"errarea/F");
256  fpurTree->Branch("meanarea",&marea_tree,"meanarea/F");
257  fpurTree->Branch("slope",&slope_tree,"slope/F");
258  fpurTree->Branch("errorslope",&errslope_tree,"errorslope/F");
259  fpurTree->Branch("chisquare",&chi_tree,"chisquare/F");
260 
261 
262  if(fFillAnaTuple)
263  purityTuple = tfs->make<TNtuple>("purityTuple","Purity Tuple","run:ev:tpc:att");
264 
265  }
266 
267  void ICARUSPurityDQM::endJob()
268  {
269  std::ofstream goodpurofinal("valore_indicativo.out",std::ios::app);
270  goodpurofinal << purityvalues2->GetMean() << std::endl;
271  //if(fPrintLevel == -1) outFile.close();
272  }
273 
274 
275  int ICARUSPurityDQM::Nothere(std::vector<int>* a, int b){
276  int Cisono=3;
277  for(int i=0; i<(int)a->size();i++)
278  {
279  if((*a)[i]==b)
280  {
281  Cisono=8;
282  }
283  }
284  return Cisono;//se "ci sono" vale 8 allora l'intero b e' contenuto nel vettore a mentre se "ci sono" vale 3 a non contiene b.
285  }
286 
287 
288  int ICARUSPurityDQM::Notheref(std::vector<float>* a, float b){
289  int Cisono=3;
290  for(int i=0; i<(int)a->size();i++)
291  {
292  if((*a)[i]==b)
293  {
294  Cisono=8;
295  }
296  }
297 
298  return Cisono;//se "ci sono" vale 8 allora l'intero b e' contenuto nel vettore a mentre se "ci sono" vale 3 a non contiene b.
299  }
300 
301  Double_t ICARUSPurityDQM::FoundMeanLog(std::vector<float>* a,float b){
302 
303 
304  float taglio=0;
305  if(a->size()>0){
306  int punto_taglio=a->size()*(b)+0.5;
307  //cout << punto_taglio << " e " << a->size() << endl;
308  if(punto_taglio==0)punto_taglio=1;
309 std::vector<float>* usedhere=new std::vector<float>;
310  for(int j=0;j<(int)a->size();j++)usedhere->push_back((*a)[j]);
311 
312  //cout << (*usedhere)[punto_taglio-1] << " PRIMA " << endl;
313 std::sort(usedhere->begin(),usedhere->end(), [](float &c, float &d){ return c<d; });
314  //cout << (*usedhere)[punto_taglio-1] << " DOPO " << endl;
315  taglio=(*usedhere)[punto_taglio-1];
316  delete usedhere;
317  }
318  return taglio;
319 
320 /*
321  int punto_taglio=a->size()*(1-b)+0.5;
322  ///////std::cout << punto_taglio << " e " << a->size() << std::endl;
323  //if(punto_taglio==0)punto_taglio=1;
324  std::vector<float>* usedhere=new std::vector<float>;
325  for(int jj=0;jj<punto_taglio+1;jj++)
326  {
327  //cout << jj << endl;
328  float maximum=0;
329  for(int j=0;j<(int)a->size();j++)
330  {
331  if((*a)[j]>maximum && Notheref(usedhere,(*a)[j])==3)
332  {
333  maximum=(*a)[j];
334  }
335  }
336  //cout << maximum << endl;
337  usedhere->push_back(maximum);
338  }
339  //cout << (*usedhere)[punto_taglio] << endl;
340  return (*usedhere)[punto_taglio-1];
341 */
342  }
343 
344  void ICARUSPurityDQM::produce(art::Event& evt)
345  {
346 
347 // //std::cout << " Inizia Purity ICARUS Ana - upgraded by C.FARNESE, WES and OLIVIA " << std::endl;
348  // code stolen from TrackAna_module.cc
349  art::ServiceHandle<geo::Geometry> geom;
350  unsigned int fDataSize;
351  std::vector<short> rawadc; //UNCOMPRESSED ADC VALUES.
352  // get all hits in the event
353  //InputTag cluster_tag { "fuzzycluster" }; //CH comment trovato con eventdump code
354 
355  //to get run and event info, you use this "eventAuxillary()" object.
356  //art::Timestamp ts = evt.time();
357 // //std::cout << "Processing for Purity " << " Run " << evt.run() << ", " << "Event " << evt.event() << " and Time " << evt.time().timeHigh() << std::endl;
358 
359  //fRun->Fill(evt.run());
360  //fRunSub->Fill(evt.run(),evt.subRun());
361  art::Handle< std::vector<raw::RawDigit> > digitVecHandle;
362  std::vector<const raw::RawDigit*> rawDigitVec;
363 
364  //setup output vector
365  std::unique_ptr< std::vector<anab::TPCPurityInfo> > outputPtrVector(new std::vector<anab::TPCPurityInfo>() );
366 
367  anab::TPCPurityInfo purity_info;
368  purity_info.Run = evt.run();
369  purity_info.Subrun = evt.subRun(); //evt.time().timeHigh()-1598580000;//evt.subRun();
370  purity_info.Event = evt.event();
371  std::ofstream purh("dump_purity_hits.out",std::ios::app);
372 
373  for(const auto& digitlabel : fDigitModuleLabel)
374  {
375  evt.getByLabel(digitlabel, digitVecHandle);
376  std::vector<const raw::RawDigit*> rawDigitVec;
377 
378  if (digitVecHandle.isValid())
379  {
380  // Sadly, the RawDigits come to us in an unsorted condition which is not optimal for
381  // what we want to do here. So we make a vector of pointers to the input raw digits and sort them
382  // Ugliness to fill the pointer vector...
383  for(size_t idx = 0; idx < digitVecHandle->size(); idx++)
384  rawDigitVec.push_back(&digitVecHandle->at(idx));
385  // Sort (use a lambda to sort by channel id)
386  std::sort(rawDigitVec.begin(),
387  rawDigitVec.end(),
388  [](const raw::RawDigit* left, const raw::RawDigit* right) {return left->Channel() < right->Channel();});
389  }
390 
391 
392  std::vector<int> *www0=new std::vector<int>;
393  std::vector<float> *sss0=new std::vector<float>;
394  std::vector<float> *hhh0=new std::vector<float>;
395  std::vector<float> *ehh0=new std::vector<float>;
396  std::vector<int> *www2=new std::vector<int>;
397  std::vector<float> *sss2=new std::vector<float>;
398  std::vector<float> *hhh2=new std::vector<float>;
399  std::vector<float> *ehh2=new std::vector<float>;
400  std::vector<int> *ccc0=new std::vector<int>;
401  std::vector<int> *ccc2=new std::vector<int>;
402  std::vector<float> *aaa0=new std::vector<float>;
403  std::vector<float> *aaa2=new std::vector<float>;
404 
405 
406  short int used[4096];
407  for(const auto& rawDigit : rawDigitVec)
408  {
409  raw::ChannelID_t channel = rawDigit->Channel();
410  ////std::cout << channel << std::endl;
411  std::vector<geo::WireID> wids = geom->ChannelToWire(channel);
412  // for now, just take the first option returned from ChannelToWire
413  geo::WireID wid = wids[0];
414  // We need to know the plane to look up parameters
415  geo::PlaneID::PlaneID_t plane = wid.Plane;
416  //size_t
417 
418  int cryostat=wid.Cryostat;
419  size_t tpc=wid.TPC;
420  size_t iWire=wid.Wire;
421  ////std::cout << channel << " INFO CHANNEL " << plane << " " << tpc << " " << cryostat << " " << iWire << std::endl;
422  fDataSize = rawDigit->Samples();
423  rawadc.resize(fDataSize);
424 
425  for(int ijk=0;ijk<4096;ijk++)used[ijk]=0;
426  //UNCOMPRESS THE DATA.
427  int pedestal = (int)rawDigit->GetPedestal();
428  float pedestal2 = rawDigit->GetPedestal();
429  float sigma_pedestal = rawDigit->GetSigma();
430  raw::Uncompress(rawDigit->ADCs(), rawadc, pedestal, rawDigit->Compression());
431  ////std::cout << pedestal2 << " " << fDataSize << " " << rawadc.size() << " " << sigma_pedestal << std::endl;
432  float massimo=0;
433  float quale_sample_massimo;
434  // if (plane==0) {
435  if ((int)plane==fplanefcl && cryostat==fcryofcl){
436  for(int volte=0;volte<fquantevoltefcl;volte++){
437  massimo=0;
438  quale_sample_massimo=-1;
439 
440  TH1F *h111 = new TH1F("h111","delta aree",20000,0,0);
441  for (unsigned int ijk=0; ijk<(fDataSize); ijk++)
442  {
443  if ((rawDigit->ADC(ijk)-pedestal2)>massimo && ijk>150 && ijk<(fDataSize-150) && used[ijk]==0)
444  {
445  massimo=(rawDigit->ADC(ijk)-pedestal2);
446  quale_sample_massimo=ijk;
447  }
448  }
449  float base_massimo_before=0;
450  float base_massimo_after=0;
451  for (unsigned int ijk=quale_sample_massimo-135; ijk<quale_sample_massimo-35; ijk++)
452  {
453  base_massimo_before+=(rawDigit->ADC(ijk)-pedestal2)*0.01;
454  }
455  for (unsigned int ijk=quale_sample_massimo+35; ijk<quale_sample_massimo+135; ijk++)
456  {
457  base_massimo_after+=(rawDigit->ADC(ijk)-pedestal2)*0.01;
458  }
459  float basebase=(base_massimo_after+base_massimo_before)*0.5;
460  h_basediff->Fill(fabs(base_massimo_after-base_massimo_before)); h_base1->Fill(base_massimo_before);
461  h_base2->Fill(base_massimo_after); h_basebase->Fill(basebase);
462  float areaarea=0;
463  for (unsigned int ijk=0; ijk<(fDataSize); ijk++)
464  {
465  if (ijk>(quale_sample_massimo-30) && ijk<(quale_sample_massimo+30)) {
466  if(fabs(base_massimo_after-base_massimo_before)<=sigma_pedestal)areaarea+=(rawDigit->ADC(ijk)-pedestal2-basebase);
467  if(fabs(base_massimo_after-base_massimo_before)>sigma_pedestal && base_massimo_after<base_massimo_before)areaarea+=(rawDigit->ADC(ijk)-pedestal2-base_massimo_after);
468  if(fabs(base_massimo_after-base_massimo_before)>sigma_pedestal && base_massimo_after>base_massimo_before)areaarea+=(rawDigit->ADC(ijk)-pedestal2-base_massimo_before);
469 
470  }
471  else{
472  if(used[ijk]==0)h111->Fill(rawDigit->ADC(ijk)-pedestal2-basebase);
473  }
474  }
475  sigma_pedestal=h111->GetRMS();
476  h111->Delete();
477  h_rms->Fill(sigma_pedestal);
478  if(fdumphitsfcl==1)purh<<evt.event()<< " " <<iWire<<" "<<tpc<<" "<<cryostat<<" "<<basebase<<" "<<base_massimo_after<<" " <<base_massimo_before<<" "<<pedestal2<<" "<<sigma_pedestal<< " " << quale_sample_massimo << " " << massimo << std::endl;
479 
480  if (massimo>(fthresholdfcl*sigma_pedestal)/* && fabs(base_massimo_after-base_massimo_before)<sigma_pedestal*/)
481  {
482  if(fdumphitsfcl==1)purh<<iWire<<" "<<tpc<<" "<<cryostat<< " CANDIDATE HIT " << quale_sample_massimo << " " << massimo << std::endl;
483 
484  if(tpc==0)www0->push_back(iWire+64);
485  if(tpc==0)sss0->push_back(quale_sample_massimo);
486  if(tpc==0)hhh0->push_back(massimo);
487  if(tpc==0)ehh0->push_back(sigma_pedestal);
488  if(tpc==0)ccc0->push_back(-1);
489  if(tpc==1)www0->push_back(iWire+64+2536);
490  if(tpc==1)sss0->push_back(quale_sample_massimo);
491  if(tpc==1)hhh0->push_back(massimo);
492  if(tpc==1)ehh0->push_back(sigma_pedestal);
493  if(tpc==1)ccc0->push_back(-1);
494  if(tpc==2)www2->push_back(iWire+64);
495  if(tpc==2)sss2->push_back(quale_sample_massimo);
496  if(tpc==2)hhh2->push_back(massimo);
497  if(tpc==2)ehh2->push_back(sigma_pedestal);
498  if(tpc==2)ccc2->push_back(-1);
499  if(tpc==3)www2->push_back(iWire+64+2536);
500  if(tpc==3)sss2->push_back(quale_sample_massimo);
501  if(tpc==3)hhh2->push_back(massimo);
502  if(tpc==3)ehh2->push_back(sigma_pedestal);
503  if(tpc==3)ccc2->push_back(-1);
504 
505  if(tpc==0)aaa0->push_back(areaarea);
506  if(tpc==1)aaa0->push_back(areaarea);
507  if(tpc==2)aaa2->push_back(areaarea);
508  if(tpc==3)aaa2->push_back(areaarea);
509  for(int ijk=0;ijk<330;ijk++)
510  {
511  int ent_value=quale_sample_massimo+ijk-165;
512  used[ent_value]=1;
513  }
514 
515  }
516  }
517  }
518  }
519 ////std::cout << "SIZE " << www0->size() << " " << www1->size() << " " << www2->size() << " " << www3->size() << std::endl;
520  for (unsigned int ijk=0; ijk<www0->size(); ijk++)
521  {
522  for (unsigned int ijk2=0; ijk2<www0->size(); ijk2++)
523  {
524  if (fabs((*www0)[ijk]-(*www0)[ijk2])<fdwclusfcl && ijk!=ijk2 && fabs((*sss0)[ijk]-(*sss0)[ijk2])<fdsclusfcl)
525  {
526  if ((*ccc0)[ijk]<0 && (*ccc0)[ijk2]<0)
527  {
528  (*ccc0)[ijk]=ijk+1;
529  (*ccc0)[ijk2]=ijk+1;
530  }
531  else if ((*ccc0)[ijk]>0 && (*ccc0)[ijk2]<0)
532  {
533  (*ccc0)[ijk2]=(*ccc0)[ijk];
534  }
535  else if ((*ccc0)[ijk]<0 && (*ccc0)[ijk2]>0)
536  {
537  (*ccc0)[ijk]=(*ccc0)[ijk2];
538  }
539  else if ((*ccc0)[ijk]>0 && (*ccc0)[ijk2]>0)
540  {
541  for (unsigned int ijk3=0; ijk3<www0->size(); ijk3++)
542  {
543  if(((*ccc0)[ijk3])==((*ccc0)[ijk2]))(*ccc0)[ijk3]=(*ccc0)[ijk];
544  }
545  }
546  }
547  }
548  }
549 
550 
551  for (unsigned int ijk=0; ijk<www2->size(); ijk++)
552  {
553  ////std::cout << " WWW2 " << ijk << std::endl;
554  for (unsigned int ijk2=0; ijk2<www2->size(); ijk2++)
555  {
556  if (fabs((*www2)[ijk]-(*www2)[ijk2])<fdwclusfcl && ijk!=ijk2 && fabs((*sss2)[ijk]-(*sss2)[ijk2])<fdsclusfcl)
557  {
558  if ((*ccc2)[ijk]<0 && (*ccc2)[ijk2]<0)
559  {
560  (*ccc2)[ijk]=ijk+1;
561  (*ccc2)[ijk2]=ijk+1;
562  }
563  else if ((*ccc2)[ijk]>0 && (*ccc2)[ijk2]<0)
564  {
565  (*ccc2)[ijk2]=(*ccc2)[ijk];
566  }
567  else if ((*ccc2)[ijk]<0 && (*ccc2)[ijk2]>0)
568  {
569  (*ccc2)[ijk]=(*ccc2)[ijk2];
570  }
571  else if ((*ccc2)[ijk]>0 && (*ccc2)[ijk2]>0)
572  {
573  for (unsigned int ijk3=0; ijk3<www2->size(); ijk3++)
574  {
575  if(((*ccc2)[ijk3])==((*ccc2)[ijk2]))(*ccc2)[ijk3]=(*ccc2)[ijk];
576  }
577  }
578  }
579  }
580  }
581 
582 
583  Int_t clusters_creation[4][30000];
584  //Int_t clusters_avewire[4][1000];
585  Int_t clusters_swire[4][30000];
586  Int_t clusters_lwire[4][30000];
587  Int_t clusters_ssample[4][30000];
588  Int_t clusters_lsample[4][30000];
589 
590  Int_t clusters_nn[1000];
591  Int_t clusters_vi[1000];
592  Int_t clusters_qq[1000];
593  //Int_t clusters_avewire[4][1000];
594  Int_t clusters_dw[1000];
595  Int_t clusters_ds[1000];
596  Double_t clusters_mw[1000];
597  Double_t clusters_ms[1000];
598  Double_t clusters_mintime[1000];
599 
600 
601  for (unsigned int ijk=0; ijk<4; ijk++) {
602  for (unsigned int ijk2=0; ijk2<30000; ijk2++) {
603  clusters_creation[ijk][ijk2]=0;
604  clusters_swire[ijk][ijk2]=100000;
605  clusters_lwire[ijk][ijk2]=0;
606  clusters_ssample[ijk][ijk2]=100000;
607  clusters_lsample[ijk][ijk2]=0;
608  if (ijk2<1000 && ijk>2) {
609  clusters_nn[ijk2]=-1;
610  clusters_vi[ijk2]=-1;
611  clusters_dw[ijk2]=-1;
612  clusters_ds[ijk2]=-1;
613  clusters_qq[ijk2]=-1;
614  clusters_mw[ijk2]=-1;
615  clusters_ms[ijk2]=-1;
616  clusters_mintime[ijk2]=-1;
617  }
618  }
619  }
620 
621  for (unsigned int ijk=0; ijk<www0->size(); ijk++)
622  {
623  if ((*ccc0)[ijk]>0) {
624  int numero=(*ccc0)[ijk];
625  clusters_creation[0][numero]+=1;
626  if((*www0)[ijk]<clusters_swire[0][numero])clusters_swire[0][numero]=(*www0)[ijk];
627  if((*www0)[ijk]>clusters_lwire[0][numero])clusters_lwire[0][numero]=(*www0)[ijk];
628  if((*sss0)[ijk]<clusters_ssample[0][numero])clusters_ssample[0][numero]=(*sss0)[ijk];
629  if((*sss0)[ijk]>clusters_lsample[0][numero])clusters_lsample[0][numero]=(*sss0)[ijk];
630  }
631  }
632  for (unsigned int ijk=0; ijk<www2->size(); ijk++)
633  {
634  if ((*ccc2)[ijk]>0) {
635  int numero=(*ccc2)[ijk];
636  clusters_creation[2][numero]+=1;
637  if((*www2)[ijk]<clusters_swire[2][numero])clusters_swire[2][numero]=(*www2)[ijk];
638  if((*www2)[ijk]>clusters_lwire[2][numero])clusters_lwire[2][numero]=(*www2)[ijk];
639  if((*sss2)[ijk]<clusters_ssample[2][numero])clusters_ssample[2][numero]=(*sss2)[ijk];
640  if((*sss2)[ijk]>clusters_lsample[2][numero])clusters_lsample[2][numero]=(*sss2)[ijk];
641  }
642  }
643 
644  int quanti_clusters=0;
645  for (unsigned int ijk=0; ijk<4; ijk++) {
646  for (unsigned int ijk2=0; ijk2<30000; ijk2++) {
647  if(clusters_creation[ijk][ijk2]>50)
648  {
649  ////std::cout<<clusters_creation[ijk][ijk2] << " CLUSTER MINE " << clusters_swire[ijk][ijk2] << " " << clusters_lwire[ijk][ijk2] << " " << clusters_ssample[ijk][ijk2] << " " << clusters_lsample[ijk][ijk2] << std::endl;
650  clusters_qq[quanti_clusters]=clusters_creation[ijk][ijk2];
651  clusters_vi[quanti_clusters]=ijk;
652  clusters_nn[quanti_clusters]=ijk2;
653  clusters_dw[quanti_clusters]=clusters_lwire[ijk][ijk2]-clusters_swire[ijk][ijk2];
654  clusters_ds[quanti_clusters]=clusters_lsample[ijk][ijk2]-clusters_ssample[ijk][ijk2];
655  clusters_mw[quanti_clusters]=(clusters_lwire[ijk][ijk2]+clusters_swire[ijk][ijk2])*0.5;
656  clusters_ms[quanti_clusters]=(clusters_lsample[ijk][ijk2]+clusters_ssample[ijk][ijk2])*0.5;
657  clusters_mintime[quanti_clusters]=(clusters_ssample[ijk][ijk2]);
658  quanti_clusters+=1;
659  }
660  }
661  }
662 
663 
664  for(int icl = 0; icl < quanti_clusters; ++icl){
665  if (clusters_qq[icl]>0) {
666  //std::cout << " CLUSTER INFO " << icl << " " << clusters_qq[icl] << " " << clusters_vi[icl] << " " << clusters_dw[icl] << " " << clusters_ds[icl] << " " << clusters_nn[icl] << std::endl;
667  int tpc_number=clusters_vi[icl];//qui andrebbe messo il numero di TPC
668 
669  if (clusters_ds[icl]>2250 && clusters_dw[icl]>100)
670  {//if analisi
671 
672  std::vector<float> *whc=new std::vector<float>;
673  std::vector<float> *shc=new std::vector<float>;
674  std::vector<float> *ahc=new std::vector<float>;
675  std::vector<float> *fahc=new std::vector<float>;
676 
677  if (clusters_vi[icl]==0) {
678  for (unsigned int ijk=0; ijk<www0->size(); ijk++) {
679  if ((*ccc0)[ijk]==clusters_nn[icl]) {
680  whc->push_back((*www0)[ijk]);
681  shc->push_back((*sss0)[ijk]);
682  ///ahc->push_back((*hhh0)[ijk]);
683  ahc->push_back((*aaa0)[ijk]);
684  fahc->push_back((*hhh0)[ijk]);
685  }
686  }
687  }
688  if (clusters_vi[icl]==2) {
689  for (unsigned int ijk=0; ijk<www2->size(); ijk++) {
690  if ((*ccc2)[ijk]==clusters_nn[icl]) {
691  whc->push_back((*www2)[ijk]);
692  shc->push_back((*sss2)[ijk]);
693  ahc->push_back((*aaa2)[ijk]);
694  //ahc->push_back((*hhh2)[ijk]);
695  fahc->push_back((*hhh2)[ijk]);
696  }
697  }
698  }
699 
700 
701  ////std::cout << " CLUSTER INFO " << icl << " " << clusters_qq[icl] << " " << whc->size() << std::endl;
702 
703 
704  if(whc->size()>100)//prima 0
705  {
706  float pendenza=0;float intercetta=0;int found_ok=0;
707  std::vector<int> *escluse=new std::vector<int>;
708 /*
709 for(int j=0;j<(int)whc->size();j++)
710 {
711 if(((*shc)[j]-clusters_mintime[icl]<200) || ((*shc)[j]-clusters_mintime[icl])>2150)escluse->push_back(j);
712 }
713 */
714  for(int j=0;j<(int)whc->size();j++)
715  {
716  if(found_ok<1)
717  {
718  Double_t wires[10000];
719  Double_t samples[10000];
720  Double_t ex[10000];
721  Double_t quale[10000];
722  Double_t ey[10000];
723  int quanti=0;
724  for(int k=0;k<(int)whc->size();k++)
725  {
726  if(Nothere(escluse,k)==3)
727  {
728  wires[quanti]=(*whc)[k]*3;
729  samples[quanti]=(*shc)[k]*0.628;
730  ex[quanti]=0;
731  ey[quanti]=0;
732  quale[quanti]=k;
733  quanti+=1;
734  }
735  }
736  ////////std::cout << quanti << std::endl;
737  TGraphErrors *gr3 = new TGraphErrors(quanti,wires,samples,ex,ey);
738  gr3->Fit("pol1","Q");
739  TF1 *fitfunc = gr3->GetFunction("pol1");
740  pendenza=fitfunc->GetParameter(1);
741  intercetta=fitfunc->GetParameter(0);
742  float distance_maximal=fdisfcl;
743  int quella_a_distance_maximal=0;
744  int found_max=0;
745  for(int jj=0;jj<quanti;jj++)
746  {
747  if((abs((pendenza)*(wires[jj])-samples[jj]+intercetta)/sqrt((pendenza)*pendenza+1))>distance_maximal)
748  {
749  found_max=1;
750  quella_a_distance_maximal=quale[jj];
751  distance_maximal=(abs(pendenza*wires[jj]-samples[jj]+intercetta)/sqrt(pendenza*pendenza+1));
752  }
753  }
754  if(found_max==1)escluse->push_back(quella_a_distance_maximal);
755  if(found_max==0)found_ok=1;
756  }
757  }
758  ////std::cout << escluse->size() << " escluse " << whc->size() << " " << found_ok << std::endl;
759  delete escluse;
760  std::vector<float> *hittime=new std::vector<float>;
761  std::vector<float> *hitwire=new std::vector<float>;
762  std::vector<float> *hitarea=new std::vector<float>;
763  std::vector<float> *hittimegood=new std::vector<float>;
764  std::vector<float> *hitareagood=new std::vector<float>;
765  std::vector<float> *hitwiregood=new std::vector<float>;
766  ///////std::cout << found_ok << std::endl;
767  ///////std::cout << pendenza << std::endl;
768  ///////std::cout << intercetta << std::endl;
769  if(found_ok==1)
770  {
771  for(int kkk=0;kkk<(int)whc->size();kkk++)
772  {
773  if((*ahc)[kkk]>0)// && (((*shc)[kkk]-clusters_mintime[icl]<200) || ((*shc)[kkk]-clusters_mintime[icl])<2150))
774  {
775  float distance=(abs(pendenza*((*whc)[kkk]*3)-(*shc)[kkk]*0.628+intercetta))/sqrt(pendenza*pendenza+1);
776  if(distance<=fdisfcl)
777  {
778  //cout << log((*ahc)[kkk]/(0.4*peach)) << endl;
779  hittime->push_back((*shc)[kkk]*0.4);
780  hitwire->push_back((*whc)[kkk]);
781  hitarea->push_back((*ahc)[kkk]);
782  h_hit_height->Fill((*fahc)[kkk]);
783  h_hit_area->Fill((*ahc)[kkk]);
784  h_hit_height_area->Fill((*fahc)[kkk],(*ahc)[kkk]);
785 
786  }
787  }
788  }
789  //float result_rms=0.14;
790  ////std::cout << result_rms << endl;
791  Double_t area[10000];
792  Double_t nologarea[10000];
793  Double_t tempo[10000];
794  Double_t ex[10000];
795  //Double_t quale[10000];
796  Double_t ey[10000];
797  Double_t ek[10000];
798  Double_t ez[10000];
799  ////std::cout << hitarea->size() << " dimensione hitarea" << std::endl;
800 
801  ////std::cout<<""<<std::endl;
802  ////std::cout<<"HERE line 802"<<std::endl;
803  ////std::cout<<""<<std::endl;
804  //std::ofstream purh("dump_purity_hits.out",std::ios::app);
805 
806  if(hitarea->size()>100)//prima 30
807  {
808  h_ratio->Fill(((float)whc->size())/((float)clusters_dw[icl]));
809  h_ratio_3->Fill(clusters_ds[icl],((float)whc->size())/((float)clusters_dw[icl]));
810 
811  ////std::cout << "RATIO INFO 1 " << (float)hitarea->size() << " " << (float)clusters_dw[icl] << " " << (((float)hitarea->size())/((float)clusters_dw[icl])) << std::endl;
812  float minimo=100000;
813  float massimo=0;
814  float wire_minimo=100000;
815  float wire_massimo=0;
816 
817  float delta_sample_selected=0;
818  float delta_wire_selected=0;
819  float wire_del_massimo=-1;
820  float wire_del_minimo=-1;
821  float sample_massimo=-1;
822  float sample_minimo=-1;
823  int quante_hit_nel_range_tempo=0;
824  purh<< evt.run() << " " << evt.subRun() << " " << evt.event() << " -1 " << " -1 " << " -1 " << " -1 " << std::endl;
825  for(int kk=0;kk<(int)hitarea->size();kk++)
826  {
827  if(fdumphitsfcl==1)purh<< evt.run() << " SELECTED " << evt.subRun() << " " << evt.event() << " " << tpc_number << " " << (*hitwire)[kk] << " " << (*hittime)[kk] << " " << (*hitarea)[kk] << std::endl;
828 
829 quante_hit_nel_range_tempo+=1;
830  if((*hittime)[kk]>massimo)
831  {
832  massimo=(*hittime)[kk];
833  wire_del_massimo=(*hitwire)[kk];
834  }
835 
836 
837  if((*hittime)[kk]<minimo)
838  {
839  minimo=(*hittime)[kk];
840  wire_del_minimo=(*hitwire)[kk];
841  }
842 
843 
844  if((*hitwire)[kk]>wire_massimo)
845  {
846  wire_massimo=(*hitwire)[kk];
847  }
848 
849 
850  if((*hitwire)[kk]<wire_minimo)
851  {
852  wire_minimo=(*hitwire)[kk];
853  }
854 
855  }
856 
857  delta_sample_selected=(massimo-minimo)/0.4;
858  delta_wire_selected=wire_del_massimo-wire_del_minimo;
859  sample_massimo=massimo/0.4;
860  sample_minimo=minimo/0.4;
861  h_ratio_after->Fill(((float)quante_hit_nel_range_tempo)/fabs(delta_wire_selected));
862  h_ratio_after_2->Fill(((float)quante_hit_nel_range_tempo)/fabs(wire_massimo-wire_minimo+1));
863  h_ratio_after_3->Fill(delta_sample_selected,((float)quante_hit_nel_range_tempo)/fabs(wire_massimo-wire_minimo+1));
864 
865  ////std::cout << "RATIO INFO 2 " << (float)hitarea->size() << " " << fabs(delta_wire_selected) << " " << (((float)hitarea->size())/fabs(delta_wire_selected)) << std::endl;
866 
867  ////std::cout << hitarea->size() << std::endl;
868  //int gruppi=hitarea->size()/50;
869  int gruppi=fgruppifcl;//originale 8
870  ////std::cout << gruppi << std::endl;
871 
872  float steptime=(massimo-minimo)/(gruppi+1);
873  ////std::cout << steptime << " steptime " << minimo << " " << massimo << std::endl;
874  float starting_value_tau=fValoretaufcl;
875 
876  ////std::cout << starting_value_tau << " VALORE INDICATIVO TAU " << std::endl;
877  //if(tpc_number==2 || tpc_number==5)starting_value_tau=6500;
878  //if(tpc_number==10 || tpc_number==13)starting_value_tau=5700;
879  for(int stp=0;stp<=gruppi;stp++)
880  {
881  std::vector<float>* hitpertaglio=new std::vector<float>;
882  ////std::cout << 500+stp*steptime << " time " << 500+(stp+1)*(steptime) << std::endl;
883  ///////////std::cout << minimo+stp*steptime << " " << minimo+(stp+1)*(steptime) << std::endl;
884  for(int kk=0;kk<(int)hitarea->size();kk++)
885  {
886  if((*hittime)[kk]>=(minimo+stp*steptime) && (*hittime)[kk]<=(minimo+(stp+1)*(steptime)))
887  hitpertaglio->push_back((*hitarea)[kk]*exp((*hittime)[kk]/starting_value_tau));
888  }
889  ///////////std::cout << hitpertaglio->size() << std::endl;
890  float tagliomax=FoundMeanLog(hitpertaglio,fmaxfcl);//0.9
891  float tagliomin=FoundMeanLog(hitpertaglio,fminfcl);//0.05
892  //float tagliomin=0;
893  //float tagliomax=1000000;
894  delete hitpertaglio;
895  ////std::cout << tagliomax << " t " << std::endl;
896  for(int kk=0;kk<(int)hitarea->size();kk++)
897  {
898  ////std::cout << (*hittime)[kk] << " " << (*hitwire)[kk] << " " << (*hitarea)[kk] << " " << (minimo+stp*steptime) << " " << (minimo+(stp+1)*steptime) << " " << (*hitarea)[kk]*exp((*hittime)[kk]/starting_value_tau) << std::endl;
899  if((*hittime)[kk]>(minimo+stp*steptime) &&
900  (*hittime)[kk]<(minimo+(stp+1)*steptime) &&
901  (*hitarea)[kk]*exp((*hittime)[kk]/starting_value_tau)<tagliomax &&
902  (*hitarea)[kk]*exp((*hittime)[kk]/starting_value_tau)>tagliomin)
903  {
904  ////std::cout << ((*hitarea)[kk]*exp((*hittime)[kk]/1400)) << " GOOD " << (*hitarea)[kk] << " " << (*hittime)[kk] << std::endl;
905  hitareagood->push_back((*hitarea)[kk]);
906  hittimegood->push_back((*hittime)[kk]);
907  hitwiregood->push_back((*hitwire)[kk]);
908  }
909  }
910  }
911  ////std::cout << hitareagood->size() << " hitareagood" << std::endl;
912 if(delta_sample_selected>1900)
913 {
914  for(int k=0;k<(int)whc->size();k++)
915  {
916  h_hittime->Fill((*shc)[k]-clusters_mintime[icl]);
917  }
918  for(int k=0;k<(int)hittime->size();k++)
919  {
920  h_hittime_2->Fill((*hittime)[k]/0.4-clusters_mintime[icl]);
921  }
922 
923  for(int k=0;k<(int)hittimegood->size();k++)
924  {
925  h_hittime_3->Fill((*hittimegood)[k]/0.4-clusters_mintime[icl]);
926  }
927 }
928 
929  for(int k=0;k<(int)hitareagood->size();k++)
930  {
931  if(fdumphitsfcl==1)purh<< evt.run() << " AFTER CLEAN " << evt.subRun() << " " << evt.event() << " " << tpc_number << " " << (*hitwiregood)[k] << " " << (*hittimegood)[k] << " " << (*hitareagood)[k] << std::endl;
932  //if((*hittimegood)[k]-600*0.4<=1000)//correzione 15/08
933  tempo[k]=(*hittimegood)[k];
934  area[k]=log((*hitareagood)[k]);
935  ////std::cout << (*hitareagood)[k] << " " << area[k] << std::endl;
936  nologarea[k]=((*hitareagood)[k]);
937  ex[k]=0;
938  ez[k]=60;
939  ey[k]=0.23;
940  }
941  ////std::cout<<""<<std::endl;
942  ////std::cout<<"HERE line 872"<<std::endl;
943  ////std::cout<<""<<std::endl;
944  TGraphErrors *gr31 = new TGraphErrors(hitareagood->size(),tempo,area,ex,ey);
945  //TGraphErrors *gr4 = new TGraphErrors(hitareagood->size(),tempo,nologarea,ex,ey);
946  gr31->Fit("pol1","Q");
947  TF1 *fit = gr31->GetFunction("pol1");
948  float slope_purity=fit->GetParameter(1);
949  //float error_slope_purity=fit->GetParError(1);
950  float intercetta_purezza=fit->GetParameter(0);
951 
952  TH1F *h111 = new TH1F("h111","delta aree",200,-10,10);
953  float sum_per_rms_test=0;
954  int quanti_in_h111=0;
955  for(int k=0;k<(int)hitareagood->size();k++)
956  {
957  h111->Fill(area[k]-slope_purity*tempo[k]-intercetta_purezza);
958  sum_per_rms_test+=(area[k]-slope_purity*tempo[k]-intercetta_purezza)*(area[k]-slope_purity*tempo[k]-intercetta_purezza);
959  if((area[k]-slope_purity*tempo[k]-intercetta_purezza)>-10 && (area[k]-slope_purity*tempo[k]-intercetta_purezza)<10)quanti_in_h111+=1;
960  }
961 
962  float error=100.;
963  if(quanti_in_h111>1){
964  h111->Fit("gaus","Q");
965  TF1 *fitg = h111->GetFunction("gaus");
966  error=fitg->GetParameter(2);
967  }
968  ////std::cout << " error " << error << std::endl;
969  //float error_2=sqrt(sum_per_rms_test/(hitareagood->size()-2));
970  ////std::cout << " error vero" << error_2 << std::endl;
971  h111->Delete();
972 
973 
974 
975  TGraphErrors *gr4 = new TGraphErrors(hitareagood->size(),tempo,nologarea,ex,ey);
976  gr4->Fit("expo","Q");
977  TF1 *fite = gr4->GetFunction("expo");
978  slope_purity=fite->GetParameter(1);
979  intercetta_purezza=fite->GetParameter(0);
980  float mean_hit_area=0;
981  float size_hit_area=hitareagood->size();
982  int quanti_in_h111e=0;
983  TH1F *h111e = new TH1F("h111e","delta aree",100,-1000.,1000.);
984  for(int k=0;k<(int)hitareagood->size();k++)
985  {
986  h111e->Fill(nologarea[k]-exp(slope_purity*tempo[k]+intercetta_purezza));
987  mean_hit_area+=nologarea[k]/size_hit_area;
988  if((nologarea[k]-exp(slope_purity*tempo[k]+intercetta_purezza))>-1000. && (nologarea[k]-exp(slope_purity*tempo[k]+intercetta_purezza))<1000)quanti_in_h111e+=1;
989  //cout << nologarea[k]-exp(slope_purity*tempo[k]+intercetta_purezza) << endl;
990  }
991 
992  float error_expo=1000.;
993  if(quanti_in_h111e>1)
994  {
995  h111e->Fit("gaus","Q");
996  TF1 *fitge = h111e->GetFunction("gaus");
997  error_expo=fitge->GetParameter(2);
998  }
999  //std::cout << " errors " << error << " " << error_expo << std::endl;
1000  h_errors->Fill(error_expo);
1001  h111e->Delete();//fitge->Delete();fite->Delete();
1002 
1003  for(int k=0;k<(int)hitareagood->size();k++)
1004  {
1005  ek[k]=error_expo;
1006  ez[k]=error_expo;
1007  ey[k]=error;
1008  }
1009  ////std::cout<<""<<std::endl;
1010  ////std::cout<<"HERE line 906"<<std::endl;
1011  ////std::cout<<""<<std::endl;
1012 
1013  TGraphErrors *gr32 = new TGraphErrors(hitareagood->size(),tempo,area,ex,ey);
1014  gr32->Fit("pol1","Q");
1015 
1016  TF1 *fit2 = gr32->GetFunction("pol1");
1017  float slope_purity_2=fit2->GetParameter(1);
1018  float error_slope_purity_2=fit2->GetParError(1);
1019  //float intercetta_purezza_2=fit2->GetParameter(0);
1020  float chiquadro=fit2->GetChisquare()/(hitareagood->size()-2);
1021  std::ofstream goodpuro("purity_results.out",std::ios::app);
1022  std::ofstream goodpuro2("purity_results2.out",std::ios::app);
1023 
1024  ////std::cout << -1/slope_purity_2 << std::endl;
1025  ////std::cout << -1/(slope_purity_2+error_slope_purity_2)+1/slope_purity_2 << std::endl;
1026  ////std::cout << 1/slope_purity_2-1/(slope_purity_2-error_slope_purity_2) << std::endl;
1027  TGraphAsymmErrors *gr41 = new TGraphAsymmErrors (hitareagood->size(),tempo,nologarea,ex,ex,ez,ek);
1028  gr41->Fit("expo","Q");
1029  TF1 *fitexo = gr41->GetFunction("expo");
1030  float slope_purity_exo=fitexo->GetParameter(1);
1031  float error_slope_purity_exo=fitexo->GetParError(1);
1032  //fRunSubPurity2->Fill(evt.run(),evt.subRun(),-slope_purity_exo*1000.);
1033  //fRunSubPurity->Fill(evt.run(),evt.subRun(),-slope_purity_2*1000.);
1034  ////std::cout << -1/slope_purity_exo << std::endl;
1035  ////std::cout << -1/(slope_purity_exo+error_slope_purity_exo)+1/slope_purity_exo << std::endl;
1036  ////std::cout << 1/slope_purity_exo-1/(slope_purity_exo-error_slope_purity_exo) << std::endl;
1037  ////std::cout << fitexo->GetChisquare()/(hitareagood->size()-2) << std::endl;
1038 
1039 
1040  if(fabs(slope_purity_2)<0.01 || fabs(slope_purity_exo)<0.01)
1041  {
1042  if(fabs(slope_purity_2)<0.01)purityvalues->Fill(-slope_purity_2*1000.);
1043  if(fabs(slope_purity_2)<0.01)goodpuro << evt.run() << " " << evt.subRun() << " " << evt.event() << " " << tpc_number << " " << slope_purity_2 << " " << error_slope_purity_2 << " " << chiquadro << " " << clusters_dw[icl] << " " << clusters_ds[icl] << std::endl;
1044 
1045  if(fabs(slope_purity_exo)<0.01)goodpuro2<< evt.run() << " " << evt.subRun() << " " << evt.event() << " " << tpc_number+fcryofcl*10 << " " << slope_purity_exo << " " << error_slope_purity_exo << " " << fitexo->GetChisquare()/(hitareagood->size()-2) << " " << clusters_dw[icl] << " " << clusters_ds[icl] << " " << clusters_mw[icl] << " " << clusters_ms[icl] << " " << delta_wire_selected << " " << delta_sample_selected << " " << sample_minimo << " " << sample_massimo << " " << wire_del_minimo << " " << wire_del_massimo << " " << wire_minimo << " " << wire_massimo << " " << whc->size() << " " << hitarea->size() << " " << quante_hit_nel_range_tempo << " " << error_expo << " " << clusters_mintime[icl] << " " << mean_hit_area << std::endl;
1046 
1047 
1048  if(fabs(slope_purity_exo)<0.01)
1049  {
1050  run_tree=evt.run();
1051  subrun_tree=evt.subRun();
1052  event_tree=evt.event();
1053  tpc_tree=tpc_number+fcryofcl*10;
1054  slope_tree=-slope_purity_exo*1000;
1055  errslope_tree=error_slope_purity_exo*1000;
1056  chi_tree=fitexo->GetChisquare()/(hitareagood->size()-2);
1057  dtime_tree=delta_sample_selected;
1058  dwire_tree=wire_massimo-wire_minimo+1;
1059  earea_tree=error_expo;
1060  marea_tree=mean_hit_area;
1061  qhits_tree=hitarea->size();
1062  fpurTree->Fill();
1063  }
1064  if(fabs(slope_purity_exo)<0.01)purityvalues2->Fill(-slope_purity_exo*1000.);
1065 
1066  if(fabs(slope_purity_exo)<0.01 && tpc_number==0)puritytpc0->Fill(-slope_purity_exo*1000.);
1067  if(fabs(slope_purity_exo)<0.01 && tpc_number==1)puritytpc1->Fill(-slope_purity_exo*1000.);
1068  if(fabs(slope_purity_exo)<0.01 && tpc_number==2)puritytpc2->Fill(-slope_purity_exo*1000.);
1069  if(fabs(slope_purity_exo)<0.01 && tpc_number==3)puritytpc3->Fill(-slope_purity_exo*1000.);
1070 
1071  }
1072  if(fabs(slope_purity_exo)<0.01)
1073  {
1074 
1075 
1076  anab::TPCPurityInfo purity_info;
1077  purity_info.Run = evt.run();
1078  purity_info.Subrun = evt.subRun();
1079  purity_info.Event = evt.event();
1080  purity_info.TPC = tpc_number;
1081  purity_info.Wires = delta_wire_selected;
1082  purity_info.Ticks = delta_sample_selected;
1083 
1084  //if(purity_info.TPC<2) purity_info.Cryostat=0;
1085  //else purity_info.Cryostat=1;
1086 
1087  purity_info.Cryostat=fcryofcl;
1088 
1089  // near/far from cathode tracks
1090 /*
1091  if(delta_wire_selected< 0){
1092  purity_info.AttenuationNEAR = slope_purity_exo*-1000.;
1093  }
1094  else purity_info.AttenuationNEAR = 0;
1095  if(delta_wire_selected>0){
1096  purity_info.AttenuationFAR = slope_purity_exo*-1000.;
1097  }
1098  else purity_info.AttenuationFAR = 0;
1099 
1100 */
1101 
1102  purity_info.Attenuation = slope_purity_exo*-1.;
1103  purity_info.FracError = error_slope_purity_exo / slope_purity_exo;
1104  //purity_info.Attenuation_2 = slope_purity_2*-1.;
1105  //purity_info.FracError_2 = error_slope_purity_2 / slope_purity_2;
1106 
1107  //purity_info.Wires = clusters_dw[icl];
1108  // purity_info.Ticks = clusters_ds[icl];
1109 
1110  if(fFillAnaTuple)
1111  purityTuple->Fill(purity_info.Run,purity_info.Event,purity_info.TPC,purity_info.Wires,purity_info.Ticks,purity_info.Attenuation);
1112 
1113  ////std::cout << "Calling again after filling attenuation … " << std::endl;
1114  purity_info.Print();
1115  outputPtrVector->push_back(purity_info);
1116  }
1117  ////std::cout << ts << " is time event " << std::endl;
1118  //goodpur << -1/slope_purity_exo << std::endl;
1119  //goodpur << -1/(slope_purity_exo+error_slope_purity_exo)+1/slope_purity_exo << std::endl;
1120  //goodpur << 1/slope_purity_exo-1/(slope_purity_exo-error_slope_purity_exo) << std::endl;
1121  //goodpur << timeevent << " is time event " << std::endl;
1122  }
1123  }
1124  ////std::cout << "Delete hit stuff." << std::endl;
1125 
1126  delete hittime;
1127  delete hitarea;
1128  delete hittimegood;
1129  delete hitareagood;
1130  delete hitwiregood;
1131  delete hitwire;
1132  }
1133 
1134  ////std::cout << "Delete cluster stuff." << std::endl;
1135 
1136  delete shc;
1137  delete ahc;
1138  delete fahc;
1139  delete whc;
1140 
1141  }//fine if ananlisi
1142  }
1143  }
1144 
1145  ////std::cout << "Delete big stuff." << std::endl;
1146 
1147  delete www0;
1148  delete sss0;
1149  delete hhh0;
1150  delete ehh0;
1151  delete ccc0;
1152  delete www2;
1153  delete sss2;
1154  delete hhh2;
1155  delete ehh2;
1156  delete ccc2;
1157  delete aaa0;
1158  delete aaa2;
1159 
1160  }
1161 
1162  //std::cout << "Checking everything in the output..." << std::endl;
1163  //std::cout << "There are " << outputPtrVector->size() << " objects in the output vector." << std::endl;
1164  /* //don't need this printed here.
1165  for (size_t i_info = 0; i_info<outputPtrVector->size(); ++i_info){
1166  auto info = outputPtrVector->at(i_info);
1167  info.Print();
1168  }
1169  */
1170 
1171  //put info onto the event
1172  evt.put(std::move(outputPtrVector));
1173  } // produces
1174 
1175 } //end namespace
1176 
1177 namespace icarus{
1178 
1179  DEFINE_ART_MODULE(ICARUSPurityDQM)
1180 
1181 }
1182 
pmtana::std
double std(const std::vector< short > &wf, const double ped_mean, size_t start, size_t nsample)
Definition: UtilFunc.cxx:42
raw::RawDigit
Collection of charge vs time digitized from a single readout channel.
Definition: RawDigit.h:69
test_msg.error
error
Definition: larcorealg/test/GeoAlgo/test_msg.py:32
icarus::ICARUSPurityDQM::produce
void produce(art::Event &evt)
read access to event
Definition: ICARUSPurityDQM_module.cc:344
Hit.h
Declaration of signal hit object.
icarus::ICARUSPurityDQM::Nothere
int Nothere(std::vector< int > *a, int b)
Definition: ICARUSPurityDQM_module.cc:275
right
walls no right
Definition: selectors.fcl:105
get
get()
Definition: PMTconfigurationExtraction_module.cc:168
geo::PlaneID::PlaneID_t
unsigned int PlaneID_t
Type for the ID number.
Definition: geo_types.h:473
raw::RawDigit::Channel
ChannelID_t Channel() const
DAQ channel this raw data was read from.
Definition: RawDigit.h:212
geo::CryostatID::Cryostat
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:212
larg4::bool
static constexpr bool
Definition: OpFastScintillation.hh:588
RawDigit.h
Definition of basic raw digits.
anab::TPCPurityInfo::Print
void Print() const
Definition: TPCPurityInfo.cc:15
geo::WireID::Wire
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
icarus::ICARUSPurityDQM::Notheref
int Notheref(std::vector< float > *a, float b)
Definition: ICARUSPurityDQM_module.cc:288
a
process_name gaushit a
Definition: decode_signalprocess_icarus.fcl:37
anab::TPCPurityInfo::Subrun
unsigned int Subrun
Definition: TPCPurityInfo.hh:19
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
anab::TPCPurityInfo::Event
unsigned int Event
Definition: TPCPurityInfo.hh:20
icarus::distance
double distance(geo::Point_t const &point, CathodeDesc_t const &cathode)
Returns the distance of a point from the cathode.
Definition: CathodeCrossingUtils.cxx:35
raw.h
Collect all the RawData header files together.
anab::TPCPurityInfo::Cryostat
unsigned int Cryostat
Definition: TPCPurityInfo.hh:22
left
walls no left
Definition: selectors.fcl:105
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
Cluster.h
Declaration of cluster object.
make_calib_df.exp
def exp
Definition: make_calib_df.py:43
icarus::ICARUSPurityDQM::FoundMeanLog
Double_t FoundMeanLog(std::vector< float > *a, float b)
Definition: ICARUSPurityDQM_module.cc:301
Wire.h
Declaration of basic channel signal object.
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
k
pdgs k
Definition: selectors.fcl:22
raw::ChannelID_t
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
geo::TPCID::TPC
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
raw::Uncompress
void Uncompress(const std::vector< short > &adc, std::vector< short > &uncompressed, raw::Compress_t compress)
Uncompresses a raw data buffer.
Definition: raw.cxx:776
icarus::ICARUSPurityDQM::ICARUSPurityDQM
ICARUSPurityDQM(fhicl::ParameterSet const &pset)
Definition: ICARUSPurityDQM_module.cc:175
icarus::ICARUSPurityDQM::fDigitModuleLabel
std::vector< art::InputTag > fDigitModuleLabel
Definition: ICARUSPurityDQM_module.cc:130
Geometry.h
art framework interface to geometry description
Generated by   doxygen 1.8.5