From: Christian Holm Christensen Date: Mon, 9 Dec 2013 13:37:22 +0000 (+0100) Subject: Added tasks and trains to investigate time-to-previous event. X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=9aba161af6b85454d292523356da2c57e5ba0c81;p=u%2Fmrichter%2FAliRoot.git Added tasks and trains to investigate time-to-previous event. Some fixes for the AOD and MultDist trains. --- diff --git a/PWGLF/FORWARD/analysis2/trains/ELossTimeTask.C b/PWGLF/FORWARD/analysis2/trains/ELossTimeTask.C new file mode 100644 index 00000000000..8a3f5b7a791 --- /dev/null +++ b/PWGLF/FORWARD/analysis2/trains/ELossTimeTask.C @@ -0,0 +1,372 @@ +#ifndef __CINT__ +# include +# include +# include +# include +# include +# include +# include +# include +#else +class AliFMDEventInspector; +class TH2; +class AliESDFMD; +#endif +#include +#include +#define NO_TASK +#define NO_SORTER +#include "EventTimeTask.C" + +/** + * Task to analyse the energy loss in the FMD rings as function of the + * time to previous event. + * + */ +struct ELossTimeTask : public AliBaseESDTask +{ + /** + * Constructor - for I/O only + */ + ELossTimeTask() + : AliBaseESDTask(), + fEventInspector(), + fFMD1i(), + fFMD2i(), + fFMD2o(), + fFMD3i(), + fFMD3o(), + fMap(0), + fDt(0) + {} + /** + * Constructor + * + * @param name Name of task + */ + ELossTimeTask(const char* name) + : AliBaseESDTask(name, "ELossTimeTask", + &(AliForwardCorrectionManager::Instance())), + fEventInspector("fmdEventInspector"), + fFMD1i(1, 'i'), + fFMD2i(2, 'i'), + fFMD2o(2, 'o'), + fFMD3i(3, 'i'), + fFMD3o(3, 'o'), + fMap(0), + fDt(0) + { + } + /** + * Book output objects. Derived class should define this to book + * output objects on the processing output list @c fList before the + * actual event processing. This is called on the master and on + * each slave. + * + * If this member function returns false, the execution is stopped + * with a fatal signal. + * + * @return true on success. + */ + virtual Bool_t Book() + { + fNeededCorrections = 0; + fExtraCorrections = 0; + + fDt = new TH1D("dt", "Time-to-last event (PS-triggered)", 60, 0, 15); + fDt->SetXTitle("log_{10}(#Deltat)"); + fDt->SetFillColor(kYellow+2); + fDt->SetLineColor(kYellow+2); + fDt->SetFillStyle(3001); + fDt->SetDirectory(0); + fList->Add(fDt); + + fFMD1i.Book(fList, fDt); + fFMD2i.Book(fList, fDt); + fFMD2o.Book(fList, fDt); + fFMD3i.Book(fList, fDt); + fFMD3o.Book(fList, fDt); + + // Possibly re-read map + ReadMap("map.root"); + + return true; + } + /** + * Process a single event + * + * @param esd Input event + * + * @return true on success + */ + virtual Bool_t Event(AliESDEvent& esd) + { + Bool_t lowFlux = kFALSE; + UInt_t triggers = 0; + UShort_t ivz = 0; + TVector3 ip; + Double_t cent = 0; + UShort_t nClusters = 0; + UInt_t found = fEventInspector.Process(&esd, triggers, lowFlux, + ivz, ip, cent, nClusters); + if (found & AliFMDEventInspector::kNoEvent) return false; + if (found & AliFMDEventInspector::kNoTriggers) return false; + if (found & AliFMDEventInspector::kNoSPD) return false; + if (found & AliFMDEventInspector::kNoFMD) return false; + if (found & AliFMDEventInspector::kNoVertex) return false; + if (found & AliFMDEventInspector::kBadVertex) return false; + + // do not process pile-up, A, C, and E events + if (triggers & AliAODForwardMult::kPileUp) return false; + if (triggers & AliAODForwardMult::kA) return false; + if (triggers & AliAODForwardMult::kC) return false; + if (triggers & AliAODForwardMult::kE) return false; + + // We want only the events found by off-line + if (!(triggers & AliAODForwardMult::kOffline)) return false; + + // Perhaps we should also insist on MB only + // if (fOnlyMB && (!(triggers & AliAODForwardMult::kInel))) return false; + + // Get FMD data + AliESDFMD* esdFMD = esd.GetFMDData(); + ULong64_t period = esd.GetPeriodNumber(); + ULong64_t orbit = esd.GetOrbitNumber(); + ULong64_t bc = esd.GetBunchCrossNumber(); + ULong64_t full = EventTimeData::EncodeFull(period, orbit, bc); + ULong64_t dt = fMap->Get(full); + Double_t logDt = TMath::Log10(25. * dt); + if (dt == EventTimeMap::kInvalidTime) { + logDt = 0; + Printf("!!! Failed to find dT for 0x%016llu", full); + } + // else + // Printf("=== 0x%016llu -> 0x%016llu", full, dt); + fDt->Fill(logDt); + + for (UShort_t d = 1; d <= 3; d++) { + UShort_t nQ = d == 1 ? 1 : 2; + for (UShort_t q = 0; q < nQ; q++) { + RingHistos* r = 0; + switch (d) { + case 1: r = &fFMD1i; break; + case 2: r = (q == 0 ? &fFMD2i : &fFMD2o); break; + case 3: r = (q == 0 ? &fFMD3i : &fFMD3o); break; + } + r->Event(*esdFMD, logDt); + } + } + return true; + } + /** + * Do the final analysis on the merged output. + * + * @return true on success + */ + virtual Bool_t Finalize() + { + fDt = static_cast(fList->FindObject("dt")); + + fFMD1i.Finalize(fList, fResults, fDt); + fFMD2i.Finalize(fList, fResults, fDt); + fFMD2o.Finalize(fList, fResults, fDt); + fFMD3i.Finalize(fList, fResults, fDt); + fFMD3o.Finalize(fList, fResults, fDt); + return true; + } + /** + * Get a reference to the event inspector. User must override this + * to return proper object + * + * @return Reference to the event inspector + */ + virtual AliFMDEventInspector& GetEventInspector() { return fEventInspector; } + /** + * Get a reference to the event inspector. User must override this + * to return proper object + * + * @return Reference to the event inspector + */ + virtual const AliFMDEventInspector& GetEventInspector() const + { + return fEventInspector; + } + /** + * Read the map from timestamp to time-to-previous event + * + * @param filename File to read the map from + * + * @return true on success, false otherwise + */ + Bool_t ReadMap(const char* filename) + { + if (gSystem->AccessPathName(filename, kReadPermission)) { + // TSystem::AccessPathName returns false if file can be accessed! + Error("ReadMap", "File \"%s\" cannot be open for reading", filename); + return false; + } + Printf("Opening \"%s\" ...", filename); + TFile* file = TFile::Open(filename, "READ"); + if (!file) { + Error("ReadMap", "Failed to open file \"%s\"", filename); + return false; + } + Printf("Opened \"%s\" ...", filename); + TObject* o = file->Get("eventTimeMap"); + if (!o) { + Error("ReadMap", "Failed to get \"eventTimeMap\" from %s", filename); + return false; + } + Printf("Got object \"eventTimeMap\" from \"%s\" ...", filename); + if (!o->IsA()->InheritsFrom(EventTimeMap::Class())) { + Error("ReadMap", "Object \"%s\" is not an EventTimeMap, but a %s", + o->GetName(), o->ClassName()); + return false; + } + Printf("Set the \"eventTimeMap\" to use"); + if (fMap) { + delete fMap; + fMap = 0; + } + fMap = static_cast(o); + file->Close(); + return true; + } + /** + * Create and connect the task + * + * @param mapfile File name of file containing timestamp map + * + * @return true on connect + */ + static Bool_t Create(const char* mapfile) + { + ELossTimeTask* task = new ELossTimeTask("elossTime"); + if (!task->ReadMap(mapfile)) return false; + task->Connect(); + return true; + } +protected: + /** + * Dummy copy constructor + * + * @param o Object to copy from + */ + ELossTimeTask(const ELossTimeTask& o) : AliBaseESDTask(o) {} + /** Our event inspector */ + AliFMDEventInspector fEventInspector; + /** + * Structure to hold per-ring histograms + */ + struct RingHistos : public AliForwardUtil::RingHistos + { + /** + * Constructor - for I/O only + */ + RingHistos() + : AliForwardUtil::RingHistos(), + fDtVsELoss(0) + {} + /** + * Constructor + * + * @param d detector number + * @param r ring identifier + */ + RingHistos(UShort_t d, Char_t r) + : AliForwardUtil::RingHistos(d,r), + fDtVsELoss(0) + {} + /** + * Book histograms + * + * @param dir Parent list to add our list to + * + * @return true on success + */ + Bool_t Book(TList* dir, TH1* dt) + { + TList* l = DefineOutputList(dir); + // Double_t dtBins[] = { 0, 1, 2e3, 5e5, 1.1e6, 5e6, 1e12, 1e20 }; + fDtVsELoss = new TH2D("dtVsELoss", + Form("#Deltat vs #Delta/#Delta_{mip} - %s", + GetName()), 450, 0, 15, + dt->GetXaxis()->GetNbins(), + dt->GetXaxis()->GetXmin(), + dt->GetXaxis()->GetXmax()); + fDtVsELoss->SetXTitle("#Delta/#Delta_{mip}"); + fDtVsELoss->SetYTitle("log_{10}(#Deltat) [ns]"); + fDtVsELoss->SetMarkerColor(Color()); + fDtVsELoss->SetDirectory(0); + l->Add(fDtVsELoss); + + return true; + } + /** + * Process a single event + * + * @param fmd FMD ESD data + * @param dt Time-to-previous event + */ + void Event(AliESDFMD& fmd, Double_t logDt) + { + const UShort_t nSec = NSector(); + const UShort_t nStr = NStrip(); + for (UShort_t s = 0; s < nSec; s++) { + for (UShort_t t = 0; t < nStr; t++) { + Float_t mult = fmd.Multiplicity(fDet,fRing,s,t); + if(mult == AliESDFMD::kInvalidMult || mult <= 0) + continue; + fDtVsELoss->Fill(mult, logDt); + } + } + } + void Finalize(const TList* input, TList* output, TH1* dt) + { + TList* in = GetOutputList(input); + TList* out = DefineOutputList(output); + + fDtVsELoss = static_cast(in->FindObject("dtVsELoss")); + + TH2* dtVsELoss = static_cast(fDtVsELoss->Clone()); + dtVsELoss->SetDirectory(0); + dtVsELoss->SetZTitle("1/N_{events}"); + dtVsELoss->Reset(); + for (UShort_t j = 1; j <= dt->GetNbinsX(); j++) { + Double_t norm = dt->GetBinContent(j); + if (norm <= 0) { + // Warning("Finalize", "Bin %d empty in dT", j); + continue; + } + for (UShort_t i = 1; i <= fDtVsELoss->GetNbinsX(); i++) { + Double_t c = fDtVsELoss->GetBinContent(i, j); + Double_t e = fDtVsELoss->GetBinError(i, j); + dtVsELoss->SetBinContent(i,j,c/norm); + dtVsELoss->SetBinError(i,j,e/norm); + } + } + out->Add(dtVsELoss); + } + /** Our histogram */ + TH2* fDtVsELoss; + ClassDef(RingHistos,1) + }; + /** Container of FMD1i histograms */ + RingHistos fFMD1i; + /** Container of FMD2i histograms */ + RingHistos fFMD2i; + /** Container of FMD2o histograms */ + RingHistos fFMD2o; + /** Container of FMD3i histograms */ + RingHistos fFMD3i; + /** Container of FMD3o histograms */ + RingHistos fFMD3o; + /** Map from timestamp to time-to-previous event*/ + EventTimeMap* fMap; + /** Distribution of log10(dt) */ + TH1* fDt; + ClassDef(ELossTimeTask,2) +}; +// +// EOF +// diff --git a/PWGLF/FORWARD/analysis2/trains/ELossTimeTrain.C b/PWGLF/FORWARD/analysis2/trains/ELossTimeTrain.C new file mode 100644 index 00000000000..2bfc4672bb0 --- /dev/null +++ b/PWGLF/FORWARD/analysis2/trains/ELossTimeTrain.C @@ -0,0 +1,291 @@ +// ELossTimeTrain.C +#ifndef NO_TRAIN +#ifndef __CINT__ +#include +#include +#else +class AliAnalysisManager; +#endif +#include "TrainSetup.C" +#include "ParUtilities.C" + +/** + * Train to record time of each event + * + */ +class ELossTimeTrain : public TrainSetup +{ +public: + /** + * Constructor + * + * @param name The name + */ + ELossTimeTrain(const char* name="eventTime") : TrainSetup(name) + { + fOptions.Add("map", "FILE", "File containg map", "map.root"); + fOptions.Set("type", "ESD"); + } + /** + * Create our tasks + * + * @param mgr Manager + */ + void CreateTasks(AliAnalysisManager*) + { + if (!fHelper->LoadLibrary("PWGLFforward2")) + Fatal("CreateTasks", "Failed to load PWGLFforward2"); + + if (!ParUtilities::MakeScriptPAR(fHelper->Mode() == Helper::kLocal, + "EventTimeTask.C", + // Gui because of CDB - sigh! + // XMLParser because of CDB + "Gui,XMLParser," + "STEERBase,CDB,ESD,AOD,ANALYSIS,OADB," + "ANALYSISalice", + fHelper)) + Fatal("","Failed to make support PAR"); + if (!fHelper->LoadLibrary("EventTimeTask")) + Fatal("CreateTasks", "Failed to load EventTimeTask"); + + if (!ParUtilities::MakeScriptPAR(fHelper->Mode() == Helper::kLocal, + "ELossTimeTask.C", + "Gui,STEERBase,CDB,ESD,AOD,ANALYSIS,OADB," + "ANALYSISalice,PWGLFforward2," + "EventTimeTask", + fHelper)) + Fatal("","Failed to make PAR"); + if (!fHelper->LoadLibrary("ELossTimeTask")) + Fatal("CreateTasks", "Failed to load ELossTimeTask"); + + TString mapfile = fOptions.Get("map"); + gROOT->ProcessLine(Form("ELossTimeTask::Create(\"%s\")", mapfile.Data())); + + fHelper->LoadAux(mapfile.Data(), true); + } + /** + * Do not create a physics selection + */ + // void CreatePhysicsSelection(Bool_t, AliAnalysisManager*) {} + /** + * Do not create a centrality selection + */ + void CreateCentralitySelection(Bool_t, AliAnalysisManager*) {} + /** + * Do not create an output handler + */ + AliVEventHandler* CreateOutputHandler(UShort_t) { return 0; } + /** + * The train class name + * + * @return Train class name + */ + const char* ClassName() const { return "ELossTimeTrain"; } + /** + * Overloaded to create new dNdeta.C and dndeta.sh in the output + * directory + * + * @param asShellScript + */ + void SaveSetup(Bool_t asShellScript) + { + TrainSetup::SaveSetup(asShellScript); + + SaveDraw(); + } + void SaveDraw() + { + std::ofstream o("Draw.C"); + o << "// Written by " << ClassName() << "\n" + << "void Draw(const char* fileName=\"AnalysisResults.root\")\n" + << "{\n" + << " gSystem->AddIncludePath(\"-DNO_TRAIN -DSUMMARY\");\n" + << " const char* fwd = \"$ALICE_ROOT/PWGLF/FORWARD/analysis2\";\n" + << " gSystem->AddIncludePath(Form(\"-I%s/scripts\", fwd));\n" + << " gROOT->SetMacroPath(Form(\"%s/trains:%s\", fwd,\n" + << " gROOT->GetMacroPath()));\n" + << " gROOT->LoadMacro(\"ELossTimeTrain.C+\");\n" + << " ELossTimeSummary s;\n" + << " s.Run(fileName);\n" + << "}\n" + << std::endl; + o.close(); + } + void PostShellCode(std::ostream& f) + { + f << " echo \"=== Draw results ...\"\n" + << " aliroot -l -b -q ${prefix}Draw.C\\(\\\"AnalysisResults.root\\\"\\)\n" + << std::endl; + } +}; +#endif +#ifdef SUMMARY +# include +# include + +struct ELossTimeSummary : public SummaryDrawer +{ + enum EFlags { + kEventInspector = 0x001, + }; + + void Run(const char* fname, UShort_t flags=0x01) + { + // --- Open the file ----------------------------------------------- + TString filename(fname); + TFile* file = TFile::Open(filename.Data(), "READ"); + if (!file) { + Error("Run", "Failed to open \"%s\"", filename.Data()); + return; + } + fPause = flags & kPause; + + // --- Make our canvas --------------------------------------------- + TString pdfName(filename); + pdfName.ReplaceAll(".root", ".pdf"); + CreateCanvas(pdfName, flags & kLandscape); + + // --- Make title page ------------------------------------------- + TCollection* c = GetCollection(file, "elossTimeSums"); + DrawTitlePage(c); + + if (flags & kEventInspector) DrawEventInspector(c); + + TH1* dt = GetH1(c, "dt"); + DrawInPad(fBody, 0, dt, "", kLogy); + PrintCanvas("#Deltat"); + + const char* rings[] = { "FMD1i", "FMD2i", "FMD2o", "FMD3o", "FMD3i", 0 }; + const char** pring = rings; + + while (*pring) { + DrawRing(c, *pring, dt); + pring++; + } + + CloseCanvas(); + } + void DrawTitlePage(TCollection* c) + { + fBody->cd(); + + Double_t y = .7; + TLatex* ltx = new TLatex(.5, y, "#Deltat vs #Delta/#Delta_{mip}"); + ltx->SetTextSize(0.07); + ltx->SetTextFont(62); + ltx->SetTextAlign(22); + ltx->SetNDC(); + ltx->Draw(); + + Double_t save = fParName->GetTextSize(); + fParName->SetTextSize(0.03); + fParVal->SetTextSize(0.03); + y = .6; + + TCollection* fc = GetCollection(c, "fmdEventInspector"); + UShort_t sys; + UShort_t sNN; + ULong_t runNo; + GetParameter(fc, "sys", sys); + GetParameter(fc, "sNN", sNN); + GetParameter(fc, "runNo", runNo); + + DrawParameter(y, "Run #", Form("%lu", runNo)); + TString tS; SysString(sys, tS); DrawParameter(y, "System", tS); + TString tE; SNNString(sNN, tE); DrawParameter(y, "#sqrt{s_{NN}}", tE); + + PrintCanvas("Title page"); + fParName->SetTextSize(save); + fParVal->SetTextSize(save); + } + void DrawRing(TCollection* c, const char* ring, TH1* dt) + { + TCollection* lring = GetCollection(c, ring); + if (!lring) return; + + TH2* h2 = GetH2(lring, "dtVsELoss"); + if (!h2) return; + + THStack* stack = new THStack(ring, ring); + // stack->SetTitle(ring); + + THStack* ratios = new THStack(Form("Ratios for %s",ring), ring); + // stack->SetTitle(ring); + + Printf(ring); + Int_t j = 2; + TH1* first = 0; + Double_t lfirst = 0; + Double_t rmax = 0; + Double_t max = 3; + for (Int_t i = 1; i <= h2->GetNbinsY(); i++) { + TH1* h = h2->ProjectionX(Form("%s_%03d", ring, i), i,i); + Double_t logDt = h2->GetYaxis()->GetBinCenter(i); + + Int_t nFill = h->GetEntries(); + if (nFill <= 1000) continue; + Double_t norm = dt->GetBinContent(i); + if (norm <= 1e-6) { + Warning("", "Normalization=%f<1e-6 but got " + "%d>1000 entries for log10(dt)=%5.3f", norm, nFill, logDt); + continue; + } + if (!first && logDt > TMath::Log10(25.)) { + lfirst = logDt; + first = h; + } + // Info("", "Normalization is %f", norm); + h->Sumw2(); + h->Scale(1. / norm); + h->SetTitle(Form("log_{10}(#Deltat)=%5.3f", logDt)); + + Float_t r, g, b; + TColor::HSV2RGB((j-1)*45, 1, .8, r, g, b); + Int_t col = TColor::GetColor(r, g, b); + j++; + h->SetLineColor(col); + h->SetLineStyle(j % 3+1); + h->SetLineWidth(2); + // h->SetFillColor(col); + // h->SetFillStyle(3002); + stack->Add(h); + + if (h == first) continue; + TH1* rh = static_cast(h->Clone(Form("ratio%s", h->GetName()))); + // rh->SetTitle(Form("log_{10}(#Deltat)=%5.3f", logDt)); + rh->Divide(first); + for (Int_t k = 1; k <= rh->GetNbinsX(); k++) { + if (rh->GetXaxis()->GetBinCenter(k) > max) break; + rmax = TMath::Max(rmax, rh->GetBinContent(k)); + } + + ratios->Add(rh); + } + Double_t savX = fParVal->GetX(); + Double_t savY = fParVal->GetY(); + fParVal->SetX(0.12); + fParVal->SetY(0.12); + fBody->Divide(1,2,0,0); + DrawInPad(fBody,1,stack,"nostack hist", kLogy|kLegend); + stack->GetXaxis()->SetRangeUser(-.1,max); + stack->GetYaxis()->SetTitle("1/N_{ev} dN/d(#Delta/#Delta_{mip})"); + + fParVal->SetX(0.6); + fParVal->SetY(0.4); + DrawInPad(fBody,2,ratios,"nostack hist", kLegend); + ratios->GetXaxis()->SetRangeUser(-.1, max); + ratios->GetXaxis()->SetTitle("#Delta/#Delta_{mip}"); + ratios->GetYaxis() + ->SetTitle(Form("X/(1/N_{ev}dN/d(#Delta/#Delta_{mip}))|_{%5.3f}",lfirst)); + Printf("Max: %f (%f)", ratios->GetMaximum(), rmax); + ratios->SetMaximum(rmax*1.2); + + + fParVal->SetX(savX); + fParVal->SetY(savY); + PrintCanvas(ring); + } +}; + +#endif +// EOF diff --git a/PWGLF/FORWARD/analysis2/trains/EventTimeTask.C b/PWGLF/FORWARD/analysis2/trains/EventTimeTask.C new file mode 100644 index 00000000000..95d49de063b --- /dev/null +++ b/PWGLF/FORWARD/analysis2/trains/EventTimeTask.C @@ -0,0 +1,776 @@ +// EventTimeTask.C +#ifndef __CINT__ +# include +# include +#else +class TTree; +// Force load of libGui +class TGWindow; +#endif +class AliESDEvent; + +/** + * Data structure to be filled by task - one for each event + * + */ +struct EventTimeData +{ + /** + * Widths of field in the full timestamp + */ + enum { + kBCWidth = 12, + kOrbitWidth = 24, + kPeriodWidth = 28 + }; + /** Full time stamp */ + ULong64_t fFull; + /** LDC creation time - standard time_t */ + UInt_t fTime; + /** Mask of detectors present in the event - FMD is 0x1000 */ + UInt_t fDetectors; + /** Type of event - 7 is physics */ + UShort_t fType; + /** + * Create a branch in a tree + * + * @param tree Tree to create the branch in + */ + void CreateBranch(TTree* tree) + { + tree->Branch("event", &(this->fFull), + "full/l:time/i:detector:type/s"); + } + /** + * Set the address of a branch for reading back objects from the tree + * + * @param tree Tree to read back from + */ + void ReadBranch(TTree* tree) + { + tree->SetBranchAddress("event", &(this->fFull)); + } + /** + * Utility function to encode the full time stamp from components. + * + * @param period Period counter (overflow of orbit counter) + * @param orbit Orbit counter (period of 88us) + * @param bc Bunch crossing number (period of 25ns) + * + * @return Encoded full time stamp + */ + static ULong64_t EncodeFull(ULong64_t period, + ULong64_t orbit, + ULong64_t bc) + { + const ULong64_t bcMask = (1 << kBCWidth) - 1; + const ULong64_t orbitMask = (1 << kOrbitWidth) - 1; + const ULong64_t periodMask = (1 << kPeriodWidth) - 1; + ULong64_t ret = ((bc & bcMask) | + (orbit & orbitMask) << kBCWidth | + (period & periodMask) << (kBCWidth+kOrbitWidth)); + return ret; + } + /** + * Fill information from ESD into this data structure. + * + * @param esd Event + * @param dets List of active detectors in this event. + */ + void Fill(AliESDEvent* esd, UInt_t dets); +}; + +#ifndef NO_TASK +# ifndef __CINT__ +# include +# endif +inline void EventTimeData::Fill(AliESDEvent* esd, UInt_t dets) +{ + ULong64_t period = esd->GetPeriodNumber(); + ULong64_t orbit = esd->GetOrbitNumber(); + ULong64_t bc = esd->GetBunchCrossNumber(); + fType = esd->GetEventType(); + fTime = esd->GetTimeStamp();//LDC time + fDetectors = dets; // esd->GetDAQDetectorPattern(); + fFull = EncodeFull(period, orbit, bc); +} +# else +inline void EventTimeData::Fill(AliESDEvent*, UInt_t) +{ + Warning("Fill", "Calling empty method - shouldn't happen"); +} +#endif + +#ifndef NO_TASK +# ifndef __CINT__ +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# else +class AliAnalysisManager; +class TTree; +class AliTimeStamp; +class TList; +class TH2; +# endif +# include +# include + +/** + * A task to record the unique timestamp of each event. + * + * @par Input: ESD + * @par Output: A tree with a single branch + */ +class EventTimeTask : public AliAnalysisTaskSE +{ +public: + enum { + kListSlot = 1, + kTreeSlot = 2 + }; + /** + * Default CTOR - for I/O only. + */ + EventTimeTask() + : AliAnalysisTaskSE(), + fTree(0), + fHistograms(0), + fDetVsType(0) + { + } + /** + * User constructor + * + * @param name Name of task + */ + EventTimeTask(const char* name) + : AliAnalysisTaskSE(name), + fTree(0), + fHistograms(0), + fDetVsType(0) + { + DefineOutput(kListSlot, TList::Class()); + DefineOutput(kTreeSlot, TTree::Class()); + fBranchNames = "ESD:AliESDRun.,AliESDHeader."; + } + /** + * Create user output objects. + * + * Called on each slave at start of job + */ + void UserCreateOutputObjects() + { + Printf("Creating tree and histogram"); + fHistograms = new TList(); + fHistograms->SetOwner(); + fHistograms->SetName("L"); + + fDetVsType = new TH2D("detVsType", "Detector vs type", + 16, -.5, 15.5, 31, -.5, 30.5); + fDetVsType->SetXTitle("Type"); + fDetVsType->SetYTitle("Detector"); + fDetVsType->SetDirectory(0); + fHistograms->Add(fDetVsType); + Printf("Histogram (%d,%f,%f)x(%d,%f,%f)", + fDetVsType->GetXaxis()->GetNbins(), + fDetVsType->GetXaxis()->GetXmin(), + fDetVsType->GetXaxis()->GetXmax(), + fDetVsType->GetYaxis()->GetNbins(), + fDetVsType->GetYaxis()->GetXmin(), + fDetVsType->GetYaxis()->GetXmax()); + + // TDirectory* savdir = gDirectory; + // Printf("Opening file at slot %d", kTreeSlot); + // OpenFile(kTreeSlot); + Printf("Make tree and disassociate from file"); + fTree = new TTree("T", "T"); + fTree->SetDirectory(0); + Printf("Create branch"); + fData.CreateBranch(fTree); + // savdir->cd(); + + + PostData(kListSlot, fHistograms); + PostData(kTreeSlot, fTree); + } + /** + * Analyse a single event + * + */ + void UserExec(Option_t*) + { + static Bool_t first = true; + LoadBranches(); + + AliESDEvent* esd = dynamic_cast(InputEvent()); + if (!esd) return; + if (!esd->GetHeader()) return; + + if (first) { + LoadTriggerConfig(esd->GetRunNumber()); + first = false; + } + ULong64_t mask = esd->GetTriggerMask(); + UInt_t dets = 0; + for (UShort_t i = 0; i < fDets.size(); i++) { + if ((1 << i) & mask) dets |= fDets[i]; + } + // Printf("Event mask 0x%016llx -> 0x%08x", mask, dets); + fData.Fill(esd, dets); + fTree->Fill(); + + UInt_t type = esd->GetEventType(); + UInt_t detectors = dets; + for (UInt_t i = 0; i < 31; i++) { + if ((1 << i) & detectors) fDetVsType->Fill(type, i); + } + + PostData(kListSlot, fHistograms); + PostData(kTreeSlot, fTree); + } + void LoadTriggerConfig(Int_t runNo) + { + Printf("Loading trigger configuration for run %d", runNo); + // --- Connect to CDB -------------------------------------------- + AliCDBManager* cdb = AliCDBManager::Instance(); + cdb->SetDefaultStorageFromRun(runNo); + cdb->SetRun(runNo); + + // --- Get entry ------------------------------------------------- + AliCDBEntry* entry = cdb->Get("GRP/CTP/Config"); + if (!entry || !entry->GetObject()) { + Warning("LoadTriggerConfig", "Couldn't get trigger configuration"); + return; + } + AliTriggerConfiguration* config = + static_cast(entry->GetObject()); + + // --- Get the classes, and resize cache ------------------------- + const TObjArray& clss = config->GetClasses(); + fDets.resize(clss.GetEntries()); + + // --- Loop over configurations ---------------------------------- + TIter next(&clss); + AliTriggerClass* cls = 0; + while ((cls = static_cast(next()))) { + Int_t mask = cls->GetMask(); + AliTriggerCluster* cluster = cls->GetCluster(); + if (!cluster) { + Warning("LoadTriggerConfig", + "Failed to get trigger cluster for %s", cls->GetName()); + continue; + } + TString names = cluster->GetDetectorsInCluster(); + if (names.IsNull()) { + Warning("LoadTriggerConfig", "No detectors for cluster %s", + cls->GetName()); + continue; + } + UInt_t dets = AliDAQ::DetectorPattern(names); + UInt_t id = UInt_t(TMath::Log2(mask)); + fDets[id] = dets; + Printf("Trigger mask 0x%08x (%3d): 0x%08x (%s)", + mask, id, dets, names.Data()); + } + // cdb->Destroy(); + } + /** + * Register with manager and connect output containers + * + * @param mgr Manager + */ + void Connect(AliAnalysisManager* mgr) + { + mgr->AddTask(this); + mgr->ConnectInput(this, 0, mgr->GetCommonInputContainer()); + AliAnalysisDataContainer* contTree = + mgr->CreateContainer("T", TTree::Class(), + AliAnalysisManager::kOutputContainer, + "time.root"); + AliAnalysisDataContainer* contList = + mgr->CreateContainer("L", TList::Class(), + AliAnalysisManager::kOutputContainer, + "hist.root"); + mgr->ConnectOutput(this, kTreeSlot, contTree); + mgr->ConnectOutput(this, kListSlot, contList); + } + /** + * Create an instance of this task, and register it and it's + * outputs. + */ + static void Create() + { + EventTimeTask* task = new EventTimeTask("eventTime"); + task->Connect(AliAnalysisManager::GetAnalysisManager()); + } + TTree* fTree; // Our tree + EventTimeData fData; // Our data + TList* fHistograms; // List + TH2D* fDetVsType; // Histogram + std::vector fDets; // Per-trigger-bit detector mask + + ClassDef(EventTimeTask,3); +}; +#endif // NO_TASK + +#ifndef NO_MAP +#include +#include + +typedef std::pair EventTimeMapPair; + +/** + * A map of event time-stamp to distance to previous event + * + */ +struct EventTimeMap : public TObject +{ + /** Map type */ + typedef std::map Map; + /** Key type */ + typedef Map::key_type Key; + /** Mapped value type */ + typedef Map::mapped_type Value; + /** Element type */ + typedef Map::value_type Pair; + /** Iterator type */ + typedef Map::iterator Iterator; + /** Constant iterator type */ + typedef Map::const_iterator ConstIterator; + /** + * Constructor + */ + EventTimeMap() : TObject(), fMap() {} + /** + * Destructor + */ + virtual ~EventTimeMap() {} + /** + * Get name of this object - always the same + * + * @return The string "eventTimeMap" + */ + const char* GetName() const { return "eventTimeMap"; } + /** + * Element access. If the key @a k doesn't already exist, it is + * created + * + * @param k Key + * + * @return A pair of key and value + */ + Value& operator[](const Key& k) + { + return fMap[k]; + } + /** + * Find the element whos key is @a k + * + * @param k Key to look for + * + * @return Iterator pointing to element, or end of container if not found + */ + Iterator Find(const Key& k) + { + return fMap.find(k); + } + /** + * Find the element whos key is @a k + * + * @param k Key to look for + * + * @return Iterator pointing to element, or end of container if not found + */ + ConstIterator Find(const Key& k) const + { + return fMap.find(k); + } + /** + * Get forward iterator pointing beginning of the container + * + * @return Iterator to start of container + */ + Iterator Begin() + { + return fMap.begin(); + } + /** + * Get forward iterator pointing beginning of the container + * + * @return Iterator to start of container + */ + ConstIterator Begin() const + { + return fMap.begin(); + } + /** + * Get forward iterator pointing just beyond the end of the container + * + * @return Iterator just beyond container + */ + Iterator End() + { + return fMap.end(); + } + /** + * Get forward iterator pointing just beyond the end of the container + * + * @return Iterator just beyond container + */ + ConstIterator End() const + { + return fMap.end(); + } + enum { + kInvalidTime = 0xFFFFFFFFFFFFFFFF + }; + /** + * Get the time difference to previous event from a event with a + * given time stamp. + * + * @param timestamp Time stamp of the event + * + * @return time difference or kInvalidTime if not found + */ + Value Get(const Key& timestamp) const + { + ConstIterator i = Find(timestamp); + if (i == End()) return kInvalidTime; + return i->second; + } + /** + * Set the time to previous event for a given event time stamp + * + * @param timestamp Event time stamp + * @param diff Time to previous event + */ + void Set(const Key& timestamp, const Value& diff) + { + this->operator[](timestamp) = diff; + } + ULong64_t Size() const { return fMap.size(); } + /** Our embeded map */ + Map fMap; + ClassDef(EventTimeMap,1) +}; +#endif // NO_MAP + +#ifndef NO_SORTER +# ifndef __CINT__ +# include +# include +# include +# include +# include +# include +# else +class TFile; +class TTree; +class TCanvas; +# endif + +/** + * A class to sort the tree and generate our timestamp->dT map. + */ +struct EventTimeSorter +{ + TTree* fTree; + EventTimeData fData; + ULong64_t fMin; + ULong64_t fMax; + + /** + * Constructor + */ + EventTimeSorter() : fTree(0), fData(), fMin(0xFFFFFFFFFFFFFFFF), fMax(0) {} + /** + * Progress meter + * + * @param cur Current entry + * @param total Total number of entries + */ + void Progress(Long64_t cur, Long64_t total) const + { + static UShort_t old = 0; + if (cur == 0) old = 0; + UShort_t now = UShort_t(100 * (cur + 1 == total ? 1 : + Double_t(cur)/total)); + if (now > old) + std::cout << "\rLooping over " << total << " entries: ... " + << now << '%' << std::flush; + if (now >= 100) std::cout << std::endl; + old = now; + } + /** + * Connect a tree + * + * @param inputName + * @param treeName + * + * @return + */ + Bool_t OpenInput(const char* inputName, const char* treeName) + { + CloseInput(); + + // --- Get input ------------------------------------------------- + TFile* input = TFile::Open(inputName, "READ"); + if (!input) { + Error("Run", "Failed to open \"%s\"", inputName); + return false; + } + + fTree = static_cast(input->Get(treeName)); + if (!fTree) { + Error("Run", "Couldn't get tree \"%s\" from \"%s\"", treeName,inputName); + return false; + } + + // --- Set branch address --------------------------------------- + fData.ReadBranch(fTree); + + return true; + } + /** + * Disconnect tree + */ + void CloseInput() + { + if (!fTree) return; + TFile* file = fTree->GetCurrentFile(); + if (file) file->Close(); + fTree = 0; + } + /** + * Run the sorter + * + * @param inputName Input file name + * @param outputName Output file name + * @param treeName Tree name + * + * @return true on success + */ + Bool_t Run(const char* inputName, const char* outputName, + const char* treeName="T") + { + if (!OpenInput(inputName, treeName)) return false; + + // --- Loop over the data ---------------------------------------- + Long64_t nEnt = fTree->GetEntries(); + Long64_t n = 0; + ULong64_t* values = new ULong64_t[nEnt]; + UInt_t* dets = new UInt_t[nEnt]; + UShort_t* types = new UShort_t[nEnt]; + UInt_t* times = new UInt_t[nEnt]; + for (Long64_t iEnt = 0; iEnt < nEnt; iEnt++) { + Progress(iEnt, nEnt); + + fTree->GetEntry(iEnt); + + if (!(fData.fDetectors & 0x1000)) + // No FMD read-out + continue; + if (fData.fType != 7) { + // Ignore all but physics events + Warning("Run", "Non-PHYSICS event (%d) with FMD in it", fData.fType); + // continue; + } + + // Store values + Int_t j = n; + values[j] = fData.fFull; + dets[j] = fData.fDetectors; + types[j] = fData.fType; + times[j] = fData.fTime; + n++; + } + + // --- Now sort all values --------------------------------------- + TArrayL64 index(n); + TMath::Sort(n, values, index.fArray, false); + + // Maps the unique time to the distance to previous event + EventTimeMap eventTimeMap; + ULong64_t key = values[index[0]]; + ULong64_t prev = 0; + ULong64_t dt = key-prev; + ULong64_t min = 0xFFFFFFFFFFFFFFFF; + ULong64_t max = 0; + eventTimeMap[key] = dt; + for (Long64_t i = 1; i < n; i++) { + Long64_t j = index[i]; + Long64_t k = index[i-1]; + key = values[j]; + prev = values[k]; + dt = key - prev; + if (dt <= 0) { + Printf("0x%016llx==0x%016llx -> dt=%10llu [%10lld %10lld] " + "(det: 0x%08x 0x%08x type: 0x%04x 0x%04x time: %08d %08d)", + key, prev, dt, j, k, dets[j], dets[k], + types[j], types[k], times[j], times[k]); + // continue; + } + eventTimeMap[key] = dt; + min = TMath::Min(dt, min); + max = TMath::Max(dt, max); + } + std::cout << "Range is [" << min << "," << max << ']' << std::endl; + + TFile* output = TFile::Open(outputName, "RECREATE"); + if (!output) { + Error("Run", "Failed to create output file \"%s\"", outputName); + return false; + } + fMin = min; + fMax = max; + eventTimeMap.Write(); + output->Write(); + output->Close(); + + delete [] values; + + CloseInput(); + + return true; + } + Bool_t Test(const char* inputName, const char* outputName, + const char* treeName="T") + { + if (!OpenInput(inputName, treeName)) return false; + + // --- Get our map ----------------------------------------------- + TFile* output = TFile::Open(outputName, "UPDATE"); + if (!output) { + Error("Test", "Failed to open \"%s\"", outputName); + return false; + } + + EventTimeMap* eventTimeMap = + static_cast(output->Get("eventTimeMap")); + if (!eventTimeMap) { + Error("Test", "Failed to get \"%s\" from \"%s\"", + "eventTimeMap", outputName); + return false; + } + + // --- Histograms -------------------------------------------------- + ULong64_t mmin = TMath::Min(25*fMin, 900000ULL); + ULong64_t mmax = TMath::Min(25*fMax, 110000000ULL); + TH1* diff = new TH1D("diff", "Time-to-last-event (10#mus bins w/FMD)", + (mmax-mmin)/10000, mmin, mmax); + diff->SetStats(0); + diff->SetXTitle("#Deltat [ns]"); + diff->SetFillColor(kRed+1); + diff->SetFillStyle(3001); + diff->SetDirectory(0); + + + TH1* ldiff = new TH1D("ldiff", "log(#Deltat) - Events w/FMD", + 300, 0, 15); + ldiff->SetStats(0); + ldiff->SetXTitle("log_{10}(#Deltat) [ns]"); + ldiff->SetFillColor(kGreen+1); + ldiff->SetFillStyle(3001); + ldiff->SetDirectory(0); + + // --- Loop over the data ---------------------------------------- + Long64_t nEnt = fTree->GetEntries(); + Long64_t nZero = 0; + Long64_t nMiss = 0; + Long64_t nGood = 0; + Long64_t nNoFMD = 0; + for (Long64_t iEnt = 0; iEnt < /*10*/ nEnt; iEnt++) { + Progress(iEnt, nEnt); + fTree->GetEntry(iEnt); + + if (!(fData.fDetectors & 0x1000)) { + // No FMD read-out + nNoFMD++; + continue; + } + if (fData.fType != 7) { + // Ignore all but physics events + Warning("Run", "Non-PHYSICS event (%d) with FMD in it", fData.fType); + // continue; + } + + // Look-up the timestamp + ULong64_t value = fData.fFull; + ULong64_t dT = eventTimeMap->Get(value); + if (dT == EventTimeMap::kInvalidTime) { + Warning("Test", "Value %llu not found", value); + ldiff->Fill(1); + nMiss++; + continue; + } + if (dT <= 0) { +#if 0 + Warning("Test", "Impossible dt=%llu found for %llu (0x%0x %2d)", + dT, value, fData.fDetectors, fData.fType); +#endif + ldiff->Fill(1); + nZero++; + continue; + } + nGood++; + Double_t dt = 25.*dT; + diff->Fill(dt); + Double_t logDt = TMath::Log10(dt); + ldiff->Fill(logDt); + } + CloseInput(); + Printf("missing: %llu, bad: %llu, good: %llu, no FMD: %llu, all: %llu", + nMiss,nZero,nGood,nNoFMD,nEnt); + + // --- Draw results ---------------------------------------------- + TCanvas* c = new TCanvas("c", "C"); + c->SetTopMargin(0.01); + c->SetRightMargin(0.01); + c->Divide(2,1); // ,0,0); + TVirtualPad* p = c->cd(2); + // p->SetRightMargin(0.10); + p->SetLogy(); + ldiff->DrawCopy(); + + p = c->cd(1); + p->SetLogy(); + p->SetLogx(); + diff->DrawCopy(); + + c->Print("dt.png"); + c->Print("dt.C"); + + // --- Write histogram to file ----------------------------------- + output->cd(); + diff->Write(); + output->Write(); + output->Close(); + + return true; + } +}; +#endif // NO_SORTER +#ifdef __MAKECINT__ +# ifndef NO_MAP +# pragma link C++ class std::pair+; +# pragma link C++ class std::map+; +# endif +# ifndef NO_TASK +# pragma link C++ class std::vector+; +# endif +#endif + +// EOF + + diff --git a/PWGLF/FORWARD/analysis2/trains/EventTimeTrain.C b/PWGLF/FORWARD/analysis2/trains/EventTimeTrain.C new file mode 100644 index 00000000000..1c9f75ae60d --- /dev/null +++ b/PWGLF/FORWARD/analysis2/trains/EventTimeTrain.C @@ -0,0 +1,104 @@ +// EventTimeTrain.C +#ifndef __CINT__ +#include +#include +#else +class AliAnalysisManager; +#endif +#include "TrainSetup.C" +#include "ParUtilities.C" + +/** + * Train to record time of each event + * + */ +class EventTimeTrain : public TrainSetup +{ +public: + /** + * Constructor + * + * @param name The name + */ + EventTimeTrain(const char* name="eventTime") : TrainSetup(name) + { + fOptions.Set("type", "ESD"); + } + /** + * Create our tasks + * + * @param mgr Manager + */ + void CreateTasks(AliAnalysisManager*) + { + if (!ParUtilities::MakeScriptPAR(fHelper->Mode() == Helper::kLocal, + "EventTimeTask.C", + // Gui because of CDB + // XMLParser because of CDB + // CDB because look-up of trigger config + "Gui,XMLParser," + "STEERBase,CDB,ESD,AOD,ANALYSIS,OADB," + "ANALYSISalice", + fHelper)) + Fatal("","Failed to make PAR"); + fHelper->LoadLibrary("EventTimeTask"); + gROOT->ProcessLine("EventTimeTask::Create()"); + } + /** + * Do not create a physics selection + */ + void CreatePhysicsSelection(Bool_t, AliAnalysisManager*) {} + /** + * Do not create a centrality selection + */ + void CreateCentralitySelection(Bool_t, AliAnalysisManager*) {} + /** + * Do not create an output handler + */ + AliVEventHandler* CreateOutputHandler(UShort_t) { return 0; } + /** + * The train class name + * + * @return Train class name + */ + const char* ClassName() const { return "EventTimeTrain"; } + /** + * Overloaded to create new dNdeta.C and dndeta.sh in the output + * directory + * + * @param asShellScript + */ + void SaveSetup(Bool_t asShellScript) + { + TrainSetup::SaveSetup(asShellScript); + + SaveSort(); + } + void SaveSort() + { + std::ofstream o("Sort.C"); + o << "// Written by " << ClassName() << "\n" + << "void Sort(const char* prefix=\"\",\n" + << " const char* fileName=\"time.root\",\n" + << " const char* outName=\"map.root\",\n" + << " const char* treeName=\"T\")\n" + << "{\n" + << " gSystem->AddIncludePath(\"-DNO_TASK -I$ALICE_ROOT/include\");\n" + << " TString mac(\"EventTimeTask/EventTimeTask.C+g\");\n" + << " if (prefix && prefix[0] != '\\0') mac.Prepend(prefix);\n" + << " gROOT->LoadMacro(mac);\n" + << " EventTimeSorter s;\n" + << " if (!s.Run(fileName,outName,treeName)) return;\n" + << " s.Test(fileName,outName,treeName);\n" + << "}\n" + << std::endl; + o.close(); + } + void PostShellCode(std::ostream& f) + { + f << " echo \"=== Sort results ...\"\n" + << " aliroot -l -b -q ${prefix}Sort.C\\(\\\"${prefix}\\\"\\)\n" + << std::endl; + } +}; +// EOF diff --git a/PWGLF/FORWARD/analysis2/trains/MakeAODTrain.C b/PWGLF/FORWARD/analysis2/trains/MakeAODTrain.C index 3f97b380e4f..56ee98b7c3c 100644 --- a/PWGLF/FORWARD/analysis2/trains/MakeAODTrain.C +++ b/PWGLF/FORWARD/analysis2/trains/MakeAODTrain.C @@ -264,7 +264,7 @@ protected: << "// 0x200 Pause\n" << "//\n" << "void Summarize(const char* filename=\"forward.root\",\n" - << " UShort_t what=0x0FF)\n" + << " UShort_t what=0x1FF)\n" << "{\n" << " const char* fwd=\"$ALICE_ROOT/PWGLF/FORWARD/analysis2\";\n" << " gROOT->LoadMacro(Form(\"%s/DrawAODSummary.C\",fwd));\n" diff --git a/PWGLF/FORWARD/analysis2/trains/MakeMultDistsTrain.C b/PWGLF/FORWARD/analysis2/trains/MakeMultDistsTrain.C index f8d702b427e..60532c56948 100644 --- a/PWGLF/FORWARD/analysis2/trains/MakeMultDistsTrain.C +++ b/PWGLF/FORWARD/analysis2/trains/MakeMultDistsTrain.C @@ -33,7 +33,7 @@ public: fOptions.Add("trig", "TYPE", "Trigger type", "V0AND"); fOptions.Add("vzMin", "CENTIMETER", "Min Ip Z", -4); fOptions.Add("vzMax", "CENTIMETER", "Max Ip Z", +4); - fOptions.Add("phi-acc", "Use stored phi acceptance", true); + fOptions.Add("phi-acc", "Use stored phi acceptance", false); fOptions.Add("asymmetric", "Make asymmetric (+/-) bins", false); } protected: