4502f240c60f200b6d518791a53a57125f1e7ac5
[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 Int_t AliITSVertexerZ::GetPeakRegion(TH1F*h, Int_t &binmin, Int_t &binmax) const {
161   // Finds a region around a peak in the Z histogram
162   // Case of 2 peaks is treated 
163   Int_t imax=h->GetNbinsX();
164   Float_t maxval=0;
165   Int_t bi1=h->GetMaximumBin();
166   Int_t bi2=0;
167   for(Int_t i=imax;i>=1;i--){
168     if(h->GetBinContent(i)>maxval){
169       maxval=h->GetBinContent(i);
170       bi2=i;
171     }
172   }
173   Int_t npeaks=0;
174
175   if(bi1==bi2){
176     binmin=bi1-3;
177     binmax=bi1+3;
178     npeaks=1;
179   }else{
180     TH1F *copy = new TH1F(*h);
181     copy->SetBinContent(bi1,0.);
182     copy->SetBinContent(bi2,0.);
183     Int_t l1=TMath::Max(bi1-3,1);
184     Int_t l2=TMath::Min(bi1+3,h->GetNbinsX());
185     Float_t cont1=copy->Integral(l1,l2);
186     Int_t ll1=TMath::Max(bi2-3,1);
187     Int_t ll2=TMath::Min(bi2+3,h->GetNbinsX());
188     Float_t cont2=copy->Integral(ll1,ll2);
189     if(cont1>cont2){
190       binmin=l1;
191       binmax=l2;
192       npeaks=1;
193     }
194     if(cont2>cont1){
195       binmin=ll1;
196       binmax=ll2;
197       npeaks=1;
198     }
199     if(cont1==cont2){
200       binmin=l1;
201       binmax=ll2;
202       if(bi2-bi1==1) npeaks=1;
203       else npeaks=2;
204     }  
205     delete copy;
206   }    
207   return npeaks;
208 }
209 //______________________________________________________________________
210 AliESDVertex* AliITSVertexerZ::FindVertexForCurrentEvent(Int_t evnumber){
211   // Defines the AliESDVertex for the current event
212   VertexZFinder(evnumber);
213   Int_t ntrackl=0;
214   for(Int_t iteraz=0;iteraz<fMaxIter;iteraz++){
215     if(fCurrentVertex) ntrackl=fCurrentVertex->GetNContributors();
216     if(!fCurrentVertex || ntrackl==0 || ntrackl==-1){
217       Float_t diffPhiMaxOrig=fDiffPhiMax;
218       fDiffPhiMax=GetPhiMaxIter(iteraz);
219       VertexZFinder(evnumber);
220       fDiffPhiMax=diffPhiMaxOrig;
221     }
222   }
223   FindMultiplicity(evnumber);
224   return fCurrentVertex;
225 }  
226
227
228
229
230 //______________________________________________________________________
231 void AliITSVertexerZ::VertexZFinder(Int_t evnumber){
232   // Defines the AliESDVertex for the current event
233   fCurrentVertex = 0;
234   AliRunLoader *rl =AliRunLoader::GetRunLoader();
235   AliITSLoader* itsLoader = (AliITSLoader*)rl->GetLoader("ITSLoader");
236   AliITSgeom* geom = itsLoader->GetITSgeom();
237   itsLoader->LoadRecPoints();
238   rl->GetEvent(evnumber);
239
240   AliITSDetTypeRec detTypeRec;
241
242   TTree *tR = itsLoader->TreeR();
243   detTypeRec.SetTreeAddressR(tR);
244   TClonesArray *itsRec  = 0;
245   // lc and gc are local and global coordinates for layer 1
246   Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
247   Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
248   // lc2 and gc2 are local and global coordinates for layer 2
249   Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
250   Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
251
252   itsRec = detTypeRec.RecPoints();
253   TBranch *branch;
254   branch = tR->GetBranch("ITSRecPoints");
255
256   Int_t nrpL1 = 0;
257   Int_t nrpL2 = 0;
258   // By default fFirstL1=0 and fLastL1=79
259   // This loop counts the number of recpoints on layer1 (central modules)
260   for(Int_t module= fFirstL1; module<=fLastL1;module++){
261     // Keep only central modules
262     if(module%4==0 || module%4==3)continue;
263     //   cout<<"Procesing module "<<module<<" ";
264     branch->GetEvent(module);
265     //    cout<<"Number of clusters "<<clusters->GetEntries()<<endl;
266     nrpL1+= itsRec->GetEntries();
267     detTypeRec.ResetRecPoints();
268   }
269   //By default fFirstL2=80 and fLastL2=239
270   //This loop counts the number of RP on layer 2
271   for(Int_t module= fFirstL2; module<=fLastL2;module++){
272     branch->GetEvent(module);
273     nrpL2+= itsRec->GetEntries();
274     detTypeRec.ResetRecPoints();
275   }
276   if(nrpL1 == 0 || nrpL2 == 0){
277     ResetHistograms();
278     itsLoader->UnloadRecPoints();
279     fCurrentVertex = new AliESDVertex(0.,5.3,-2);
280     return;
281   }
282   // The vertex finding is attempted only if the number of RP is !=0 on
283   // both layers
284   Float_t *xc1 = new Float_t [nrpL1]; // coordinates of the L1 Recpoints
285   Float_t *yc1 = new Float_t [nrpL1];
286   Float_t *zc1 = new Float_t [nrpL1];
287   Float_t *phi1 = new Float_t [nrpL1];
288   Float_t *xc2 = new Float_t [nrpL2]; // coordinates of the L1 Recpoints
289   Float_t *yc2 = new Float_t [nrpL2];
290   Float_t *zc2 = new Float_t [nrpL2];
291   Float_t *phi2 = new Float_t [nrpL2];
292   Int_t ind = 0;// running index for RP
293   // Force a coarse bin size of 200 microns if the number of clusters on layer 2
294   // is low
295   if(nrpL2<fPPsetting[0])SetBinWidthCoarse(fPPsetting[1]);
296   // By default nbinfine = (10+10)/0.0005=40000
297   Int_t nbinfine = static_cast<Int_t>((fHighLim-fLowLim)/fStepFine);
298   // By default nbincoarse=(10+10)/0.01=2000
299   Int_t nbincoarse = static_cast<Int_t>((fHighLim-fLowLim)/fStepCoarse);
300   // Set stepverycoarse = 3*fStepCoarse
301   Int_t nbinvcoarse = static_cast<Int_t>((fHighLim-fLowLim)/(3.*fStepCoarse));
302   if(fZCombc)delete fZCombc;
303   fZCombc = new TH1F("fZCombc","Z",nbincoarse,fLowLim,fLowLim+nbincoarse*fStepCoarse);
304   if(fZCombv)delete fZCombv;
305   fZCombv = new TH1F("fZCombv","Z",nbinvcoarse,fLowLim,fLowLim+nbinvcoarse*3.*fStepCoarse);
306   if(fZCombf)delete fZCombf;
307   fZCombf = new TH1F("fZCombf","Z",nbinfine,fLowLim,fLowLim+nbinfine*fStepFine);
308   // Loop on modules of layer 1 (restricted to central modules)
309   for(Int_t module= fFirstL1; module<=fLastL1;module++){
310     if(module%4==0 || module%4==3)continue;
311     branch->GetEvent(module);
312     Int_t nrecp1 = itsRec->GetEntries();
313     for(Int_t j=0;j<nrecp1;j++){
314       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
315       // Local coordinates of this recpoint
316       lc[0]=recp->GetDetLocalX();
317       lc[2]=recp->GetDetLocalZ();
318       geom->LtoG(module,lc,gc);
319       // Global coordinates of this recpoints
320       gc[0]-=fX0; // Possible beam offset in the bending plane
321       gc[1]-=fY0; //   "               "
322       xc1[ind]=gc[0];
323       yc1[ind]=gc[1];
324       zc1[ind]=gc[2];
325       // azimuthal angle is computed in the interval 0 --> 2*pi
326       phi1[ind] = TMath::ATan2(gc[1],gc[0]);
327       if(phi1[ind]<0)phi1[ind]=2*TMath::Pi()+phi1[ind];
328       ind++;
329     }
330     detTypeRec.ResetRecPoints();
331   }
332   ind = 0; // the running index is reset for Layer 2
333   for(Int_t module= fFirstL2; module<=fLastL2;module++){
334     branch->GetEvent(module);
335     Int_t nrecp2 = itsRec->GetEntries();
336     for(Int_t j=0;j<nrecp2;j++){
337       AliITSRecPoint *recp = (AliITSRecPoint*)itsRec->At(j);
338       lc[0]=recp->GetDetLocalX();
339       lc[2]=recp->GetDetLocalZ();
340       geom->LtoG(module,lc,gc);
341       gc[0]-=fX0;
342       gc[1]-=fY0;
343       xc2[ind]=gc[0];
344       yc2[ind]=gc[1];
345       zc2[ind]=gc[2];
346       phi2[ind] = TMath::ATan2(gc[1],gc[0]);
347       if(phi2[ind]<0)phi2[ind]=2*TMath::Pi()+phi2[ind];
348       ind++;
349     }
350     detTypeRec.ResetRecPoints();
351   }
352  
353   // Int_t nolines=0;
354   for(Int_t i=0;i<nrpL1;i++){ // loop on L1 RP
355     Float_t r1=TMath::Sqrt(xc1[i]*xc1[i]+yc1[i]*yc1[i]); // radius L1 RP
356     for(Int_t j=0;j<nrpL2;j++){ // loop on L2 RP
357       Float_t diff = TMath::Abs(phi2[j]-phi1[i]); // diff in azimuth
358       if(diff>TMath::Pi())diff=2.*TMath::Pi()-diff; //diff<pi
359       if(diff<fDiffPhiMax){ // cut on 10 milliradians by def.
360         Float_t r2=TMath::Sqrt(xc2[j]*xc2[j]+yc2[j]*yc2[j]); // radius L2 RP
361         Float_t zr0=(r2*zc1[i]-r1*zc2[j])/(r2-r1); //Z @ null radius
362         fZCombf->Fill(zr0);
363         fZCombc->Fill(zr0);
364         fZCombv->Fill(zr0);
365         Double_t pA[3];
366         Double_t pB[3];
367         pA[0]=xc1[i];
368         pA[1]=yc1[i];
369         pA[2]=zc1[i];
370         pB[0]=xc2[j];
371         pB[1]=yc2[j];
372         pB[2]=zc2[j];
373         //      MakeTracklet(pA,pB,nolines);
374       }
375     }
376   }
377   delete [] xc1;
378   delete [] yc1;
379   delete [] zc1;
380   delete [] phi1;
381   delete [] xc2;
382   delete [] yc2;
383   delete [] zc2;
384   delete [] phi2;
385   Double_t contents = fZCombc->GetEntries()- fZCombc->GetBinContent(0)-fZCombc->GetBinContent(nbincoarse+1);
386   if(contents<1.){
387     //    Warning("FindVertexForCurrentEvent","Insufficient number of rec. points\n");
388     ResetHistograms();
389     itsLoader->UnloadRecPoints();
390     fCurrentVertex = new AliESDVertex(0.,5.3,-1);
391     return;
392   }
393
394   TH1F *hc = fZCombc;
395
396   
397   if(hc->GetBinContent(hc->GetMaximumBin())<3)hc = fZCombv;
398   Int_t binmin,binmax;
399   Int_t nPeaks=GetPeakRegion(hc,binmin,binmax);   
400   if(nPeaks==2)AliWarning("2 peaks found");
401   Int_t ii1=1+static_cast<Int_t>((hc->GetBinLowEdge(binmin)-fLowLim)/fStepFine);
402   Int_t ii2=1+static_cast<Int_t>((hc->GetBinLowEdge(binmax)+hc->GetBinWidth(binmax)-fLowLim)/fStepFine);
403   Float_t centre = 0.;
404   Int_t nn=0;
405   for(Int_t ii=ii1;ii<=ii2;ii++){
406     centre+= fZCombf->GetBinCenter(ii)*fZCombf->GetBinContent(ii);
407     nn+=static_cast<Int_t>(fZCombf->GetBinContent(ii));
408   }
409   if (nn) centre/=nn;
410
411   // n1 is the bin number of fine histogram containing the point located 1 coarse bin less than "centre"
412   Int_t n1 = 1+static_cast<Int_t>((centre-hc->GetBinWidth(1)-fLowLim)/fStepFine);
413   // n2 is the bin number of fine histogram containing the point located 1 coarse bin more than "centre"
414   Int_t n2 = 1+static_cast<Int_t>((centre+hc->GetBinWidth(1)-fLowLim)/fStepFine);
415   if(n1<1)n1=1;
416   if(n2>nbinfine)n2=nbinfine;
417   Int_t niter = 0;
418   Bool_t goon = kTRUE;
419   Int_t num=0;
420   Bool_t last = kFALSE;
421
422   while(goon || last){
423     fZFound = 0.;
424     fZsig = 0.;
425     num=0;
426     // at the end of the loop:
427     // fZFound = N*(Average Z) - where N is the number of tracklets
428     // num=N
429     // fZsig = N*Q - where Q is the average of Z*Z
430     for(Int_t n=n1;n<=n2;n++){
431       fZFound+=fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
432       num+=static_cast<Int_t>(fZCombf->GetBinContent(n));
433       fZsig+=fZCombf->GetBinCenter(n)*fZCombf->GetBinCenter(n)*fZCombf->GetBinContent(n);
434     }
435     if(num<2){
436       fZsig = 5.3;  // Default error from the beam sigmoid
437     }
438     else{
439       Float_t radi =  fZsig/(num-1)-fZFound*fZFound/num/(num-1);
440       // radi = square root of sample variance of Z
441       if(radi>0.)fZsig=TMath::Sqrt(radi);
442       else fZsig=5.3;  // Default error from the beam sigmoid
443       // fZfound - Average Z
444       fZFound/=num;
445     }
446     if(!last){
447       // goon is true if the distance between the found Z and the lower bin differs from the distance between the found Z and
448       // the upper bin by more than tolerance (0.002)
449       goon = TMath::Abs(TMath::Abs(fZFound-fZCombf->GetBinCenter(n1))-TMath::Abs(fZFound-fZCombf->GetBinCenter(n2)))>fTolerance;
450       // a window in the fine grained histogram is centered aroung the found Z. The width is 2 bins of
451       // the coarse grained histogram
452       if(num>0){
453         n1 = 1+static_cast<Int_t>((fZFound-hc->GetBinWidth(1)-fLowLim)/fStepFine);
454         if(n1<1)n1=1;
455         n2 = 1+static_cast<Int_t>((fZFound+hc->GetBinWidth(1)-fLowLim)/fStepFine);
456         if(n2>nbinfine)n2=nbinfine;
457       }
458       else {
459         n1 = 1+static_cast<Int_t>((centre-(niter+2)*hc->GetBinWidth(1)-fLowLim)/fStepFine);
460         n2 = 1+static_cast<Int_t>((centre+(niter+2)*hc->GetBinWidth(1)-fLowLim)/fStepFine);
461         if(n1<1)n1=1;
462         if(n2>nbinfine)n2=nbinfine;
463       }
464       niter++;
465       // no more than 10 adjusting iterations
466       if(niter>=10)goon = kFALSE;
467
468       if((fZsig> 0.0001) && !goon && num>5){
469         last = kTRUE;
470         n1 = 1+static_cast<Int_t>((fZFound-fZsig-fLowLim)/fStepFine);
471         if(n1<1)n1=1;
472         n2 = 1+static_cast<Int_t>((fZFound+fZsig-fLowLim)/fStepFine);  
473         if(n2>nbinfine)n2=nbinfine;
474       }
475     }
476     else {
477       last = kFALSE;
478     }
479
480   }
481   //  if(fDebug>0)cout<<"Numer of Iterations "<<niter<<endl<<endl;
482   //  if(num!=0)fZsig/=TMath::Sqrt(num);
483   if (fZsig<=0) fZsig=5.3; // Default error from the beam sigmoid
484   fCurrentVertex = new AliESDVertex(fZFound,fZsig,num);
485   fCurrentVertex->SetTitle("vertexer: B");
486   ResetHistograms();
487   itsLoader->UnloadRecPoints();
488   return;
489 }
490
491 //_____________________________________________________________________
492 void AliITSVertexerZ::ResetHistograms(){
493   // delete TH1 data members
494   if(fZCombc)delete fZCombc;
495   if(fZCombf)delete fZCombf;
496   if(fZCombv)delete fZCombv;
497   fZCombc = 0;
498   fZCombf = 0;
499   fZCombv = 0;
500 }
501
502 //______________________________________________________________________
503 void AliITSVertexerZ::FindVertices(){
504   // computes the vertices of the events in the range FirstEvent - LastEvent
505   AliRunLoader *rl = AliRunLoader::GetRunLoader();
506   AliITSLoader* itsLoader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
507   itsLoader->ReloadRecPoints();
508   for(Int_t i=fFirstEvent;i<=fLastEvent;i++){
509     //  cout<<"Processing event "<<i<<endl;
510     rl->GetEvent(i);
511     FindVertexForCurrentEvent(i);
512     if(fCurrentVertex){
513       WriteCurrentVertex();
514     }
515     /*
516     else {
517       if(fDebug>0){
518         cout<<"Vertex not found for event "<<i<<endl;
519         cout<<"fZFound = "<<fZFound<<", fZsig= "<<fZsig<<endl;
520       }
521     }
522     */
523   }
524 }
525
526 //________________________________________________________
527 void AliITSVertexerZ::PrintStatus() const {
528   // Print current status
529   cout <<"=======================================================\n";
530   cout <<" First layer first and last modules: "<<fFirstL1<<", ";
531   cout <<fLastL1<<endl;
532   cout <<" Second layer first and last modules: "<<fFirstL2<<", ";
533   cout <<fLastL2<<endl;
534   cout <<" Max Phi difference: "<<fDiffPhiMax<<endl;
535   cout <<"Limits for Z histograms: "<<fLowLim<<"; "<<fHighLim<<endl;
536   cout <<"Bin sizes for coarse and fine z histos "<<fStepCoarse<<"; "<<fStepFine<<endl;
537   cout <<" Current Z "<<fZFound<<"; Z sig "<<fZsig<<endl;
538   cout <<" Debug flag: "<<fDebug<<endl;
539   cout <<"First event to be processed "<<fFirstEvent;
540   cout <<"\n Last event to be processed "<<fLastEvent<<endl;
541   if(fZCombc){
542     cout<<"fZCombc exists - entries="<<fZCombc->GetEntries()<<endl;
543   }
544   else{
545     cout<<"fZCombc does not exist\n";
546   }
547   if(fZCombv){
548     cout<<"fZCombv exists - entries="<<fZCombv->GetEntries()<<endl;
549   }
550   else{
551     cout<<"fZCombv does not exist\n";
552   }
553    if(fZCombf){
554     cout<<"fZCombf exists - entries="<<fZCombv->GetEntries()<<endl;
555   }
556   else{
557     cout<<"fZCombf does not exist\n";
558   }
559  
560   cout <<"=======================================================\n";
561 }
562