3afd108a4aa0ec24905c92ade4e90536ee3c22ac
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliCentralMultiplicityTask.cxx
1 //====================================================================
2 // 
3 // Base class for classes that calculate the multiplicity in the
4 // central region event-by-event
5 // 
6 // Inputs: 
7 //   - AliESDEvent 
8 //
9 // Outputs: 
10 //   - AliAODCentralMult 
11 // 
12 // Histograms 
13 //   
14 // Corrections used 
15 #include "AliCentralMultiplicityTask.h"
16 #include "AliAODForwardMult.h"
17 #include "AliForwardUtil.h"
18 #include "AliLog.h"
19 #include "AliAODHandler.h"
20 #include "AliAnalysisManager.h"
21 #include "AliESDEvent.h"
22 #include "AliMultiplicity.h"
23 #include <TROOT.h>
24 #include <TSystem.h>
25 #include <TFile.h>
26 #include <TError.h>
27 #include <TSystem.h>
28 #include <iostream>
29 #include <iomanip>
30
31 //====================================================================
32 AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name) 
33   : AliAnalysisTaskSE(name),
34     fInspector("centralEventInspector"),
35     fData(0),
36     fList(0),
37     fHits(0),
38     fAODCentral(kFALSE),
39     fManager(),
40     fUseSecondary(true),
41     fUseAcceptance(true),
42     fFirstEventSeen(false), 
43     fIvz(0),
44     fNClusterTracklet(0),
45     fClusterPerTracklet(0),
46     fNCluster(0),
47     fNTracklet(0),
48     fEtaMin(0),
49     fEtaMax(0)
50 {
51   // 
52   // Constructor 
53   //   
54   DGUARD(fDebug, 3,"Named CTOR of AliCentralMultiplicityTask: %s", name);
55   DefineOutput(1, TList::Class());
56   fBranchNames = 
57     "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
58     "SPDVertex.,PrimaryVertex.";
59 }
60 //____________________________________________________________________
61 AliCentralMultiplicityTask::AliCentralMultiplicityTask() 
62   : AliAnalysisTaskSE(),
63     fInspector(),
64     fData(0),
65     fList(0),
66     fHits(0),
67     fAODCentral(),
68     fManager(),
69     fUseSecondary(true),
70     fUseAcceptance(true),
71     fFirstEventSeen(false), 
72     fIvz(0),
73     fNClusterTracklet(0),
74     fClusterPerTracklet(0),
75     fNCluster(0),
76     fNTracklet(0),
77     fEtaMin(0),
78     fEtaMax(0)
79 {
80   // 
81   // Constructor 
82   // 
83   DGUARD(fDebug, 3,"Default CTOR of AliCentralMultiplicityTask");
84 }
85 //____________________________________________________________________
86 AliCentralMultiplicityTask::AliCentralMultiplicityTask(const AliCentralMultiplicityTask& o)
87   : AliAnalysisTaskSE(o),
88     fInspector(o.fInspector),
89     fData(o.fData),
90     fList(o.fList),
91     fHits(o.fHits),
92     fAODCentral(o.fAODCentral),
93     fManager(o.fManager),
94     fUseSecondary(o.fUseSecondary),
95     fUseAcceptance(o.fUseAcceptance),
96     fFirstEventSeen(o.fFirstEventSeen), 
97     fIvz(0),
98     fNClusterTracklet(o.fNClusterTracklet),
99     fClusterPerTracklet(o.fClusterPerTracklet),
100     fNCluster(o.fNCluster),
101     fNTracklet(o.fNTracklet),
102     fEtaMin(o.fEtaMin),
103     fEtaMax(o.fEtaMax)      
104 {
105   //
106   // Copy constructor 
107   // 
108   DGUARD(fDebug, 3,"COPY CTOR of AliCentralMultiplicityTask");
109
110 }
111 //____________________________________________________________________
112 AliCentralMultiplicityTask&
113 AliCentralMultiplicityTask::operator=(const AliCentralMultiplicityTask& o)
114 {
115   // 
116   // Assignment operator 
117   //
118   DGUARD(fDebug,3,"Assignment of AliCentralMultiplicityTask");
119   if (&o == this) return *this; 
120   fInspector         = o.fInspector;
121   fData              = o.fData;
122   fList              = o.fList;
123   fHits              = o.fHits;
124   fAODCentral        = o.fAODCentral;
125   fManager           = o.fManager;
126   fUseSecondary      = o.fUseSecondary;
127   fUseAcceptance     = o.fUseAcceptance;
128   fFirstEventSeen    = o.fFirstEventSeen;
129   fIvz               = 0; 
130   fNClusterTracklet  = o.fNClusterTracklet;
131   fClusterPerTracklet= o.fClusterPerTracklet;
132   fNCluster          = o.fNCluster;
133   fNTracklet         = o.fNTracklet;
134   fEtaMin            = o.fEtaMin;
135   fEtaMax            = o.fEtaMax;
136   return *this;
137 }
138 //____________________________________________________________________
139 Bool_t 
140 AliCentralMultiplicityTask::Configure(const char* macro)
141 {
142   // --- Configure the task ------------------------------------------
143   TString macroPath(gROOT->GetMacroPath());
144   if (!macroPath.Contains("$(ALICE_ROOT)/PWGLF/FORWARD/analysis2")) { 
145     macroPath.Append(":$(ALICE_ROOT)/PWGLF/FORWARD/analysis2");
146     gROOT->SetMacroPath(macroPath);
147   }
148   const char* config = gSystem->Which(gROOT->GetMacroPath(),macro);
149   if (!config) {
150     AliWarningF("%s not found in %s", macro, gROOT->GetMacroPath());
151     return false;
152   }
153
154   AliInfoF("Loading configuration of '%s' from %s", ClassName(), config);
155   gROOT->Macro(Form("%s((AliCentralMultiplicityTask*)%p)", config, this));
156   delete config;
157   
158   return true;
159 }
160
161 //____________________________________________________________________
162 void AliCentralMultiplicityTask::UserCreateOutputObjects() 
163 {
164   // 
165   // Create output objects 
166   // 
167   //
168   DGUARD(fDebug,1,"Create user output in AliCentralMultiplicityTask");
169
170   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
171   AliAODHandler*      ah = 
172     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
173   if (ah) 
174 {
175    //AliFatal("No AOD output handler set in analysis manager");
176         TObject* obj = &fAODCentral;
177         ah->AddBranch("AliAODCentralMult", &obj);
178  } 
179   
180     
181   fList = new TList();
182   fList->SetOwner();
183
184   fInspector.CreateOutputObjects(fList);
185
186   PostData(1,fList);  
187 }
188
189 //____________________________________________________________________
190 AliESDEvent*
191 AliCentralMultiplicityTask::GetESDEvent()
192 {
193   //
194   // Get the ESD event. IF this is the first event, initialise
195   //
196   DGUARD(fDebug,1,"Get ESD event in AliCentralMultiplicityTask");
197   if (IsZombie()) return 0;
198   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
199   if (!esd) {
200     AliWarning("No ESD event found for input event");
201     return 0;
202   }
203
204   // Load in the data needed
205   LoadBranches();
206
207   // IF we've read the first event already, just return the event 
208   if (fFirstEventSeen) return esd;
209   
210   // Read the details of the rung 
211   fInspector.ReadRunDetails(esd);
212
213   // If we weren't initialised before (i.e., in the setup), do so now. 
214   if (!GetManager().IsInit()) {
215     GetManager().Init(fInspector.GetCollisionSystem(),
216                       fInspector.GetEnergy(),
217                       fInspector.GetField());
218     //AliInfo("Manager of corrections in AliCentralMultiplicityTask init");
219   }
220   Bool_t ok = true;
221   if (/*fUseSecondary &&*/ !GetManager().HasSecondaryCorrection()) {
222     ok = false;
223     AliError("No secondary correction defined!");
224   }
225   if (/*fUseAcceptance &&*/ !GetManager().HasAcceptanceCorrection()) {
226     ok = false;
227     AliError("No acceptance correction defined!");
228   }
229   // If the corrections are not seen, make this a zombie, and prevent
230   // further execution of this task.
231   if (!ok) { 
232     AliError("Missing corrections, make this a zombie");
233     SetZombie(true);
234     esd = 0;
235     fFirstEventSeen = true;
236     return esd;
237   }
238
239   // Check for existence and get secondary map 
240   AliCentralCorrSecondaryMap* secMap = GetManager().GetSecMap();
241   const TAxis& vaxis = secMap->GetVertexAxis();
242
243   FindEtaLimits();
244
245   fNClusterTracklet = new TH2D("nClusterVsnTracklet", 
246                                "Total number of cluster vs number of tracklets",
247                                100, 0, 10000, 100, 0, 10000);
248   fNClusterTracklet->SetDirectory(0);
249   fNClusterTracklet->SetXTitle("# of free clusters");
250   fNClusterTracklet->SetYTitle("# of tracklets");
251   fNClusterTracklet->SetStats(0);
252   fList->Add(fNClusterTracklet);
253
254   Int_t    nEta = 80;
255   Double_t lEta = 2;
256   fClusterPerTracklet = new TH2D("clusterPerTracklet", 
257                                  "N_{free cluster}/N_{tracklet} vs. #eta", 
258                                  nEta,-lEta,lEta, 101, -.05, 10.05);
259   fClusterPerTracklet->SetDirectory(0);
260   fClusterPerTracklet->SetXTitle("#eta");
261   fClusterPerTracklet->SetYTitle("N_{free cluster}/N_{tracklet}");
262   fClusterPerTracklet->SetStats(0);
263   fList->Add(fClusterPerTracklet);
264
265   // Cache histograms 
266   fNCluster = new TH1D("cacheCluster", "", nEta,-lEta,lEta);
267   fNCluster->SetDirectory(0);
268   fNCluster->Sumw2();
269                        
270   fNTracklet = new TH1D("cacheTracklet", "", nEta,-lEta,lEta);
271   fNTracklet->SetDirectory(0);
272   fNTracklet->Sumw2();
273
274   // Initialize the inspecto 
275   fInspector.SetupForData(vaxis);
276   fFirstEventSeen = kTRUE;
277
278   // Print some information 
279   Print();
280
281   return esd;
282 }
283 //____________________________________________________________________
284 void
285 AliCentralMultiplicityTask::MarkEventForStore() const
286 {
287   // Make sure the AOD tree is filled 
288   DGUARD(fDebug,1,"Mark AOD event for store in AliCentralMultiplicityTask");
289   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
290   AliAODHandler*      ah = 
291     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
292   if (ah)
293  {  
294     //AliFatal("No AOD output handler set in analysis manager");
295         ah->SetFillAOD(kTRUE);
296  }
297 }
298
299 //____________________________________________________________________
300 void AliCentralMultiplicityTask::FindEtaLimits()
301 {
302   // Find our pseudo-rapidity limits 
303   // 
304   // Uses the secondary map to do so.
305   DGUARD(fDebug,1,"Find eta limits in AliCentralMultiplicityTask");
306   AliCentralCorrSecondaryMap* secMap = GetManager().GetSecMap();
307
308   const TAxis& vaxis = secMap->GetVertexAxis();
309   
310   fEtaMin.Set(vaxis.GetNbins());
311   fEtaMax.Set(vaxis.GetNbins());
312   
313   fHits = new TList;
314   fHits->SetOwner();
315   fHits->SetName("hitMaps");
316   fList->Add(fHits);
317   
318   TList* secs = new TList;
319   secs->SetOwner();
320   secs->SetName("secondaryMaps");
321   fList->Add(secs);
322   unsigned short s = 1;
323   TH2D* hCoverage = new TH2D("coverage", "#eta coverage per v_{z}", 
324                              secMap->GetCorrection(s)->GetXaxis()->GetNbins(),
325                              secMap->GetCorrection(s)->GetXaxis()->GetXmin(),
326                              secMap->GetCorrection(s)->GetXaxis()->GetXmax(),
327                              vaxis.GetNbins(),vaxis.GetXmin(),vaxis.GetXmax());
328   hCoverage->SetDirectory(0);
329   hCoverage->SetXTitle("#eta");
330   hCoverage->SetYTitle("v_{z} [cm]");
331   hCoverage->SetZTitle("n_{bins}");
332   fList->Add(hCoverage);
333   
334   fAODCentral.Init(*(secMap->GetCorrection(s)->GetXaxis()));
335   
336   for (Int_t v = 1; v <= vaxis.GetNbins(); v++) { 
337     TH2D* corr = secMap->GetCorrection(UShort_t(v));
338     TH1D* proj = corr->ProjectionX(Form("secCor%02d", v));
339     proj->Scale(1. / corr->GetNbinsY());
340     proj->SetTitle(Form("Projection of secondary correction "
341                         "for %+5.1f<v_{z}<%+5.1f",
342                         vaxis.GetBinLowEdge(v), vaxis.GetBinUpEdge(v)));
343     proj->SetYTitle("#LT 2^{nd} correction#GT");
344     proj->SetDirectory(0);
345     proj->SetMarkerStyle(20);
346     proj->SetMarkerColor(kBlue+1);
347     secs->Add(proj);
348     
349     TH2D* obg = static_cast<TH2D*>(corr->Clone(Form("secCor2DFiducial%02d",v)));
350     obg->SetDirectory(0);
351     secs->Add(obg);
352     
353     TH1D* after = static_cast<TH1D*>(proj->Clone(Form("secCorFiducial%02d",v)));
354     after->SetDirectory(0);
355     after->SetMarkerColor(kRed+1);
356     secs->Add(after);
357     
358     TH2D* data = static_cast<TH2D*>(corr->Clone(Form("hitMap%02d",v)));
359     //d->SetTitle(Form("hitMap%02d",v));
360     data->SetTitle(Form("d^{2}N/d#eta d#phi "
361                         "for %+5.1f<v_{z}<%+5.1f",
362                         vaxis.GetBinLowEdge(v), vaxis.GetBinUpEdge(v)));
363     data->GetZaxis()->SetTitle("");
364     data->SetMarkerColor(kBlack);
365     data->SetMarkerStyle(1);
366     fHits->Add(data);
367     
368     TH1D* hAcceptance = fManager.GetAcceptanceCorrection(v);
369     TH1D* accClone   = static_cast<TH1D*>(hAcceptance->Clone(Form("acceptance%02d",v)));
370     secs->Add(accClone);
371     
372     // Double_t prev = 0;
373     for (Int_t e = 1; e <= proj->GetNbinsX(); e++) { 
374       Double_t c = proj->GetBinContent(e);
375       if (c > .5 /*&& TMath::Abs(c - prev) < .1*c*/) {
376         fEtaMin[v-1] = e;
377         break;
378       }
379       // prev = c;
380       after->SetBinContent(e, 0);
381       after->SetBinError(e, 0);
382       for(Int_t nn =1; nn <=obg->GetNbinsY();nn++)
383         obg->SetBinContent(e,nn,0);
384       
385     }
386     for (Int_t e = proj->GetNbinsX(); e >= 1; e--) { 
387       Double_t c = proj->GetBinContent(e);
388       if (c > .5 /*&& TMath::Abs(c - prev) < .1*c*/) {
389         fEtaMax[v-1] = e;
390         break;
391       }
392       // prev = c;
393       after->SetBinContent(e, 0);
394       after->SetBinError(e, 0);
395       for(Int_t nn =1; nn <=obg->GetNbinsY();nn++)
396         obg->SetBinContent(e,nn,0);
397       
398     }
399     
400     for (Int_t nn = fEtaMin[v-1]; nn<=fEtaMax[v-1]; nn++) { 
401       hCoverage->SetBinContent(nn,v,1);
402     }
403     
404   }
405 }
406
407 //____________________________________________________________________
408 void AliCentralMultiplicityTask::UserExec(Option_t* /*option*/) 
409 {
410   // 
411   // Process each event 
412   // 
413   // Parameters:
414   //    option Not used
415   //  
416   DGUARD(fDebug,1,"Process event in AliCentralMultiplicityTask");
417   fAODCentral.Clear("");
418   fIvz = 0;
419
420   AliESDEvent* esd = GetESDEvent();
421   if (!esd) return;
422
423   Bool_t   lowFlux   = kFALSE;
424   UInt_t   triggers  = 0;
425   UShort_t ivz       = 0;
426   TVector3 ip;
427   Double_t cent      = -1;
428   UShort_t nClusters = 0;
429   UInt_t   found     = fInspector.Process(esd, triggers, lowFlux, 
430                                           ivz, ip, cent, nClusters);
431
432   // No event or no trigger 
433   if (found & AliFMDEventInspector::kNoEvent)    return;
434   if (found & AliFMDEventInspector::kNoTriggers) return;
435   
436   // Make sure AOD is filled
437   MarkEventForStore();
438
439   if (found == AliFMDEventInspector::kNoSPD)      return;
440   if (found == AliFMDEventInspector::kNoVertex)   return;
441   if (triggers & AliAODForwardMult::kPileUp)      return;
442   if (found == AliFMDEventInspector::kBadVertex)  return; // Out of range
443   
444   //Doing analysis
445   fIvz = ivz;
446   const AliMultiplicity* spdmult = esd->GetMultiplicity();
447
448   TH2D& aodHist = fAODCentral.GetHistogram();
449
450   ProcessESD(aodHist, spdmult);
451   CorrectData(aodHist, ivz);
452   //Producing hit maps
453   // TList* hitList = static_cast<TList*>(fList->FindObject("hitMaps"));
454   TH2D* data  = static_cast<TH2D*>(fHits->At(ivz-1));
455   if(data) data->Add(&aodHist);
456   
457   PostData(1,fList);
458 }
459 //____________________________________________________________________
460 void 
461 AliCentralMultiplicityTask::ProcessESD(TH2D& aodHist, 
462                                        const AliMultiplicity* spdmult) const
463 {
464   DGUARD(fDebug,1,"Process the ESD in AliCentralMultiplicityTask");
465   fNTracklet->Reset();
466   fNCluster->Reset();
467
468   //Filling clusters in layer 1 used for tracklets...
469   for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++) {
470     Double_t eta = spdmult->GetEta(j);
471     fNTracklet->Fill(eta);
472     aodHist.Fill(eta,spdmult->GetPhi(j));
473   }
474
475   //...and then the unused ones in layer 1 
476   for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
477     Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
478     fNCluster->Fill(eta);
479     aodHist.Fill(eta, spdmult->GetPhiSingle(j));
480   }
481   fNClusterTracklet->Fill(fNCluster->GetEntries(), 
482                           fNTracklet->GetEntries());
483   
484   fNCluster->Divide(fNTracklet);
485   for (Int_t j = 1; j <= fNCluster->GetNbinsX(); j++)  
486     fClusterPerTracklet->Fill(fNCluster->GetXaxis()->GetBinCenter(j), 
487                               fNCluster->GetBinContent(j));
488
489 }
490
491 //____________________________________________________________________
492 void 
493 AliCentralMultiplicityTask::CorrectData(TH2D& aodHist, UShort_t vtxbin) const
494 {  
495   // Corrections
496   DGUARD(fDebug,1,"Correct data in AliCentralMultiplicityTask");
497   TH1D* hAcceptance = fManager.GetAcceptanceCorrection(vtxbin);
498   TH2D* hSecMap     = fManager.GetSecMapCorrection(vtxbin);
499   
500   if (!hSecMap)     AliFatal("No secondary map!");
501   if (!hAcceptance) AliFatal("No acceptance!");
502     
503   if (fUseSecondary && hSecMap) aodHist.Divide(hSecMap);
504
505   Int_t nY = aodHist.GetNbinsY();
506   for(Int_t ix = 1; ix <= aodHist.GetNbinsX(); ix++) {
507     Float_t accCor = hAcceptance->GetBinContent(ix);
508     Float_t accErr = hAcceptance->GetBinError(ix);
509
510     Bool_t fiducial = true;
511     if (ix < fEtaMin[vtxbin-1] || ix > fEtaMax[vtxbin-1]) 
512       fiducial = false;
513     //  Bool_t etabinSeen = kFALSE;  
514     for(Int_t iy = 1; iy <= nY; iy++) {
515 #if 1
516       if (!fiducial) { 
517         aodHist.SetBinContent(ix, iy, 0);
518         aodHist.SetBinError(ix, iy, 0);
519         continue;
520       }
521 #endif  
522       // Get currrent value 
523       Float_t aodValue = aodHist.GetBinContent(ix,iy);
524       Float_t aodErr   = aodHist.GetBinError(ix,iy);
525
526 #if 0 // This is done once in the FindEtaBins function
527       // Set underflow bin
528       Float_t secCor   = 0;
529       if(hSecMap)       secCor     = hSecMap->GetBinContent(ix,iy);
530       if (secCor > 0.5) etabinSeen = kTRUE;
531 #endif
532       if (aodValue < 0.000001) { 
533         aodHist.SetBinContent(ix,iy, 0); 
534         aodHist.SetBinError(ix,iy, 0); 
535         continue; 
536       }
537       if (!fUseAcceptance) continue; 
538
539       // Acceptance correction 
540       if (accCor   < 0.000001) accCor = 1;
541       Float_t aodNew   = aodValue / accCor ;
542       Float_t error    = aodNew*TMath::Sqrt(TMath::Power(aodErr/aodValue,2) +
543                                             TMath::Power(accErr/accCor,2) );
544       aodHist.SetBinContent(ix,iy, aodNew);
545       //test
546       aodHist.SetBinError(ix,iy,error);
547       aodHist.SetBinError(ix,iy,aodErr);
548     }
549     //Filling underflow bin if we eta bin is in range
550     if (fiducial) {
551       aodHist.SetBinContent(ix,0, 1.);
552       aodHist.SetBinContent(ix,nY+1, 1.);
553     }
554     // if (etabinSeen) aodHist.SetBinContent(ix,0, 1.);
555   }  
556 }
557
558 //____________________________________________________________________
559 void AliCentralMultiplicityTask::Terminate(Option_t* /*option*/) 
560 {
561   // 
562   // End of job
563   // 
564   // Parameters:
565   //    option Not used 
566   //
567   DGUARD(fDebug,1,"Process merged output in AliCentralMultiplicityTask");
568 }
569 //____________________________________________________________________
570 void
571 AliCentralMultiplicityTask::Print(Option_t* option) const
572 {
573   // 
574   // Print information 
575   // 
576   // Parameters:
577   //    option Not used
578   //
579   std::cout << ClassName() << ": " << GetName() << "\n" 
580             << std::boolalpha
581             << "  Use secondary correction:  " << fUseSecondary << '\n'
582             << "  Use acceptance correction: " << fUseAcceptance << '\n' 
583             << "  Off-line trigger mask:  0x" 
584             << std::hex     << std::setfill('0') 
585             << std::setw (8) << fOfflineTriggerMask 
586             << std::dec     << std::setfill (' ') 
587             << std::noboolalpha << std::endl;
588   AliCentralCorrSecondaryMap* secMap = GetManager().GetSecMap();
589   if (secMap) {
590     const TAxis& vaxis = secMap->GetVertexAxis();
591     std::cout << "  Eta ranges:\n"
592             << "     Vertex        | Eta bins\n"
593               << "   bin     range   | \n"
594               << "   ----------------+-----------" << std::endl;
595     for (Int_t v = 1; v <= vaxis.GetNbins(); v++) { 
596       std::cout << "   " << std::setw(2) << v << "  " 
597                 << std::setw(5) << vaxis.GetBinLowEdge(v) << "-"
598                 << std::setw(5) << vaxis.GetBinUpEdge(v) << " | ";
599       if (fEtaMin.GetSize() <= 0) 
600         std::cout << " ? -   ?";
601       else 
602         std::cout << std::setw(3) << fEtaMin[v-1] << "-" 
603                   << std::setw(3) << fEtaMax[v-1];
604       std::cout << std::endl;
605     }
606   }
607
608   gROOT->IncreaseDirLevel();
609   fManager.Print(option);
610   fInspector.Print(option);
611   gROOT->DecreaseDirLevel();
612   
613 }
614 //====================================================================
615 AliCentralMultiplicityTask::Manager::Manager() :
616   fAcceptancePath("$ALICE_ROOT/PWGLF/FORWARD/corrections/CentralAcceptance"),
617   fSecMapPath("$ALICE_ROOT/PWGLF/FORWARD/corrections/CentralSecMap"),
618   fAcceptance(),
619   fSecmap(),
620   fAcceptanceName("centralacceptance"),
621   fSecMapName("centralsecmap"),
622   fIsInit(kFALSE)
623 {
624   //
625   // Constructor 
626   // 
627 }
628 //____________________________________________________________________
629 AliCentralMultiplicityTask::Manager::Manager(const Manager& o) 
630   :fAcceptancePath(o.fAcceptancePath),
631    fSecMapPath(o.fSecMapPath),
632    fAcceptance(o.fAcceptance),
633    fSecmap(o.fSecmap),
634    fAcceptanceName(o.fAcceptanceName),
635    fSecMapName(o.fSecMapName),
636    fIsInit(o.fIsInit)
637 {
638   //
639   // Copy Constructor 
640   // 
641 }
642 //____________________________________________________________________
643 AliCentralMultiplicityTask::Manager&
644 AliCentralMultiplicityTask::Manager::operator=(const Manager& o)
645 {
646   //
647   // Assignment operator  
648   // 
649   if (&o == this) return *this; 
650   fAcceptancePath = o.fAcceptancePath;
651   fSecMapPath     = o.fSecMapPath;
652   fAcceptance     = o.fAcceptance;
653   fSecmap         = o.fSecmap;
654   fAcceptanceName = o.fAcceptanceName;
655   fSecMapName     = o.fSecMapName;
656   fIsInit         = o.fIsInit;
657   return *this;
658 }
659
660 //____________________________________________________________________
661 const char* 
662 AliCentralMultiplicityTask::Manager::GetFullFileName(UShort_t what, 
663                                                      UShort_t sys, 
664                                                      UShort_t sNN,  
665                                                      Short_t  field) const
666 {
667   // 
668   // Get full path name to object file 
669   // 
670   // Parameters:
671   //    what   What to get 
672   //    sys    Collision system
673   //    sNN    Center of mass energy 
674   //    field  Magnetic field 
675   // 
676   // Return:
677   //    
678   //
679   return Form("%s/%s",
680               what == 0 ? GetSecMapPath() : GetAcceptancePath(), 
681               GetFileName(what, sys, sNN, field));
682 }
683
684 //____________________________________________________________________
685 const char* 
686 AliCentralMultiplicityTask::Manager::GetFileName(UShort_t  what ,
687                                                  UShort_t  sys, 
688                                                  UShort_t  sNN,
689                                                  Short_t   field) const
690 {
691   // 
692   // Get the full path name 
693   // 
694   // Parameters:
695   //    what   What to get
696   //    sys    Collision system
697   //    sNN    Center of mass energy 
698   //    field  Magnetic field 
699   // 
700   // Return:
701   //    
702   //
703   // Must be static - otherwise the data may disappear on return from
704   // this member function
705   static TString fname = "";
706   
707   switch(what) {
708   case 0:  fname = fSecMapName;     break;
709   case 1:  fname = fAcceptanceName; break;
710   default:
711     ::Error("GetFileName", 
712             "Invalid indentifier %d for central object, must be 0 or 1!", what);
713     break;
714   }
715   fname.Append(Form("_%s_%04dGeV_%c%1dkG.root",
716                     AliForwardUtil::CollisionSystemString(sys), 
717                     sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field)));
718   
719   return fname.Data();
720 }
721
722 //____________________________________________________________________
723 TH2D* 
724 AliCentralMultiplicityTask::Manager::GetSecMapCorrection(UShort_t vtxbin) const
725 {
726   // 
727   // Get the secondary map
728   // 
729   // Parameters:
730   //    vtxbin 
731   // 
732   // Return:
733   //    
734   //
735   if (!fSecmap) { 
736     ::Warning("GetSecMapCorrection","No secondary map defined");
737     return 0;
738   }
739   return fSecmap->GetCorrection(vtxbin);
740 }
741 //____________________________________________________________________
742 TH1D* 
743 AliCentralMultiplicityTask::Manager::GetAcceptanceCorrection(UShort_t vtxbin) 
744   const 
745 {
746   // 
747   // Get the acceptance correction 
748   // 
749   // Parameters:
750   //    vtxbin 
751   // 
752   // Return:
753   //    
754   //
755   if (!fAcceptance) { 
756     ::Warning("GetAcceptanceCorrection","No acceptance map defined");
757     return 0;
758   }
759   return fAcceptance->GetCorrection(vtxbin);
760 }
761
762 //____________________________________________________________________
763 void 
764 AliCentralMultiplicityTask::Manager::Init(UShort_t  sys, 
765                                           UShort_t  sNN,
766                                           Short_t   field) 
767 {
768   // 
769   // Initialize 
770   // 
771   // Parameters:
772   //    sys    Collision system (1: pp, 2: PbPb, 3: pPb)
773   //    sNN    Center of mass energy per nucleon pair [GeV]
774   //    field  Magnetic field [kG]
775   //
776   if(fIsInit) ::Warning("Init","Already initialised - overriding...");
777   
778   TFile fsec(GetFullFileName(0,sys,sNN,field));
779   fSecmap = 
780     dynamic_cast<AliCentralCorrSecondaryMap*>(fsec.Get(fSecMapName.Data()));  
781   if(!fSecmap) {
782     ::Error("Init", "no central Secondary Map found!") ;
783     return;
784   }
785   TFile facc(GetFullFileName(1,sys,sNN,field));
786   fAcceptance = 
787     dynamic_cast<AliCentralCorrAcceptance*>(facc.Get(fAcceptanceName.Data()));
788   if(!fAcceptance) {
789     ::Error("Init", "no central Acceptance found!") ;
790     return;
791   }
792   
793   if(fSecmap && fAcceptance) {
794     fIsInit = kTRUE;
795     ::Info("Init", 
796            "Central Manager initialised for %s, energy %dGeV, field %dkG",
797            sys == 1 ? "pp" : sys == 2 ? "PbPb" : sys == 3 ? "pPb" : "unknown", sNN,field);
798   }  
799 }
800 //____________________________________________________________________
801 Bool_t
802 AliCentralMultiplicityTask::Manager::WriteFile(UShort_t what, 
803                                                UShort_t sys, 
804                                                UShort_t sNN, 
805                                                Short_t  fld, 
806                                                TObject* obj, 
807                                                Bool_t   full) const
808 {
809   // 
810   // Write correction output to (a temporary) file 
811   // 
812   // Parameters: 
813   //   What     What to write 
814   //   sys      Collision system (1: pp, 2: PbPb, 3: pPb)
815   //   sNN      Center of mass energy per nucleon (GeV)
816   //   fld      Field (kG)
817   //   obj      Object to write 
818   //   full     if true, write to full path, otherwise locally
819   // 
820   // Return: 
821   //   true on success. 
822   TString ofName;
823   if (!full)
824     ofName = GetFileName(what, sys, sNN, fld);
825   else 
826     ofName = GetFullFileName(what, sys, sNN, fld);
827   if (ofName.IsNull()) { 
828     AliErrorGeneral("Manager",Form("Unknown object type %d", what));
829     return false;
830   }
831   TFile* output = TFile::Open(ofName, "RECREATE");
832   if (!output) { 
833     AliErrorGeneral("Manager",Form("Failed to open file %s", ofName.Data()));
834     return false;
835   }
836   
837   TString oName(GetObjectName(what));
838   Int_t ret = obj->Write(oName);
839   if (ret <= 0) { 
840     AliErrorGeneral("Manager",Form("Failed to write %p to %s/%s (%d)", 
841                                    obj, ofName.Data(), oName.Data(), ret));
842     return false;
843   }
844
845   ret = output->Write();
846   if (ret < 0) { 
847     AliErrorGeneral("Manager",
848                     Form("Failed to write %s to disk (%d)", ofName.Data(),ret));
849     return false;
850   }
851   // output->ls();
852   output->Close();
853   
854 #if 0
855   TString cName(obj->IsA()->GetName());
856   AliInfoGeneral("Manager",
857                  Form("Wrote %s object %s to %s\n",
858                       cName.Data(),oName.Data(), ofName.Data()));
859   if (!full) { 
860     TString dName(GetFileDir(what));
861     AliInfoGeneral("Manager",
862                    Form("\n  %s should be copied to %s\n"
863                         "Do for example\n\t"
864                         "aliroot $ALICE_ROOT/PWGLF/FORWARD/analysis2/scripts/"
865                         "MoveCorrections.C\\(%d\\)\nor\n\t"
866                         "cp %s %s/", 
867                         ofName.Data(),dName.Data(), 
868                         what, ofName.Data(), 
869                         gSystem->ExpandPathName(dName.Data())));
870
871
872   }
873 #endif
874   return true;
875 }
876
877 //____________________________________________________________________
878 void 
879 AliCentralMultiplicityTask::Manager::Print(Option_t* option) const
880 {
881   // 
882   // Print information to standard output 
883   //
884   std::cout << " AliCentralMultiplicityTask::Manager\n" 
885             << std::boolalpha 
886             << "  Initialized:     " << fIsInit << '\n'
887             << "  Acceptance path: " << fAcceptancePath << '\n'
888             << "  Acceptance name: " << fAcceptanceName << '\n'
889             << "  Acceptance:      " << fAcceptance << '\n'
890             << "  Secondary path:  " << fSecMapPath << '\n'
891             << "  Secondary name:  " << fSecMapName << '\n'
892             << "  Secondary map:   " << fSecmap 
893             << std::noboolalpha << std::endl;
894   if (fAcceptance) fAcceptance->Print(option);
895   if (fSecmap)     fSecmap->Print(option);
896 }
897
898 //
899 // EOF
900 //