Adding the AliAnalysisGUI class which is the main class that controls the GUI.
[u/mrichter/AliRoot.git] / ITS / AliITSVertexerZ.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15 #include "AliITSVertexerZ.h"
16 #include<TBranch.h>
17 #include<TH1.h>
18 #include <TString.h>
19 #include<TTree.h>
20 #include "AliITSLoader.h"
21 #include "AliITSgeom.h"
22 #include "AliITSDetTypeRec.h"
23 #include "AliITSRecPoint.h"
24
25 /////////////////////////////////////////////////////////////////
26 // this class implements a fast method to determine
27 // the Z coordinate of the primary vertex
28 // for p-p collisions it seems to give comparable or better results
29 // with respect to what obtained with AliITSVertexerPPZ
30 // It can be used successfully with Pb-Pb collisions
31 ////////////////////////////////////////////////////////////////
32
33 ClassImp(AliITSVertexerZ)
34
35
36
37 //______________________________________________________________________
38 AliITSVertexerZ::AliITSVertexerZ():AliITSVertexer(),
39 fFirstL1(0),
40 fLastL1(0),
41 fFirstL2(0),
42 fLastL2(0),
43 fDiffPhiMax(0),
44 fX0(0.),
45 fY0(0.),
46 fZFound(0),
47 fZsig(0.),
48 fZCombc(0),
49 fZCombv(0),
50 fZCombf(0),
51 fLowLim(0.),
52 fHighLim(0.),
53 fStepCoarse(0),
54 fStepFine(0),
55 fTolerance(0.),
56 fMaxIter(0){
57   // Default constructor
58   SetDiffPhiMax(0);
59   SetFirstLayerModules(0);
60   SetSecondLayerModules(0);
61   SetLowLimit(0.);
62   SetHighLimit(0.);
63   SetBinWidthCoarse(0.);
64   SetBinWidthFine(0.);
65   SetTolerance(0.);
66   SetPPsetting(0.,0.);
67   ConfigIterations();
68 }
69
70 //______________________________________________________________________
71 AliITSVertexerZ::AliITSVertexerZ(TString fn, Float_t x0, Float_t y0):AliITSVertexer(fn),
72 fFirstL1(0),
73 fLastL1(0),
74 fFirstL2(0),
75 fLastL2(0),
76 fDiffPhiMax(0),
77 fX0(x0),
78 fY0(y0),
79 fZFound(0),
80 fZsig(0.),
81 fZCombc(0),
82 fZCombv(0),
83 fZCombf(0),
84 fLowLim(0.),
85 fHighLim(0.),
86 fStepCoarse(0),
87 fStepFine(0),
88 fTolerance(0.),
89 fMaxIter(0) {
90   // Standard Constructor
91   SetDiffPhiMax();
92   SetFirstLayerModules();
93   SetSecondLayerModules();
94   SetLowLimit();
95   SetHighLimit();
96   SetBinWidthCoarse();
97   SetBinWidthFine();
98   SetTolerance();
99   SetPPsetting();
100   ConfigIterations();
101
102 }
103
104 //______________________________________________________________________
105 AliITSVertexerZ::AliITSVertexerZ(const AliITSVertexerZ &vtxr) : AliITSVertexer(vtxr),
106 fFirstL1(vtxr.fFirstL1),
107 fLastL1(vtxr.fLastL1),
108 fFirstL2(vtxr.fFirstL2),
109 fLastL2(vtxr.fLastL2),
110 fDiffPhiMax(vtxr.fDiffPhiMax),
111 fX0(vtxr.fX0),
112 fY0(vtxr.fY0),
113 fZFound(vtxr.fZFound),
114 fZsig(vtxr.fZsig),
115 fZCombc(vtxr.fZCombc),
116 fZCombv(vtxr.fZCombv),
117 fZCombf(vtxr.fZCombf),
118 fLowLim(vtxr.fLowLim),
119 fHighLim(vtxr.fHighLim),
120 fStepCoarse(vtxr.fStepCoarse),
121 fStepFine(vtxr.fStepFine),
122 fTolerance(vtxr.fTolerance),
123 fMaxIter(vtxr.fMaxIter) {
124   // Copy constructor
125
126 }
127
128 //______________________________________________________________________
129 AliITSVertexerZ& AliITSVertexerZ::operator=(const AliITSVertexerZ&  vtxr ){
130   // Assignment operator
131
132   this->~AliITSVertexerZ();
133   new(this) AliITSVertexerZ(vtxr);
134   return *this;
135 }
136
137 //______________________________________________________________________
138 AliITSVertexerZ::~AliITSVertexerZ() {
139   // Destructor
140   delete fZCombc;
141   delete fZCombf;
142   delete fZCombv;
143 }
144
145 //______________________________________________________________________
146 void AliITSVertexerZ::ConfigIterations(Int_t noiter,Float_t *ptr){
147   // configure the iterative procedure to gain efficiency for
148   // pp events with very low multiplicity
149   Float_t defaults[5]={0.05,0.1,0.2,0.3,0.5};
150   fMaxIter=noiter;
151   if(noiter>5){
152     Error("ConfigIterations","Maximum number of iterations is 5\n");
153     fMaxIter=5;
154   }
155   for(Int_t j=0;j<5;j++)fPhiDiffIter[j]=defaults[j];
156   if(ptr)for(Int_t j=0;j<fMaxIter;j++)fPhiDiffIter[j]=ptr[j];
157 }
158
159 //______________________________________________________________________
160 AliESDVertex* AliITSVertexerZ::FindVertexForCurrentEvent(Int_t evnumber){
161   // Defines the AliESDVertex for the current event
162   VertexZFinder(evnumber);
163   Int_t ntrackl=0;
164   for(Int_t iteraz=0;iteraz<fMaxIter;iteraz++){
165     if(fCurrentVertex) ntrackl=fCurrentVertex->GetNContributors();
166     if(!fCurrentVertex || ntrackl==0 || ntrackl==-1){
167       Float_t diffPhiMaxOrig=fDiffPhiMax;
168       fDiffPhiMax=GetPhiMaxIter(iteraz);
169       VertexZFinder(evnumber);
170       fDiffPhiMax=diffPhiMaxOrig;
171     }
172   }
173   FindMultiplicity(evnumber);
174   return fCurrentVertex;
175 }  
176
177
178
179
180 //______________________________________________________________________
181 void AliITSVertexerZ::VertexZFinder(Int_t evnumber){
182   // Defines the AliESDVertex for the current event
183   fCurrentVertex = 0;
184   AliRunLoader *rl =AliRunLoader::GetRunLoader();
185   AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
186   AliITSgeom* geom = itsLoader->GetITSgeom();
187   itsLoader->LoadRecPoints();
188   rl->GetEvent(evnumber);
189
190   AliITSDetTypeRec detTypeRec;
191
192   TTree *tR = itsLoader->TreeR();
193   detTypeRec.SetTreeAddressR(tR);
194   TClonesArray *itsRec  = 0;
195   // lc and gc are local and global coordinates for layer 1
196   Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
197   Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
198   // lc2 and gc2 are local and global coordinates for layer 2
199   Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
200   Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
201
202   itsRec = detTypeRec.RecPoints();
203   TBranch *branch;
204   branch = tR->GetBranch("ITSRecPoints");
205
206   Int_t nrpL1 = 0;
207   Int_t nrpL2 = 0;
208   // By default fFirstL1=0 and fLastL1=79
209   // This loop counts the number of recpoints on layer1 (central modules)
210   for(Int_t module= fFirstL1; module<=fLastL1;module++){
211     // Keep only central modules
212     if(module%4==0 || module%4==3)continue;
213     //   cout<<"Procesing module "<<module<<" ";
214     branch->GetEvent(module);
215     //    cout<<"Number of clusters "<<clusters->GetEntries()<<endl;
216     nrpL1+= itsRec->GetEntries();
217     detTypeRec.ResetRecPoints();
218   }
219   //By default fFirstL2=80 and fLastL2=239
220   //This loop counts the number of RP on layer 2
221   for(Int_t module= fFirstL2; module<=fLastL2;module++){
222     branch->GetEvent(module);
223     nrpL2+= itsRec->GetEntries();
224     detTypeRec.ResetRecPoints();
225   }
226   if(nrpL1 == 0 || nrpL2 == 0){
227     ResetHistograms();
228     itsLoader->UnloadRecPoints();
229     fCurrentVertex = new AliESDVertex(0.,5.3,-2);
230     return;
231   }
232   // The vertex finding is attempted only if the number of RP is !=0 on
233   // both layers
234   Float_t *xc1 = new Float_t [nrpL1]; // coordinates of the L1 Recpoints
235   Float_t *yc1 = new Float_t [nrpL1];
236   Float_t *zc1 = new Float_t [nrpL1];
237   Float_t *phi1 = new Float_t [nrpL1];
238   Float_t *xc2 = new Float_t [nrpL2]; // coordinates of the L1 Recpoints
239   Float_t *yc2 = new Float_t [nrpL2];
240   Float_t *zc2 = new Float_t [nrpL2];
241   Float_t *phi2 = new Float_t [nrpL2];
242   Int_t ind = 0;// running index for RP
243   // Force a coarse bin size of 200 microns if the number of clusters on layer 2
244   // is low
245   if(nrpL2<fPPsetting[0])SetBinWidthCoarse(fPPsetting[1]);
246   // By default nbinfine = (10+10)/0.0005=40000
247   Int_t nbinfine = static_cast<Int_t>((fHighLim-fLowLim)/fStepFine);
248   // By default nbincoarse=(10+10)/0.01=2000
249   Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse);
250   // Set stepverycoarse = 3*fStepCoarse
251   Int_t nbinvcoarse = static_cast<Int_t>((fHighLim-fLowLim)/(3.*fStepCoarse));
252   if(fZCombc)delete fZCombc;
253   fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
254   if(fZCombv)delete fZCombv;
255   fZCombv = new TH1F("fZCombv","Z",nbinvcoarse,fLowLim,fLowLim+nbinvcoarse*3.*fStepCoarse);
256   if(fZCombf)delete fZCombf;
257   fZCombf = new TH1F("fZCombf","Z",nbinfine,fLowLim,fLowLim+nbinfine*fStepFine);
258   // Loop on modules of layer 1 (restricted to central modules)
259   for(Int_t module= fFirstL1; module<=fLastL1;module++){
260     if(module%4==0 || module%4==3)continue;
261     branch->GetEvent(module);
262     Int_t nrecp1 = itsRec->GetEntries();
263     for(Int_t j=0;j<nrecp1;j++){
264       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
265       // Local coordinates of this recpoint
266       lc[0]=recp->GetDetLocalX();
267       lc[2]=recp->GetDetLocalZ();
268       geom->LtoG(module,lc,gc);
269       // Global coordinates of this recpoints
270       gc[0]-=fX0; // Possible beam offset in the bending plane
271       gc[1]-=fY0; //   "               "
272       xc1[ind]=gc[0];
273       yc1[ind]=gc[1];
274       zc1[ind]=gc[2];
275       // azimuthal angle is computed in the interval 0 --> 2*pi
276       phi1[ind] = TMath::ATan2(gc[1],gc[0]);
277       if(phi1[ind]<0)phi1[ind]=2*TMath::Pi()+phi1[ind];
278       ind++;
279     }
280     detTypeRec.ResetRecPoints();
281   }
282   ind = 0; // the running index is reset for Layer 2
283   for(Int_t module= fFirstL2; module<=fLastL2;module++){
284     branch->GetEvent(module);
285     Int_t nrecp2 = itsRec->GetEntries();
286     for(Int_t j=0;j<nrecp2;j++){
287       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
288       lc[0]=recp->GetDetLocalX();
289       lc[2]=recp->GetDetLocalZ();
290       geom->LtoG(module,lc,gc);
291       gc[0]-=fX0;
292       gc[1]-=fY0;
293       xc2[ind]=gc[0];
294       yc2[ind]=gc[1];
295       zc2[ind]=gc[2];
296       phi2[ind] = TMath::ATan2(gc[1],gc[0]);
297       if(phi2[ind]<0)phi2[ind]=2*TMath::Pi()+phi2[ind];
298       ind++;
299     }
300     detTypeRec.ResetRecPoints();
301   }
302  
303   // Int_t nolines=0;
304   for(Int_t i=0;i<nrpL1;i++){ // loop on L1 RP
305     Float_t r1=TMath::Sqrt(xc1[i]*xc1[i]+yc1[i]*yc1[i]); // radius L1 RP
306     for(Int_t j=0;j<nrpL2;j++){ // loop on L2 RP
307       Float_t diff = TMath::Abs(phi2[j]-phi1[i]); // diff in azimuth
308       if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff; //diff<pi
309       if(diff<fDiffPhiMax){ // cut on 10 milliradians by def.
310         Float_t r2=TMath::Sqrt(xc2[j]*xc2[j]+yc2[j]*yc2[j]); // radius L2 RP
311         Float_t zr0=(r2*zc1[i]-r1*zc2[j])/(r2-r1); //Z @ null radius
312         fZCombf->Fill(zr0);
313         fZCombc->Fill(zr0);
314         fZCombv->Fill(zr0);
315         Double_t pA[3];
316         Double_t pB[3];
317         pA[0]=xc1[i];
318         pA[1]=yc1[i];
319         pA[2]=zc1[i];
320         pB[0]=xc2[j];
321         pB[1]=yc2[j];
322         pB[2]=zc2[j];
323         //      MakeTracklet(pA,pB,nolines);
324       }
325     }
326   }
327   delete [] xc1;
328   delete [] yc1;
329   delete [] zc1;
330   delete [] phi1;
331   delete [] xc2;
332   delete [] yc2;
333   delete [] zc2;
334   delete [] phi2;
335   Double_t contents = fZCombc->GetEntries()- fZCombc->GetBinContent(0)-fZCombc->GetBinContent(nbincoarse+1);
336   if(contents<1.){
337     //    Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
338     ResetHistograms();
339     itsLoader->UnloadRecPoints();
340     fCurrentVertex = new AliESDVertex(0.,5.3,-1);
341     return;
342   }
343   /*
344   else {
345     if(fDebug>0)cout<<"Number of entries in hist. "<<fZCombc->GetEntries()<<endl;
346   }
347   */
348
349   TH1F *hc = fZCombc;
350   Bool_t goon = kFALSE;
351   Int_t cnt = 0;
352   Int_t bi;
353
354   do {
355     goon = kFALSE;
356     cnt++;
357     bi = hc->GetMaximumBin();   // bin with maximal content on coarse histogram
358     if(hc->GetBinContent(bi)<3){
359       if(cnt==1)goon = kTRUE;
360       hc = fZCombv;
361     }
362   } while(goon);
363
364
365   Float_t centre = hc->GetBinCenter(bi);  // z value of the bin with maximal content
366   
367   Int_t bifine=static_cast<Int_t>((hc->GetBinCenter(bi)-fZCombf->GetBinLowEdge(0))/fStepFine);
368   Int_t nbinsfine=static_cast<Int_t>(3*hc->GetBinWidth(bi)/fStepFine);
369   // evaluation of the centroid
370   Int_t ii1=TMath::Max(bifine-nbinsfine,1);
371   Int_t ii2=TMath::Min(bifine+nbinsfine,fZCombf->GetNbinsX());
372   centre = 0.;
373   Int_t nn=0;
374   for(Int_t ii=ii1;ii<ii2;ii++){
375     centre+= fZCombf->GetBinCenter(ii)*fZCombf->GetBinContent(ii);
376     nn+=static_cast<Int_t>(fZCombf->GetBinContent(ii));
377   }
378   if (nn) centre/=nn;
379   /*
380   if(fDebug>0){
381     cout<<"Value of center "<<centre<<endl;
382     cout<<"Population of 3 central bins: "<<hc->GetBinContent(bi-1)<<", ";
383     cout<<hc->GetBinContent(bi)<<", ";
384     cout<<hc->GetBinContent(bi+1)<<endl;
385   }
386   */
387   // n1 is the bin number of fine histogram containing the point located 1 coarse bin less than "centre"
388   Int_t n1 = static_cast<Int_t>((centre-hc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fStepFine);
389   // n2 is the bin number of fine histogram containing the point located 1 coarse bin more than "centre"
390   Int_t n2 = static_cast<Int_t>((centre+hc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fStepFine);
391   if(n1<1)n1=1;
392   if(n2>nbinfine)n2=nbinfine;
393   Int_t niter = 0;
394   goon = kTRUE;
395   Int_t num=0;
396   Bool_t last = kFALSE;
397
398   while(goon || last){
399     fZFound = 0.;
400     fZsig = 0.;
401     num=0;
402     // at the end of the loop:
403     // fZFound = N*(Average Z) - where N is the number of tracklets
404     // num=N
405     // fZsig = N*Q - where Q is the average of Z*Z
406     for(Int_t n=n1;n<=n2;n++){
407       fZFound+=fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
408       num+=static_cast<Int_t>(fZCombf->GetBinContent(n));
409       fZsig+=fZCombf->GetBinCenter(n)*fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
410     }
411     if(num<2){
412       fZsig = 5.3;  // Default error from the beam sigmoid
413     }
414     else{
415       Float_t radi =  fZsig/(num-1)-fZFound*fZFound/num/(num-1);
416       // radi = square root of sample variance of Z
417       if(radi>0.)fZsig=TMath::Sqrt(radi);
418       else fZsig=5.3;  // Default error from the beam sigmoid
419       // fZfound - Average Z
420       fZFound/=num;
421     }
422     if(!last){
423       // goon is true if the distance between the found Z and the lower bin differs from the distance between the found Z and
424       // the upper bin by more than tolerance (0.002)
425       goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance;
426       // a window in the fine grained histogram is centered aroung the found Z. The width is 2 bins of
427       // the coarse grained histogram
428       if(num>0){
429         n1 = static_cast<Int_t>((fZFound-hc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fStepFine);
430         if(n1<1)n1=1;
431         n2 = static_cast<Int_t>((fZFound+hc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fStepFine);
432         if(n2>nbinfine)n2=nbinfine;
433         /*
434           if(fDebug>0){
435           cout<<"Search restricted to n1 = "<<n1<<", n2= "<<n2<<endl;
436           cout<<"z1= "<<fZCombf->GetBinCenter(n1)<<", n2 = "<<fZCombf->GetBinCenter(n2)<<endl;
437           }
438         */
439       }
440       else {
441         n1 = static_cast<Int_t>((centre-(niter+2)*hc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fStepFine);
442         n2 = static_cast<Int_t>((centre+(niter+2)*hc->GetBinWidth(bi)-fZCombf->GetBinLowEdge(0))/fStepFine);
443         if(n1<1)n1=1;
444         if(n2>nbinfine)n2=nbinfine;
445       }
446       niter++;
447       // no more than 10 adjusting iterations
448       if(niter>=10){
449         goon = kFALSE;
450         //      Warning("FindVertexForCurrentEvent","The procedure does not converge\n");
451       }
452
453       if((fZsig> 0.0001) && !goon && num>5){
454         last = kTRUE;
455         n1 = static_cast<Int_t>((fZFound-fZsig-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
456         if(n1<1)n1=1;
457         n2 = static_cast<Int_t>((fZFound+fZsig-fZCombf->GetBinLowEdge(0))/fZCombf->GetBinWidth(0));
458         if(n2>nbinfine)n2=nbinfine;
459         /*
460         if(fDebug>0){
461           cout<<"FINAL Search restricted to n1 = "<<n1<<", n2= "<<n2<<endl;
462           cout<<"z1= "<<fZCombf->GetBinCenter(n1)<<", n2 = "<<fZCombf->GetBinCenter(n2)<<endl;
463         }
464         */
465       }
466     }
467     else {
468       last = kFALSE;
469     }
470
471   }
472   //  if(fDebug>0)cout<<"Numer of Iterations "<<niter<<endl<<endl;
473   //  if(num!=0)fZsig/=TMath::Sqrt(num);
474   if (fZsig<=0) fZsig=5.3; // Default error from the beam sigmoid
475   fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
476   fCurrentVertex->SetTitle("vertexer: B");
477   ResetHistograms();
478   itsLoader->UnloadRecPoints();
479   return;
480 }
481
482 //_____________________________________________________________________
483 void AliITSVertexerZ::ResetHistograms(){
484   // delete TH1 data members
485   if(fZCombc)delete fZCombc;
486   if(fZCombf)delete fZCombf;
487   if(fZCombv)delete fZCombv;
488   fZCombc = 0;
489   fZCombf = 0;
490   fZCombv = 0;
491 }
492
493 //______________________________________________________________________
494 void AliITSVertexerZ::FindVertices(){
495   // computes the vertices of the events in the range FirstEvent - LastEvent
496   AliRunLoader *rl = AliRunLoader::GetRunLoader();
497   AliITSLoader* itsLoader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
498   itsLoader->ReloadRecPoints();
499   for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
500     //  cout<<"Processing event "<<i<<endl;
501     rl->GetEvent(i);
502     FindVertexForCurrentEvent(i);
503     if(fCurrentVertex){
504       WriteCurrentVertex();
505     }
506     /*
507     else {
508       if(fDebug>0){
509         cout<<"Vertex not found for event "<<i<<endl;
510         cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
511       }
512     }
513     */
514   }
515 }
516
517 //________________________________________________________
518 void AliITSVertexerZ::PrintStatus() const {
519   // Print current status
520   cout <<"=======================================================\n";
521   cout <<" First layer first and last modules: "<<fFirstL1<<", ";
522   cout <<fLastL1<<endl;
523   cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
524   cout <<fLastL2<<endl;
525   cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
526   cout <<"Limits for Z histograms: "<<fLowLim<<"; "<<fHighLim<<endl;
527   cout <<"Bin sizes for coarse and fine z histos "<<fStepCoarse<<"; "<<fStepFine<<endl;
528   cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
529   cout <<" Debug flag: "<<fDebug<<endl;
530   cout <<"First event to be processed "<<fFirstEvent;
531   cout <<"\n Last event to be processed "<<fLastEvent<<endl;
532   if(fZCombc){
533     cout<<"fZCombc exists - entries="<<fZCombc->GetEntries()<<endl;
534   }
535   else{
536     cout<<"fZCombc does not exist\n";
537   }
538   if(fZCombv){
539     cout<<"fZCombv exists - entries="<<fZCombv->GetEntries()<<endl;
540   }
541   else{
542     cout<<"fZCombv does not exist\n";
543   }
544    if(fZCombf){
545     cout<<"fZCombf exists - entries="<<fZCombv->GetEntries()<<endl;
546   }
547   else{
548     cout<<"fZCombf does not exist\n";
549   }
550  
551   cout <<"=======================================================\n";
552 }
553