Fixes from Marek to bring the secondary correction object
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliCentralMCCorrectionsTask.cxx
1 // 
2 // Calculate the corrections in the central regions
3 // 
4 // Inputs: 
5 //   - AliESDEvent 
6 //
7 // Outputs: 
8 //   - AliAODCentralMult 
9 // 
10 // Histograms 
11 //   
12 // Corrections used 
13 // 
14 #include "AliCentralMCCorrectionsTask.h"
15 #include "AliTriggerAnalysis.h"
16 #include "AliPhysicsSelection.h"
17 #include "AliLog.h"
18 #include "AliHeader.h"
19 #include "AliGenEventHeader.h"
20 #include "AliESDEvent.h"
21 #include "AliAODHandler.h"
22 #include "AliMultiplicity.h"
23 #include "AliInputEventHandler.h"
24 #include "AliStack.h"
25 #include "AliMCEvent.h"
26 #include "AliAODForwardMult.h"
27 #include "AliCentralCorrSecondaryMap.h"
28 #include "AliCentralCorrAcceptance.h"
29 #include "AliForwardUtil.h"
30 #include "AliMultiplicity.h"
31 #include <TH1.h>
32 #include <TH2D.h>
33 #include <TDirectory.h>
34 #include <TList.h>
35 #include <TROOT.h>
36 #include <iostream>
37
38 //====================================================================
39 namespace {
40   const char* GetEventName(Bool_t tr, Bool_t vtx) 
41   {
42     return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx" : ""));
43   }
44 }
45
46 //====================================================================
47 AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask()
48   : AliAnalysisTaskSE(),
49     fInspector(),
50     fTrackDensity(),
51     fVtxBins(0),
52     fFirstEvent(true),
53     fHEvents(0), 
54     fHEventsTr(0), 
55     fHEventsTrVtx(0),
56     fVtxAxis(),
57     fEtaAxis(),
58     fList(),
59     fNPhiBins(20),
60     fEffectiveCorr(true)
61 {
62   // 
63   // Constructor 
64   // 
65   // Parameters:
66   //    name Name of task 
67   //
68   DGUARD(fDebug,0,"Default construction of AliCentralMCCorrectionsTask");
69 }
70
71 //____________________________________________________________________
72 AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const char* name)
73   : AliAnalysisTaskSE(name),
74     fInspector("eventInspector"), 
75     fTrackDensity("trackDensity"),
76     fVtxBins(0),
77     fFirstEvent(true),
78     fHEvents(0), 
79     fHEventsTr(0), 
80     fHEventsTrVtx(0),
81     fVtxAxis(10,-10,10), 
82     fEtaAxis(200,-4,6),
83     fList(),
84     fNPhiBins(20),
85     fEffectiveCorr(true)
86 {
87   // 
88   // Constructor 
89   // 
90   // Parameters:
91   //    name Name of task 
92   //
93   DGUARD(fDebug,0,"Named construction of AliCentralMCCorrectionsTask: %s",name);
94   DefineOutput(1, TList::Class());
95   DefineOutput(2, TList::Class());
96 }
97
98 //____________________________________________________________________
99 AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const AliCentralMCCorrectionsTask& o)
100   : AliAnalysisTaskSE(o),
101     fInspector(o.fInspector),
102     fTrackDensity(),
103     fVtxBins(0),
104     fFirstEvent(o.fFirstEvent),
105     fHEvents(o.fHEvents), 
106     fHEventsTr(o.fHEventsTr), 
107     fHEventsTrVtx(o.fHEventsTrVtx),
108     fVtxAxis(10,-10,10), 
109     fEtaAxis(200,-4,6),
110     fList(o.fList),
111     fNPhiBins(o.fNPhiBins),
112     fEffectiveCorr(o.fEffectiveCorr)
113 {
114   // 
115   // Copy constructor 
116   // 
117   // Parameters:
118   //    o Object to copy from 
119   //
120   DGUARD(fDebug,0,"Copy construction of AliCentralMCCorrectionsTask");
121 }
122
123 //____________________________________________________________________
124 AliCentralMCCorrectionsTask&
125 AliCentralMCCorrectionsTask::operator=(const AliCentralMCCorrectionsTask& o)
126 {
127   // 
128   // Assignment operator 
129   // 
130   // Parameters:
131   //    o Object to assign from 
132   // 
133   // Return:
134   //    Reference to this object 
135   //
136   DGUARD(fDebug,3,"Assignment of AliCentralMCCorrectionsTask");
137   if (&o == this) return *this; 
138   fInspector         = o.fInspector;
139   fTrackDensity      = o.fTrackDensity;
140   fVtxBins           = o.fVtxBins;
141   fFirstEvent        = o.fFirstEvent;
142   fHEvents           = o.fHEvents;
143   fHEventsTr         = o.fHEventsTr;
144   fHEventsTrVtx      = o.fHEventsTrVtx;
145   SetVertexAxis(o.fVtxAxis);
146   SetEtaAxis(o.fEtaAxis);
147   fNPhiBins          = o.fNPhiBins;
148   fEffectiveCorr     = o.fEffectiveCorr;
149
150   return *this;
151 }
152
153 //____________________________________________________________________
154 void
155 AliCentralMCCorrectionsTask::Init()
156 {
157   // 
158   // Initialize the task 
159   // 
160   //
161   DGUARD(fDebug,1,"Initialize AliCentralMCCorrectionsTask");
162 }
163
164 //____________________________________________________________________
165 void
166 AliCentralMCCorrectionsTask::SetVertexAxis(Int_t nBin, Double_t min, 
167                                            Double_t max)
168 {
169   // 
170   // Set the vertex axis to use
171   // 
172   // Parameters:
173   //    nBins Number of bins
174   //    vzMin Least @f$z@f$ coordinate of interation point
175   //    vzMax Largest @f$z@f$ coordinate of interation point
176   //
177   DGUARD(fDebug,3,"Set vertex axis AliCentralMCCorrectionsTask [%d,%f,%f]",
178          nBin, min, max);
179   if (max < min) { 
180     Double_t tmp = min;
181     min          = max;
182     max          = tmp;
183   }
184   if (min < -10) 
185     AliWarning(Form("Minimum vertex %f < -10, make sure you want this",min));
186   if (max > +10) 
187     AliWarning(Form("Minimum vertex %f > +10, make sure you want this",max));
188   fVtxAxis.Set(nBin, min, max);
189 }
190 //____________________________________________________________________
191 void
192 AliCentralMCCorrectionsTask::SetVertexAxis(const TAxis& axis)
193 {
194   // 
195   // Set the vertex axis to use
196   // 
197   // Parameters:
198   //    axis Axis
199   //
200   SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
201 }
202
203 //____________________________________________________________________
204 void
205 AliCentralMCCorrectionsTask::SetEtaAxis(Int_t nBin, Double_t min, Double_t max)
206 {
207   // 
208   // Set the eta axis to use
209   // 
210   // Parameters:
211   //    nBins Number of bins
212   //    vzMin Least @f$\eta@f$ 
213   //    vzMax Largest @f$\eta@f$ 
214   //
215   DGUARD(fDebug,3,"Set eta axis AliCentralMCCorrectionsTask [%d,%f,%f]",
216          nBin, min, max);
217   if (max < min) { 
218     Double_t tmp = min;
219     min          = max;
220     max          = tmp;
221   }
222   if (min < -4) 
223     AliWarning(Form("Minimum eta %f < -4, make sure you want this",min));
224   if (max > +6) 
225     AliWarning(Form("Minimum eta %f > +6, make sure you want this",max));
226   fEtaAxis.Set(nBin, min, max);
227 }
228 //____________________________________________________________________
229 void
230 AliCentralMCCorrectionsTask::SetEtaAxis(const TAxis& axis)
231 {
232   // 
233   // Set the eta axis to use
234   // 
235   // Parameters:
236   //    axis Axis
237   //
238   SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
239 }
240
241 //____________________________________________________________________
242 void
243 AliCentralMCCorrectionsTask::DefineBins(TList* l)
244 {
245   DGUARD(fDebug,1,"Define bins in AliCentralMCCorrectionsTask");
246   if (!fVtxBins) fVtxBins = new TObjArray(fVtxAxis.GetNbins(), 1);
247   if (fVtxBins->GetEntries() > 0) return;
248
249   fVtxBins->SetOwner();
250   for (Int_t i = 1; i <= fVtxAxis.GetNbins(); i++) { 
251     Double_t low  = fVtxAxis.GetBinLowEdge(i);
252     Double_t high = fVtxAxis.GetBinUpEdge(i);
253     VtxBin*  bin  = new VtxBin(low, high, fEtaAxis, fNPhiBins);
254     fVtxBins->AddAt(bin, i);
255     bin->DefineOutput(l);
256   }
257 }
258
259 //____________________________________________________________________
260 void
261 AliCentralMCCorrectionsTask::UserCreateOutputObjects()
262 {
263   // 
264   // Create output objects 
265   // 
266   //
267   DGUARD(fDebug,1,"Create user output for AliCentralMCCorrectionsTask");
268   fList = new TList;
269   fList->SetOwner();
270   fList->SetName(Form("%sSums", GetName()));
271
272   DefineBins(fList);
273
274   fHEvents = new TH1I(GetEventName(false,false),
275                       "Number of all events", 
276                       fVtxAxis.GetNbins(), 
277                       fVtxAxis.GetXmin(), 
278                       fVtxAxis.GetXmax());
279   fHEvents->SetXTitle("v_{z} [cm]");
280   fHEvents->SetYTitle("# of events");
281   fHEvents->SetFillColor(kBlue+1);
282   fHEvents->SetFillStyle(3001);
283   fHEvents->SetDirectory(0);
284   fList->Add(fHEvents);
285
286   fHEventsTr = new TH1I(GetEventName(true, false), 
287                         "Number of triggered events",
288                         fVtxAxis.GetNbins(), 
289                         fVtxAxis.GetXmin(), 
290                         fVtxAxis.GetXmax());
291   fHEventsTr->SetXTitle("v_{z} [cm]");
292   fHEventsTr->SetYTitle("# of events");
293   fHEventsTr->SetFillColor(kRed+1);
294   fHEventsTr->SetFillStyle(3001);
295   fHEventsTr->SetDirectory(0);
296   fList->Add(fHEventsTr);
297
298   fHEventsTrVtx = new TH1I(GetEventName(true,true),
299                            "Number of events w/trigger and vertex", 
300                            fVtxAxis.GetNbins(), 
301                            fVtxAxis.GetXmin(), 
302                            fVtxAxis.GetXmax());
303   fHEventsTrVtx->SetXTitle("v_{z} [cm]");
304   fHEventsTrVtx->SetYTitle("# of events");
305   fHEventsTrVtx->SetFillColor(kBlue+1);
306   fHEventsTrVtx->SetFillStyle(3001);
307   fHEventsTrVtx->SetDirectory(0);
308   fList->Add(fHEventsTrVtx);
309
310   // Copy axis objects to output 
311   TH1* vtxAxis = new TH1D("vtxAxis", "Vertex axis", 
312                           fVtxAxis.GetNbins(), 
313                           fVtxAxis.GetXmin(), 
314                           fVtxAxis.GetXmax());
315   TH1* etaAxis = new TH1D("etaAxis", "Eta axis", 
316                           fEtaAxis.GetNbins(), 
317                           fEtaAxis.GetXmin(), 
318                           fEtaAxis.GetXmax());
319   fList->Add(vtxAxis);
320   fList->Add(etaAxis);
321
322   AliInfo(Form("Initialising sub-routines: %p, %p", 
323                &fInspector, &fTrackDensity));
324   fInspector.DefineOutput(fList);
325   fInspector.Init(fVtxAxis);
326   fTrackDensity.DefineOutput(fList);
327
328   PostData(1, fList);
329 }
330 //____________________________________________________________________
331 void
332 AliCentralMCCorrectionsTask::UserExec(Option_t*)
333 {
334   // 
335   // Process each event 
336   // 
337   // Parameters:
338   //    option Not used
339   //  
340
341   DGUARD(fDebug,1,"AliCentralMCCorrectionsTask process an event");
342   // Get the input data - MC event
343   AliMCEvent*  mcEvent = MCEvent();
344   if (!mcEvent) { 
345     AliWarning("No MC event found");
346     return;
347   }
348
349   // Get the input data - ESD event
350   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
351   if (!esd) { 
352     AliWarning("No ESD event found for input event");
353     return;
354   }
355
356   //--- Read run information -----------------------------------------
357   if (fFirstEvent && esd->GetESDRun()) {
358     fInspector.ReadRunDetails(esd);
359     
360     AliInfo(Form("Initializing with parameters from the ESD:\n"
361                  "         AliESDEvent::GetBeamEnergy()   ->%f\n"
362                  "         AliESDEvent::GetBeamType()     ->%s\n"
363                  "         AliESDEvent::GetCurrentL3()    ->%f\n"
364                  "         AliESDEvent::GetMagneticField()->%f\n"
365                  "         AliESDEvent::GetRunNumber()    ->%d\n",
366                  esd->GetBeamEnergy(),
367                  esd->GetBeamType(),
368                  esd->GetCurrentL3(),
369                  esd->GetMagneticField(),
370                  esd->GetRunNumber()));
371
372     Print();
373     fFirstEvent = false;
374   }
375
376   // Some variables 
377   UInt_t   triggers; // Trigger bits
378   Bool_t   lowFlux;  // Low flux flag
379   UShort_t iVz;      // Vertex bin from ESD
380   Double_t vZ;       // Z coordinate from ESD
381   Double_t cent;     // Centrality 
382   UShort_t iVzMc;    // Vertex bin from MC
383   Double_t vZMc;     // Z coordinate of IP vertex from MC
384   Double_t b;        // Impact parameter
385   Double_t cMC;      // Centrality estimate from b
386   Int_t    nPart;    // Number of participants 
387   Int_t    nBin;     // Number of binary collisions 
388   Double_t phiR;     // Reaction plane from MC
389   UShort_t nClusters;// Number of SPD clusters 
390   // Process the data 
391   UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ, 
392                                      cent, nClusters);
393   fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc, 
394                        b, cMC, nPart, nBin, phiR);
395
396   Bool_t isInel   = triggers & AliAODForwardMult::kInel;
397   Bool_t hasVtx   = retESD == AliFMDMCEventInspector::kOk;
398
399   // Fill the event count histograms 
400   if (isInel)           fHEventsTr->Fill(vZMc);
401   if (isInel && hasVtx) fHEventsTrVtx->Fill(vZMc);
402   fHEvents->Fill(vZMc);
403
404   // Now find our vertex bin object 
405   VtxBin* bin = 0;
406   if (iVzMc > 0 && iVzMc <= fVtxAxis.GetNbins()) 
407     bin = static_cast<VtxBin*>(fVtxBins->At(iVzMc));
408   if (!bin) { 
409     // AliError(Form("No vertex bin object @ %d (%f)", iVzMc, vZMc));
410     return;
411   }
412
413   // Now process our input data and store in own ESD object 
414   fTrackDensity.Calculate(*mcEvent, vZMc, *bin->fHits, bin->fPrimary);
415
416   // Get the ESD object
417   const AliMultiplicity* spdmult = esd->GetMultiplicity();
418
419   // Count number of tracklets per bin 
420   for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++) 
421     bin->fClusters->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
422   //...and then the unused clusters in layer 1 
423   for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
424      Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
425      bin->fClusters->Fill(eta, spdmult->GetPhiSingle(j));
426   }
427
428   // Count events  
429   bin->fCounts->Fill(0.5);
430   
431   PostData(1, fList);
432 }
433
434 //____________________________________________________________________
435 void
436 AliCentralMCCorrectionsTask::Terminate(Option_t*)
437 {
438   // 
439   // End of job
440   // 
441   // Parameters:
442   //    option Not used 
443   //
444   DGUARD(fDebug,1,"AliCentralMCCorrectionsTask analyse merged output");
445
446   fList = dynamic_cast<TList*>(GetOutputData(1));
447   if (!fList) {
448     AliError("No output list defined");
449     return;
450   }
451
452   DefineBins(fList);
453
454   // Output list 
455   TList* output = new TList;
456   output->SetOwner();
457   output->SetName(Form("%sResults", GetName()));
458
459   // --- Fill correction object --------------------------------------
460   AliCentralCorrSecondaryMap* corr = new AliCentralCorrSecondaryMap;
461   corr->SetVertexAxis(fVtxAxis);
462   // corr->SetEtaAxis(fEtaAxis);
463
464  AliCentralCorrAcceptance* acorr = new AliCentralCorrAcceptance;
465   acorr->SetVertexAxis(fVtxAxis);
466
467   TIter     next(fVtxBins);
468   VtxBin*   bin = 0;
469   UShort_t  iVz = 1;
470   while ((bin = static_cast<VtxBin*>(next()))) 
471     bin->Finish(fList, output, iVz++, fEffectiveCorr, corr,acorr);
472
473   output->Add(corr);
474  output->Add(acorr);
475
476   PostData(2, output);
477 }
478
479 //____________________________________________________________________
480 void
481 AliCentralMCCorrectionsTask::Print(Option_t* option) const
482 {
483   std::cout << ClassName() << "\n"
484             << "  Vertex bins:      " << fVtxAxis.GetNbins() << '\n'
485             << "  Vertex range:     [" << fVtxAxis.GetXmin() 
486             << "," << fVtxAxis.GetXmax() << "]\n"
487             << "  Eta bins:         " << fEtaAxis.GetNbins() << '\n'
488             << "  Eta range:        [" << fEtaAxis.GetXmin() 
489             << "," << fEtaAxis.GetXmax() << "]\n"
490             << "  # of phi bins:    " << fNPhiBins 
491             << std::endl;
492   gROOT->IncreaseDirLevel();
493   fInspector.Print(option);
494   fTrackDensity.Print(option);
495   gROOT->DecreaseDirLevel();
496 }
497
498 //====================================================================
499 const char*
500 AliCentralMCCorrectionsTask::VtxBin::BinName(Double_t low, 
501                                              Double_t high) 
502 {
503   static TString buf;
504   buf = Form("vtx%+05.1f_%+05.1f", low, high);
505   buf.ReplaceAll("+", "p");
506   buf.ReplaceAll("-", "m");
507   buf.ReplaceAll(".", "d");
508   return buf.Data();
509 }
510
511
512 //____________________________________________________________________
513 AliCentralMCCorrectionsTask::VtxBin::VtxBin()
514   : fHits(0), 
515     fClusters(0),
516     fPrimary(0),
517     fCounts(0)
518 {
519 }
520 //____________________________________________________________________
521 AliCentralMCCorrectionsTask::VtxBin::VtxBin(Double_t     low, 
522                                             Double_t     high, 
523                                             const TAxis& axis,
524                                             UShort_t     nPhi)
525   : TNamed(BinName(low, high), 
526            Form("%+5.1fcm<v_{z}<%+5.1fcm", low, high)),
527     fHits(0), 
528     fClusters(0),
529     fPrimary(0),
530     fCounts(0)
531 {
532   fPrimary = new TH2D("primary", "Primaries", 
533                       axis.GetNbins(), axis.GetXmin(), axis.GetXmax(), 
534                       nPhi, 0, 2*TMath::Pi());
535   fPrimary->SetXTitle("#eta");
536   fPrimary->SetYTitle("#varphi [radians]");
537   fPrimary->Sumw2();
538   fPrimary->SetDirectory(0);
539
540   fHits = static_cast<TH2D*>(fPrimary->Clone("hits"));
541   fHits->SetTitle("Hits");
542   fHits->SetDirectory(0);
543
544   fClusters = static_cast<TH2D*>(fPrimary->Clone("clusters"));
545   fClusters->SetTitle("Clusters");
546   fClusters->SetDirectory(0);
547
548   fCounts = new TH1D("counts", "Counts", 1, 0, 1);
549   fCounts->SetXTitle("Events");
550   fCounts->SetYTitle("# of Events");
551   fCounts->SetDirectory(0);
552 }
553
554 //____________________________________________________________________
555 AliCentralMCCorrectionsTask::VtxBin::VtxBin(const VtxBin& o)
556   : TNamed(o),
557     fHits(0),
558     fClusters(0),
559     fPrimary(0), 
560     fCounts(0)
561 {
562   if (o.fHits) {
563     fHits = static_cast<TH2D*>(o.fHits->Clone());
564     fHits->SetDirectory(0);
565   }
566   if (o.fClusters) {
567     fClusters = static_cast<TH2D*>(o.fClusters->Clone());
568     fClusters->SetDirectory(0);
569   }
570   if (o.fPrimary) {
571     fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
572     fPrimary->SetDirectory(0);
573   }
574   if (o.fCounts) {
575     fCounts = static_cast<TH1D*>(o.fCounts->Clone());
576     fCounts->SetDirectory(0);
577   }
578 }
579
580 //____________________________________________________________________
581 AliCentralMCCorrectionsTask::VtxBin&
582 AliCentralMCCorrectionsTask::VtxBin::operator=(const VtxBin& o)
583 {
584   if (&o == this) return *this; 
585   TNamed::operator=(o);
586   fHits     = 0;
587   fPrimary  = 0;
588   fClusters = 0;
589   fCounts   = 0;
590   if (o.fHits) {
591     fHits = static_cast<TH2D*>(o.fHits->Clone());
592     fHits->SetDirectory(0);
593   }
594   if (o.fClusters) {
595     fClusters = static_cast<TH2D*>(o.fClusters->Clone());
596     fClusters->SetDirectory(0);
597   }
598   if (o.fPrimary) {
599     fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
600     fPrimary->SetDirectory(0);
601   }
602   if (o.fCounts) {
603     fCounts = static_cast<TH1D*>(o.fCounts->Clone());
604     fCounts->SetDirectory(0);
605   }
606   return *this;
607 }
608
609 //____________________________________________________________________
610 void
611 AliCentralMCCorrectionsTask::VtxBin::DefineOutput(TList* l)
612 {
613   TList* d = new TList;
614   d->SetName(GetName());
615   d->SetOwner();
616   l->Add(d);
617
618   d->Add(fHits);
619   d->Add(fClusters);
620   d->Add(fPrimary);
621   d->Add(fCounts);
622 }
623 //____________________________________________________________________
624 void
625 AliCentralMCCorrectionsTask::VtxBin::Finish(const TList* input, 
626                                             TList* output, 
627                                             UShort_t iVz, 
628                                             Bool_t effectiveCorr,
629                                             AliCentralCorrSecondaryMap* map,
630                                             AliCentralCorrAcceptance* acorr)
631 {
632   TList* out = new TList;
633   out->SetName(GetName());
634   out->SetOwner();
635   output->Add(out);
636   
637   TList* l = static_cast<TList*>(input->FindObject(GetName()));
638   if (!l) { 
639     AliError(Form("List %s not found in %s", GetName(), input->GetName()));
640     return;
641   }
642
643
644   TH2D*   hits  = static_cast<TH2D*>(l->FindObject("hits"));
645   TH2D*   clus  = static_cast<TH2D*>(l->FindObject("clusters"));
646   TH2D*   prim  = static_cast<TH2D*>(l->FindObject("primary"));
647   if (!hits || !prim) {
648     AliError(Form("Missing histograms: %p, %p", hits, prim));
649     return;
650   }
651
652   TH2D* h = 0;
653   if (effectiveCorr) h = static_cast<TH2D*>(clus->Clone("bgCorr"));
654   else               h = static_cast<TH2D*>(hits->Clone("bgCorr"));
655   h->SetDirectory(0);
656   h->Divide(prim);
657
658   TH1D* acc = new TH1D(Form("SPDacc_vrtbin_%d",iVz),
659                           "Acceptance correction for SPD" ,
660                           fPrimary->GetXaxis()->GetNbins(), 
661                           fPrimary->GetXaxis()->GetXmin(), 
662                           fPrimary->GetXaxis()->GetXmax());
663   TH1F* accden = static_cast<TH1F*>(acc->Clone(Form("%s_den",
664                                                     acc->GetName())));
665
666   for(Int_t xx = 1; xx <=h->GetNbinsX(); xx++) {
667     for(Int_t yy = 1; yy <=h->GetNbinsY(); yy++) {
668       if(TMath::Abs(h->GetXaxis()->GetBinCenter(xx)) > 1.9) {
669         h->SetBinContent(xx,yy,0.); 
670         h->SetBinError(xx,yy,0.); 
671       }
672       if(h->GetBinContent(xx,yy) > 0.9) 
673         acc->Fill(h->GetXaxis()->GetBinCenter(xx));
674       else {
675         h->SetBinContent(xx,yy,0.); 
676         h->SetBinError(xx,yy,0.); 
677       }
678       accden->Fill(h->GetXaxis()->GetBinCenter(xx));
679         
680     }
681   }
682   acc->Divide(accden);
683
684   map->SetCorrection(iVz, h);
685   acorr->SetCorrection(iVz, acc);
686   out->Add(h);
687   out->Add(acc);
688 }
689
690 //
691 // EOF
692 //