1 //====================================================================
3 // Base class for classes that calculate the multiplicity in the
4 // central region event-by-event
10 // - AliAODCentralMult
15 #include "AliCentralMultiplicityTask.h"
16 #include "AliAODForwardMult.h"
17 #include "AliForwardUtil.h"
19 #include "AliAODHandler.h"
20 #include "AliAnalysisManager.h"
21 #include "AliESDEvent.h"
22 #include "AliMultiplicity.h"
30 //====================================================================
31 AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name)
32 : AliAnalysisTaskSE(name),
33 fInspector("centralEventInspector"),
40 fFirstEventSeen(false),
43 fClusterPerTracklet(0),
52 DefineOutput(1, TList::Class());
54 "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
55 "SPDVertex.,PrimaryVertex.";
57 //____________________________________________________________________
58 AliCentralMultiplicityTask::AliCentralMultiplicityTask()
59 : AliAnalysisTaskSE(),
67 fFirstEventSeen(false),
70 fClusterPerTracklet(0),
80 //____________________________________________________________________
81 AliCentralMultiplicityTask::AliCentralMultiplicityTask(const AliCentralMultiplicityTask& o)
82 : AliAnalysisTaskSE(o),
83 fInspector(o.fInspector),
86 fAODCentral(o.fAODCentral),
88 fUseSecondary(o.fUseSecondary),
89 fUseAcceptance(o.fUseAcceptance),
90 fFirstEventSeen(o.fFirstEventSeen),
92 fNClusterTracklet(o.fNClusterTracklet),
93 fClusterPerTracklet(o.fClusterPerTracklet),
94 fNCluster(o.fNCluster),
95 fNTracklet(o.fNTracklet),
103 //____________________________________________________________________
104 AliCentralMultiplicityTask&
105 AliCentralMultiplicityTask::operator=(const AliCentralMultiplicityTask& o)
108 // Assignment operator
110 if (&o == this) return *this;
111 fInspector = o.fInspector;
114 fAODCentral = o.fAODCentral;
115 fManager = o.fManager;
116 fUseSecondary = o.fUseSecondary;
117 fUseAcceptance = o.fUseAcceptance;
118 fFirstEventSeen = o.fFirstEventSeen;
120 fNClusterTracklet = o.fNClusterTracklet;
121 fClusterPerTracklet= o.fClusterPerTracklet;
122 fNCluster = o.fNCluster;
123 fNTracklet = o.fNTracklet;
128 //____________________________________________________________________
129 void AliCentralMultiplicityTask::UserCreateOutputObjects()
132 // Create output objects
136 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
138 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
139 if (!ah) AliFatal("No AOD output handler set in analysis manager");
142 TObject* obj = &fAODCentral;
143 ah->AddBranch("AliAODCentralMult", &obj);
150 fInspector.DefineOutput(fList);
155 //____________________________________________________________________
157 AliCentralMultiplicityTask::GetESDEvent()
160 // Get the ESD event. IF this is the first event, initialise
162 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
164 AliWarning("No ESD event found for input event");
168 // IF we've read the first event already, just return the event
169 if (fFirstEventSeen) return esd;
171 // Read the details of the rung
172 fInspector.ReadRunDetails(esd);
174 // If we weren't initialised before (i.e., in the setup), do so now.
175 if (!GetManager().IsInit()) {
176 GetManager().Init(fInspector.GetCollisionSystem(),
177 fInspector.GetEnergy(),
178 fInspector.GetField());
179 AliInfo("Manager of corrections in AliCentralMultiplicityTask init");
182 // Check for existence and get secondary map
183 AliCentralCorrSecondaryMap* secMap = GetManager().GetSecMap();
184 if (!secMap) AliFatal("No secondary map defined!");
185 const TAxis& vaxis = secMap->GetVertexAxis();
189 fNClusterTracklet = new TH2D("nClusterVsnTracklet",
190 "Total number of cluster vs number of tracklets",
191 100, 0, 100, 100, 0, 100);
192 fNClusterTracklet->SetDirectory(0);
193 fNClusterTracklet->SetXTitle("# of free clusters");
194 fNClusterTracklet->SetYTitle("# of tracklets");
195 fNClusterTracklet->SetStats(0);
196 fList->Add(fNClusterTracklet);
200 fClusterPerTracklet = new TH2D("clusterPerTracklet",
201 "N_{free cluster}/N_{tracklet} vs. #eta",
202 nEta,-lEta,lEta, 101, -.05, 10.05);
203 fClusterPerTracklet->SetDirectory(0);
204 fClusterPerTracklet->SetXTitle("#eta");
205 fClusterPerTracklet->SetYTitle("N_{free cluster}/N_{tracklet}");
206 fClusterPerTracklet->SetStats(0);
207 fList->Add(fClusterPerTracklet);
210 fNCluster = new TH1D("cacheCluster", "", nEta,-lEta,lEta);
211 fNCluster->SetDirectory(0);
214 fNTracklet = new TH1D("cacheTracklet", "", nEta,-lEta,lEta);
215 fNTracklet->SetDirectory(0);
218 // Initialize the inspecto
219 fInspector.Init(vaxis);
220 fFirstEventSeen = kTRUE;
222 // Print some information
227 //____________________________________________________________________
229 AliCentralMultiplicityTask::MarkEventForStore() const
231 // Make sure the AOD tree is filled
232 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
234 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
236 AliFatal("No AOD output handler set in analysis manager");
238 ah->SetFillAOD(kTRUE);
241 //____________________________________________________________________
242 void AliCentralMultiplicityTask::FindEtaLimits()
244 AliCentralCorrSecondaryMap* secMap = GetManager().GetSecMap();
246 const TAxis& vaxis = secMap->GetVertexAxis();
248 fEtaMin.Set(vaxis.GetNbins());
249 fEtaMax.Set(vaxis.GetNbins());
251 TList* hits = new TList;
253 hits->SetName("hitMaps");
256 TList* secs = new TList;
258 secs->SetName("secondaryMaps");
260 TH2D* hCoverage = new TH2D("coverage", "#eta coverage per v_{z}",
261 secMap->GetCorrection(UShort_t(1))->GetXaxis()->GetNbins(),
262 secMap->GetCorrection(UShort_t(1))->GetXaxis()->GetXmin(),
263 secMap->GetCorrection(UShort_t(1))->GetXaxis()->GetXmax(),
264 vaxis.GetNbins(),vaxis.GetXmin(),vaxis.GetXmax());
265 hCoverage->SetDirectory(0);
266 hCoverage->SetXTitle("#eta");
267 hCoverage->SetYTitle("v_{z} [cm]");
268 hCoverage->SetZTitle("n_{bins}");
269 fList->Add(hCoverage);
271 for (Int_t v = 1; v <= vaxis.GetNbins(); v++) {
272 TH2D* corr = secMap->GetCorrection(UShort_t(v));
273 TH1D* proj = corr->ProjectionX(Form("secCor%02d", v));
274 proj->Scale(1. / corr->GetNbinsY());
275 proj->SetTitle(Form("Projection of secondary correction "
276 "for %+5.1f<v_{z}<%+5.1f",
277 vaxis.GetBinLowEdge(v), vaxis.GetBinUpEdge(v)));
278 proj->SetYTitle("#LT 2^{nd} correction#GT");
279 proj->SetDirectory(0);
280 proj->SetMarkerStyle(20);
281 proj->SetMarkerColor(kBlue+1);
284 TH2D* obg = static_cast<TH2D*>(corr->Clone(Form("secCor2DFiducial%02d",v)));
285 obg->SetDirectory(0);
288 TH1D* after = static_cast<TH1D*>(proj->Clone(Form("secCorFiducial%02d",v)));
289 after->SetDirectory(0);
290 after->SetMarkerColor(kRed+1);
293 TH2D* data = static_cast<TH2D*>(corr->Clone(Form("hitMap%02d",v)));
294 //d->SetTitle(Form("hitMap%02d",v));
295 data->SetTitle(Form("d^{2}N/d#eta d#phi "
296 "for %+5.1f<v_{z}<%+5.1f",
297 vaxis.GetBinLowEdge(v), vaxis.GetBinUpEdge(v)));
298 data->GetZaxis()->SetTitle("");
299 data->SetMarkerColor(kBlack);
300 data->SetMarkerStyle(1);
303 TH1D* hAcceptance = fManager.GetAcceptanceCorrection(v);
304 TH1D* accClone = static_cast<TH1D*>(hAcceptance->Clone(Form("acceptance%02d",v)));
307 // Double_t prev = 0;
308 for (Int_t e = 1; e <= proj->GetNbinsX(); e++) {
309 Double_t c = proj->GetBinContent(e);
310 if (c > .5 /*&& TMath::Abs(c - prev) < .1*c*/) {
315 after->SetBinContent(e, 0);
316 after->SetBinError(e, 0);
317 for(Int_t nn =1; nn <=obg->GetNbinsY();nn++)
318 obg->SetBinContent(e,nn,0);
321 for (Int_t e = proj->GetNbinsX(); e >= 1; e--) {
322 Double_t c = proj->GetBinContent(e);
323 if (c > .5 /*&& TMath::Abs(c - prev) < .1*c*/) {
328 after->SetBinContent(e, 0);
329 after->SetBinError(e, 0);
330 for(Int_t nn =1; nn <=obg->GetNbinsY();nn++)
331 obg->SetBinContent(e,nn,0);
335 for (Int_t nn = fEtaMin[v-1]; nn<=fEtaMax[v-1]; nn++) {
336 hCoverage->SetBinContent(nn,v,1);
342 //____________________________________________________________________
343 void AliCentralMultiplicityTask::UserExec(Option_t* /*option*/)
346 // Process each event
351 fAODCentral.Clear("");
354 AliESDEvent* esd = GetESDEvent();
356 Bool_t lowFlux = kFALSE;
361 UShort_t nClusters = 0;
362 UInt_t found = fInspector.Process(esd, triggers, lowFlux,
363 ivz, vz, cent, nClusters);
365 // No event or no trigger
366 if (found & AliFMDEventInspector::kNoEvent) return;
367 if (found & AliFMDEventInspector::kNoTriggers) return;
369 // Make sure AOD is filled
372 if (found == AliFMDEventInspector::kNoSPD) return;
373 if (found == AliFMDEventInspector::kNoVertex) return;
374 if (triggers & AliAODForwardMult::kPileUp) return;
375 if (found == AliFMDEventInspector::kBadVertex) return; // Out of range
379 const AliMultiplicity* spdmult = esd->GetMultiplicity();
381 TH2D& aodHist = fAODCentral.GetHistogram();
383 ProcessESD(aodHist, spdmult);
384 CorrectData(aodHist, ivz);
386 TList* hitList = static_cast<TList*>(fList->FindObject("hitMaps"));
389 data = static_cast<TH2D*>(hitList->At(ivz-1));
395 //____________________________________________________________________
397 AliCentralMultiplicityTask::ProcessESD(TH2D& aodHist,
398 const AliMultiplicity* spdmult) const
403 //Filling clusters in layer 1 used for tracklets...
404 for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++) {
405 Double_t eta = spdmult->GetEta(j);
406 fNTracklet->Fill(eta);
407 aodHist.Fill(eta,spdmult->GetPhi(j));
410 //...and then the unused ones in layer 1
411 for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
412 Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
413 fNCluster->Fill(eta);
414 aodHist.Fill(eta, spdmult->GetPhiSingle(j));
416 fNClusterTracklet->Fill(fNCluster->GetEntries(),
417 fNTracklet->GetEntries());
419 fNCluster->Divide(fNTracklet);
420 for (Int_t j = 1; j <= fNCluster->GetNbinsX(); j++)
421 fClusterPerTracklet->Fill(fNCluster->GetXaxis()->GetBinCenter(j),
422 fNCluster->GetBinContent(j));
426 //____________________________________________________________________
428 AliCentralMultiplicityTask::CorrectData(TH2D& aodHist, UShort_t vtxbin) const
431 TH1D* hAcceptance = fManager.GetAcceptanceCorrection(vtxbin);
432 TH2D* hSecMap = fManager.GetSecMapCorrection(vtxbin);
434 if (!hSecMap) AliFatal("No secondary map!");
435 if (!hAcceptance) AliFatal("No acceptance!");
437 if (fUseSecondary && hSecMap) aodHist.Divide(hSecMap);
439 for(Int_t nx = 1; nx <= aodHist.GetNbinsX(); nx++) {
440 Float_t accCor = hAcceptance->GetBinContent(nx);
441 Float_t accErr = hAcceptance->GetBinError(nx);
443 Bool_t fiducial = true;
444 if (nx < fEtaMin[vtxbin-1] || nx > fEtaMax[vtxbin-1])
446 // Bool_t etabinSeen = kFALSE;
447 for(Int_t ny = 1; ny <= aodHist.GetNbinsY(); ny++) {
450 aodHist.SetBinContent(nx, ny, 0);
451 aodHist.SetBinError(nx, ny, 0);
455 // Get currrent value
456 Float_t aodValue = aodHist.GetBinContent(nx,ny);
457 Float_t aodErr = aodHist.GetBinError(nx,ny);
459 #if 0 // This is done once in the FindEtaBins function
462 if(hSecMap) secCor = hSecMap->GetBinContent(nx,ny);
463 if (secCor > 0.5) etabinSeen = kTRUE;
465 if (aodValue < 0.000001) {
466 aodHist.SetBinContent(nx,ny, 0);
467 aodHist.SetBinError(nx,ny, 0);
470 if (!fUseAcceptance) continue;
472 // Acceptance correction
473 if (accCor < 0.000001) accCor = 1;
474 Float_t aodNew = aodValue / accCor ;
475 Float_t error = aodNew*TMath::Sqrt(TMath::Power(aodErr/aodValue,2) +
476 TMath::Power(accErr/accCor,2) );
477 aodHist.SetBinContent(nx,ny, aodNew);
479 aodHist.SetBinError(nx,ny,error);
480 aodHist.SetBinError(nx,ny,aodErr);
482 //Filling underflow bin if we eta bin is in range
483 if (fiducial) aodHist.SetBinContent(nx,0, 1.);
484 // if (etabinSeen) aodHist.SetBinContent(nx,0, 1.);
488 //____________________________________________________________________
489 void AliCentralMultiplicityTask::Terminate(Option_t* /*option*/)
498 //____________________________________________________________________
500 AliCentralMultiplicityTask::Print(Option_t* option) const
508 std::cout << ClassName() << ": " << GetName() << "\n"
510 << " Use secondary correction: " << fUseSecondary << '\n'
511 << " Use acceptance correction: " << fUseAcceptance << '\n'
512 << " Off-line trigger mask: 0x"
513 << std::hex << std::setfill('0')
514 << std::setw (8) << fOfflineTriggerMask
515 << std::dec << std::setfill (' ')
516 << std::noboolalpha << std::endl;
517 AliCentralCorrSecondaryMap* secMap = GetManager().GetSecMap();
519 const TAxis& vaxis = secMap->GetVertexAxis();
520 std::cout << " Eta ranges:\n"
521 << " Vertex | Eta bins\n"
523 << " ----------------+-----------" << std::endl;
524 for (Int_t v = 1; v <= vaxis.GetNbins(); v++) {
525 std::cout << " " << std::setw(2) << v << " "
526 << std::setw(5) << vaxis.GetBinLowEdge(v) << "-"
527 << std::setw(5) << vaxis.GetBinUpEdge(v) << " | "
528 << std::setw(3) << fEtaMin[v-1] << "-"
529 << std::setw(3) << fEtaMax[v-1] << std::endl;
533 gROOT->IncreaseDirLevel();
534 fManager.Print(option);
535 fInspector.Print(option);
536 gROOT->DecreaseDirLevel();
539 //====================================================================
540 AliCentralMultiplicityTask::Manager::Manager() :
541 fAcceptancePath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralAcceptance"),
542 fSecMapPath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralSecMap"),
545 fAcceptanceName("centralacceptance"),
546 fSecMapName("centralsecmap"),
553 //____________________________________________________________________
554 AliCentralMultiplicityTask::Manager::Manager(const Manager& o)
555 :fAcceptancePath(o.fAcceptancePath),
556 fSecMapPath(o.fSecMapPath),
557 fAcceptance(o.fAcceptance),
559 fAcceptanceName(o.fAcceptanceName),
560 fSecMapName(o.fSecMapName),
567 //____________________________________________________________________
568 AliCentralMultiplicityTask::Manager&
569 AliCentralMultiplicityTask::Manager::operator=(const Manager& o)
572 // Assignment operator
574 if (&o == this) return *this;
575 fAcceptancePath = o.fAcceptancePath;
576 fSecMapPath = o.fSecMapPath;
577 fAcceptance = o.fAcceptance;
579 fAcceptanceName = o.fAcceptanceName;
580 fSecMapName = o.fSecMapName;
585 //____________________________________________________________________
587 AliCentralMultiplicityTask::Manager::GetFullFileName(UShort_t what,
593 // Get full path name to object file
597 // sys Collision system
598 // sNN Center of mass energy
599 // field Magnetic field
605 what == 0 ? GetSecMapPath() : GetAcceptancePath(),
606 GetFileName(what, sys, sNN, field));
609 //____________________________________________________________________
611 AliCentralMultiplicityTask::Manager::GetFileName(UShort_t what ,
617 // Get the full path name
621 // sys Collision system
622 // sNN Center of mass energy
623 // field Magnetic field
628 // Must be static - otherwise the data may disappear on return from
629 // this member function
630 static TString fname = "";
633 case 0: fname = fSecMapName; break;
634 case 1: fname = fAcceptanceName; break;
636 ::Error("GetFileName",
637 "Invalid indentifier %d for central object, must be 0 or 1!", what);
640 fname.Append(Form("_%s_%04dGeV_%c%1dkG.root",
641 AliForwardUtil::CollisionSystemString(sys),
642 sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field)));
647 //____________________________________________________________________
649 AliCentralMultiplicityTask::Manager::GetSecMapCorrection(UShort_t vtxbin) const
652 // Get the secondary map
661 ::Warning("GetSecMapCorrection","No secondary map defined");
664 return fSecmap->GetCorrection(vtxbin);
666 //____________________________________________________________________
668 AliCentralMultiplicityTask::Manager::GetAcceptanceCorrection(UShort_t vtxbin)
672 // Get the acceptance correction
681 ::Warning("GetAcceptanceCorrection","No acceptance map defined");
684 return fAcceptance->GetCorrection(vtxbin);
687 //____________________________________________________________________
689 AliCentralMultiplicityTask::Manager::Init(UShort_t sys,
697 // sys Collision system (1: pp, 2: PbPb, 3: pPb)
698 // sNN Center of mass energy per nucleon pair [GeV]
699 // field Magnetic field [kG]
701 if(fIsInit) ::Warning("Init","Already initialised - overriding...");
703 TFile fsec(GetFullFileName(0,sys,sNN,field));
705 dynamic_cast<AliCentralCorrSecondaryMap*>(fsec.Get(fSecMapName.Data()));
707 ::Error("Init", "no central Secondary Map found!") ;
710 TFile facc(GetFullFileName(1,sys,sNN,field));
712 dynamic_cast<AliCentralCorrAcceptance*>(facc.Get(fAcceptanceName.Data()));
714 ::Error("Init", "no central Acceptance found!") ;
718 if(fSecmap && fAcceptance) {
721 "Central Manager initialised for %s, energy %dGeV, field %dkG",
722 sys == 1 ? "pp" : sys == 2 ? "PbPb" : sys == 3 ? "pPb" : "unknown", sNN,field);
725 //____________________________________________________________________
727 AliCentralMultiplicityTask::Manager::WriteFile(UShort_t what,
735 // Write correction output to (a temporary) file
738 // What What to write
739 // sys Collision system (1: pp, 2: PbPb, 3: pPb)
740 // sNN Center of mass energy per nucleon (GeV)
742 // obj Object to write
743 // full if true, write to full path, otherwise locally
749 ofName = GetFileName(what, sys, sNN, fld);
751 ofName = GetFullFileName(what, sys, sNN, fld);
752 if (ofName.IsNull()) {
753 AliErrorGeneral("Manager",Form("Unknown object type %d", what));
756 TFile* output = TFile::Open(ofName, "RECREATE");
758 AliErrorGeneral("Manager",Form("Failed to open file %s", ofName.Data()));
762 TString oName(GetObjectName(what));
763 Int_t ret = obj->Write(oName);
765 AliErrorGeneral("Manager",Form("Failed to write %p to %s/%s (%d)",
766 obj, ofName.Data(), oName.Data(), ret));
770 ret = output->Write();
772 AliErrorGeneral("Manager",
773 Form("Failed to write %s to disk (%d)", ofName.Data(),ret));
779 TString cName(obj->IsA()->GetName());
780 AliInfoGeneral("Manager",
781 Form("Wrote %s object %s to %s\n",
782 cName.Data(),oName.Data(), ofName.Data()));
784 TString dName(GetFileDir(what));
785 AliInfoGeneral("Manager",
786 Form("%s should be copied to %s\n"
788 "aliroot $ALICE_ROOT/PWG2/FORWARD/analysis2/scripts/"
789 "MoveCorrections.C\\(%d\\)\nor\n\t"
791 ofName.Data(),dName.Data(),
793 gSystem->ExpandPathName(dName.Data())));
800 //____________________________________________________________________
802 AliCentralMultiplicityTask::Manager::Print(Option_t* option) const
805 // Print information to standard output
807 std::cout << " AliCentralMultiplicityTask::Manager\n"
809 << " Initialized: " << fIsInit << '\n'
810 << " Acceptance path: " << fAcceptancePath << '\n'
811 << " Acceptance name: " << fAcceptanceName << '\n'
812 << " Acceptance: " << fAcceptance << '\n'
813 << " Secondary path: " << fSecMapPath << '\n'
814 << " Secondary name: " << fSecMapName << '\n'
815 << " Secondary map: " << fSecmap
816 << std::noboolalpha << std::endl;
817 if (fAcceptance) fAcceptance->Print(option);
818 if (fSecmap) fSecmap->Print(option);