// ROOT #include "TFile.h" #include "TList.h" #include "TH1.h" #include "TH2.h" #include "TH3.h" #include "TRandom.h" // analysis framework #include "AliAnalysisManager.h" #include "AliInputEventHandler.h" #include "AliVEvent.h" #include "AliVTrdTrack.h" // MC stuff #include "AliMCEvent.h" #include "AliGenPythiaEventHeader.h" // ESD stuff #include "AliESDEvent.h" #include "AliESDInputHandler.h" #include "AliESDtrack.h" #include "AliESDTrdTrack.h" #include "AliESDTrdTracklet.h" #include "AliESDTrdTrigger.h" // AOD stuff #include "AliAODEvent.h" #include "AliAODJet.h" #include "AliAODTrack.h" // jet tasks #include "AliAnalysisTaskJetServices.h" #include "AliAnalysisHelperJetTasks.h" #include "AliAnalysisTaskJetsTriggerTRD.h" AliAnalysisTaskJetsTriggerTRD::AliAnalysisTaskJetsTriggerTRD(const char *name) : AliAnalysisTaskSE(name), fTriggerMask(0), fTrdTrg(), fOutputList(), fHist(), fShortTaskId("jets_trg_trd"), fNoTriggers(kTrgLast), fNoJetPtBins(80), fJetPtBinMax(400), fAvgXsection(0.), fAvgTrials(0.), fPtHard(0.), fNTrials(0), fGlobalEfficiencyGTU(.8), fGtuLabel(-1) // -3 for hw, -1821 for re-simulation { // default ctor DefineOutput(1, TList::Class()); } AliAnalysisTaskJetsTriggerTRD::~AliAnalysisTaskJetsTriggerTRD() { // dtor } void AliAnalysisTaskJetsTriggerTRD::UserCreateOutputObjects() { // create user output objects // setup list OpenFile(1); fOutputList = new TList(); fOutputList->SetOwner(); // setup histograms TH1 *histStat = AddHistogram(ID(kHistStat), "event statistics;;counts", kStatLast-1, .5, kStatLast-.5); histStat->GetXaxis()->SetBinLabel(ID(kStatSeen)); histStat->GetXaxis()->SetBinLabel(ID(kStatTrg)); histStat->GetXaxis()->SetBinLabel(ID(kStatUsed)); histStat->GetXaxis()->SetBinLabel(ID(kStatEvCuts)); if (HasMC()) { fNoTriggers = kTrgMCLast; AddHistogram(ID(kHistXsection), "xsection stats", 1, 0., 1.); AddHistogram(ID(kHistPtHard), "pt hard;#hat{p}_{T} (GeV/c);counts", fNoJetPtBins, 0., fJetPtBinMax); AddHistogram(ID(kHistJetPtMC), "leading jet spectrum (MC, |#eta| < 0.5);p_{T}^{jet} (GeV/c);counts", fNoJetPtBins, 0., fJetPtBinMax); } AddHistogram(ID(kHistJetEtaAvg), "average eta", 100, -1., 1.); AddHistogram(ID(kHistNoJets), "number of jets;N^{jet};counts", 400, -.5, 399.5); AddHistogram(ID(kHistTrackGTU), "GTU track p_{T};p_{T} (GeV/c);counts", 100, 0., 25.); AddHistogram(ID(kHistTrackEffGTU), "found GTU tracks;p_{T} (GeV/c);#eta;#varphi", 100, 0., 25., 100, -1., 1., 100, 0., 2.*TMath::Pi()); AddHistogram(ID(kHistTrackEffMC), "wanted GTU tracks;p_{T} (GeV/c);#eta;#varphi", 100, 0., 25., 100, -1., 1., 100, 0., 2.*TMath::Pi()); AddHistogram(ID(kHistNPtMin), "rejection;p_{T}^{min};N_{trk};trigger", 100, 0., 10., 20, 0., 20., fNoTriggers - 1, .5, fNoTriggers - .5); // leading jet AddHistogram(ID(kHistLeadJetPt), "leading jet spectrum (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);counts", fNoJetPtBins, 0., fJetPtBinMax, fNoTriggers - 1, .5, fNoTriggers - .5); AddHistogram(ID(kHistLeadJetPtEta), "leading jet #eta (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger;#eta", fNoJetPtBins, 0., fJetPtBinMax, fNoTriggers - 1, .5, fNoTriggers - .5, 40, -.5, .5); AddHistogram(ID(kHistLeadJetPtPhi), "leading jet #phi (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger;#varphi", fNoJetPtBins, 0., fJetPtBinMax, fNoTriggers - 1, .5, fNoTriggers - .5, 40, 0., 2. * TMath::Pi()); AddHistogram(ID(kHistLeadJetEtaPhi), "leading jet #eta - #varphi (|#eta| < 0.5);#eta;trigger;#varphi", 40, -.5, .5, fNoTriggers - 1, .5, fNoTriggers - .5, 40, 0., 2. * TMath::Pi()); AddHistogram(ID(kHistLeadJetPtTrackPt), "leading jet pt vs track pt;p_{T} (GeV/c);trigger;p_{T}^{jet,ch} (GeV/c)", 100., 0., 20., fNoTriggers - 1, .5, fNoTriggers - .5, fNoJetPtBins, 0., fJetPtBinMax); AddHistogram(ID(kHistLeadJetPtZ), "leading jet pt vs z;z;trigger;p_{T}^{jet,ch} (GeV/c)", 100., 0., 1., fNoTriggers - 1, .5, fNoTriggers - .5, fNoJetPtBins, 0., fJetPtBinMax); AddHistogram(ID(kHistLeadJetPtXi), "leading jet pt vs #xi;#xi;trigger;p_{T}^{jet,ch} (GeV/c)", 100., 0., 10., fNoTriggers - 1, .5, fNoTriggers - .5, fNoJetPtBins, 0., fJetPtBinMax); // inclusive jets AddHistogram(ID(kHistJetPt), "jet spectrum (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger", fNoJetPtBins, 0., fJetPtBinMax, fNoTriggers - 1, .5, fNoTriggers - .5); AddHistogram(ID(kHistJetPtEta), "jet #eta (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger;#eta", fNoJetPtBins, 0., fJetPtBinMax, fNoTriggers - 1, .5, fNoTriggers - .5, 40, -.5, .5); AddHistogram(ID(kHistJetPtPhi), "jet #phi (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger;#varphi", fNoJetPtBins, 0., fJetPtBinMax, fNoTriggers - 1, .5, fNoTriggers - .5, 40, 0., 2. * TMath::Pi()); AddHistogram(ID(kHistJetEtaPhi), "jet #eta - #varphi (|#eta| < 0.5);#eta;trigger;#varphi", 40, -.5, .5, fNoTriggers - 1, .5, fNoTriggers - .5, 40, 0., 2. * TMath::Pi()); AddHistogram(ID(kHistJetPtITS), "jet spectrum (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger", fNoJetPtBins, 0., fJetPtBinMax, fNoTriggers - 1, .5, fNoTriggers - .5); AddHistogram(ID(kHistJetPt3x3), "jet spectrum (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger", fNoJetPtBins, 0., fJetPtBinMax, fNoTriggers - 1, .5, fNoTriggers - .5); AddHistogram(ID(kHistJetPtTrackPt), "jet pt vs track pt;p_{T} (GeV/c);trigger;p_{T}^{jet,ch} (GeV/c)", 100., 0., 20., fNoTriggers - 1, .5, fNoTriggers - .5, fNoJetPtBins, 0., fJetPtBinMax); AddHistogram(ID(kHistJetPtZ), "jet pt vs z;z;trigger;p_{T}^{jet,ch} (GeV/c)", 100., 0., 1., fNoTriggers - 1, .5, fNoTriggers - .5, fNoJetPtBins, 0., fJetPtBinMax); AddHistogram(ID(kHistJetPtXi), "jet pt vs #xi;#xi;trigger;p_{T}^{jet,ch} (GeV/c)", 100., 0., 10., fNoTriggers - 1, .5, fNoTriggers - .5, fNoJetPtBins, 0., fJetPtBinMax); for (Int_t iHist = kHistLeadJetPt; iHist <= kHistJetPtXi; ++iHist) { TH1 *h = GetHistogram(Hist_t (iHist)); h->GetYaxis()->SetBinLabel(ID(kTrgMinBias)); h->GetYaxis()->SetBinLabel(ID(kTrgInt7)); h->GetYaxis()->SetBinLabel(ID(kTrgInt8)); h->GetYaxis()->SetBinLabel(ID(kTrgEMC7)); h->GetYaxis()->SetBinLabel(ID(kTrgEMC8)); h->GetYaxis()->SetBinLabel(ID(kTrgPbPb)); h->GetYaxis()->SetBinLabel(ID(kTrgCentral)); h->GetYaxis()->SetBinLabel(ID(kTrgSemiCentral)); h->GetYaxis()->SetBinLabel(ID(kTrgInt7WUHJT)); h->GetYaxis()->SetBinLabel(ID(kTrgInt8WUHJT)); h->GetYaxis()->SetBinLabel(ID(kTrgEMC7WUHJT)); h->GetYaxis()->SetBinLabel(ID(kTrgEMC8WUHJT)); h->GetYaxis()->SetBinLabel(ID(kTrgEMCEJE)); h->GetYaxis()->SetBinLabel(ID(kTrgEMCEGA)); h->GetYaxis()->SetBinLabel(ID(kTrgInt7_WU)); h->GetYaxis()->SetBinLabel(ID(kTrgInt7_WUHJT)); h->GetYaxis()->SetBinLabel(ID(kTrgEMCEJE_WU)); h->GetYaxis()->SetBinLabel(ID(kTrgEMCEJE_WUHJT)); if (HasMC()) { h->GetYaxis()->SetBinLabel(ID(kTrgMC3x3Vtx)); h->GetYaxis()->SetBinLabel(ID(kTrgMC3x3TRD)); h->GetYaxis()->SetBinLabel(ID(kTrgMC3x3TRDeff)); h->GetYaxis()->SetBinLabel(ID(kTrgMC3x3TRDeffmap)); } } AddHistogram(ID(kHistJetPtNoTracks3), "number of tracks above 3 GeV;p_{T}^{jet,ch};no. of tracks", fNoJetPtBins, 0., fJetPtBinMax, 40, -.5, 39.5); PostData(1, fOutputList); } Bool_t AliAnalysisTaskJetsTriggerTRD::Notify() { // actions to be taken upon notification about input file change // ??? check ??? // we should only count the cross section we see in the analysis, // i.e. fXsection / nTrials if (HasMC()) { fAvgXsection = 0.; Float_t xSection = 0.; Float_t nTrials = 0.; Float_t nEntries = 0.; TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree(); if (tree) { TFile *curfile = tree->GetCurrentFile(); if (!curfile) { AliError("No current file"); return kFALSE; } if (!AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(), xSection, nTrials)) { AliError("retrieval of cross section failed"); } nEntries = (Float_t) tree->GetTree()->GetEntries(); if (nEntries > 0.) { fAvgTrials = nTrials / nEntries; if (nTrials > 0.) fAvgXsection = xSection / (nTrials / nEntries); } } printf("n_trials = %f\n", nTrials); } return AliAnalysisTaskSE::Notify(); } void AliAnalysisTaskJetsTriggerTRD::UserExec(Option_t * /* option */) { // actual work // setup pointers to input data (null if unavailable) // mcEvent: MC input // esdEvent: ESD input // outEvent: AOD output // aodEvent: AOD input if available, otherwise AOD output AliMCEvent *mcEvent = this->MCEvent(); AliESDEvent *esdEvent = dynamic_cast(this->InputEvent()); // could also be AOD input AliAODEvent* outEvent = this->AODEvent(); AliAODEvent *aodEvent = outEvent; if (dynamic_cast(this->InputEvent())) aodEvent = (AliAODEvent*) (this->InputEvent()); if ((fDebug > 0) && esdEvent) printf("event: %s-%06i\n", CurrentFileName(), esdEvent->GetEventNumberInFile()); Int_t nTracksMC = mcEvent ? mcEvent->GetNumberOfTracks() : 0; // no. of MC tracks // Int_t nTracks = InputEvent()->GetNumberOfTracks(); // no. of global tracks Int_t nTracksTRD = InputEvent()->GetNumberOfTrdTracks(); // no. of GTU tracks TList partMC; TList partEsd; TList partGtu; Float_t leadingJetPtMC = 0.; // leading jet energy from MC information Float_t leadingJetPtRec = 0.; // leading jet energy from AOD information // record number of sampled events and detect trigger contributions FillH1(kHistStat, kStatSeen); if (!DetectTriggers()) { AliError("Failed to detect the triggers"); return; } // only continue for events from interesting triggers if (fTriggerMask == 0) return; FillH1(kHistStat, kStatTrg); // no further technical requirements for the event at the moment FillH1(kHistStat, kStatUsed); if (HasMC()) { if (fAvgXsection == 0.) { AliError("zero cross section"); return; } FillH1(kHistXsection, .5, fAvgXsection); } // apply event cuts const AliVVertex *vtx = InputEvent()->GetPrimaryVertex(); if (!vtx || (vtx->GetNContributors() < 3.) || (vtx->GetZ() > 10.)) return; FillH1(kHistStat, kStatEvCuts); // extract MC information if (mcEvent) { // check for PYTHIA event header AliGenPythiaEventHeader *pythiaHeader = dynamic_cast (mcEvent->GenEventHeader()); if (!pythiaHeader) { AliWarning("MC event without PYTHIA event header!\n"); } else { fPtHard = pythiaHeader->GetPtHard(); Int_t nTrials = pythiaHeader->Trials(); fNTrials += nTrials; FillH1(kHistPtHard, fPtHard); // loop over jets from PYTHIA Float_t eta1 = -10.; Float_t eta2 = -10.; for (Int_t iJet = 0; iJet < pythiaHeader->NTriggerJets(); iJet++) { Float_t p[4]; pythiaHeader->TriggerJet(iJet, p); TLorentzVector pJet(p); Float_t pt = pJet.Pt(); // only consider jets with |eta| < 0.5 Float_t eta = pJet.Eta(); if (TMath::Abs(eta) > 0.5) continue; if (pt > leadingJetPtMC) leadingJetPtMC = pt; // for eta averge determination consider only jets above 60 GeV if (pt < 60.) continue; if (eta1 < -1.) eta1 = eta; else if (eta2 < -1.) eta2 = eta; } // fill histogram for leading jet pt spectrum FillH1(kHistJetPtMC, leadingJetPtMC); // fill histogram for eta average if ((eta1 > -1.) && (eta2 > -1.)) FillH1(kHistJetEtaAvg, (eta1 + eta2)/2.); } for (Int_t iTrack = 0; iTrack < nTracksMC; ++iTrack) { AliVParticle *part = mcEvent->GetTrack(iTrack); if (AcceptTrackMC(iTrack)) { FillH3(kHistTrackEffMC, part->Pt(), part->Eta(), part->Phi()); } } } // loop over GTU tracks for (Int_t iTrack = 0; iTrack < nTracksTRD; ++iTrack) { AliVTrdTrack *trk = InputEvent()->GetTrdTrack(iTrack); // printf("trk %p has pt %5.2f and label %i\n", // trk, trk->Pt(), trk->GetLabel()); if ((fGtuLabel != -1) && (trk->GetLabel() != fGtuLabel)) continue; FillH1(kHistTrackGTU, TMath::Abs(trk->Pt())); partGtu.Add(trk); } partGtu.Sort(kSortAscending); Int_t nTracksPerStack[90] = { 0 }; Int_t nTracksPerStackMax = 0; TIter nextPartGtu(&partGtu); while (AliVTrdTrack *trdTrack = (AliVTrdTrack*) nextPartGtu()) { // count no. of tracks in stack, // check whether this number was reached before, // if not store pt^min(n), // i.e. pt of current track because of sorting Int_t sec = trdTrack->GetSector(); Int_t stack = trdTrack->GetStack(); if ((sec > -1) && (sec < 18) && (stack > -1) && (stack < 5)) { ++nTracksPerStack[5*sec + stack]; if (nTracksPerStack[5*sec + stack] > nTracksPerStackMax) { ++nTracksPerStackMax; for (Int_t iTrigger = 1; iTrigger < fNoTriggers; ++iTrigger) if (IsTrigger(Trigger_t (iTrigger))) FillH3(kHistNPtMin, TMath::Abs(trdTrack->Pt()), nTracksPerStackMax, iTrigger); } if (HasMC()) { Int_t label = trdTrack->GetLabel(); if (label > -1) { if (AcceptTrackMC(label)) { AliVParticle *part = MCEvent()->GetTrack(label); FillH3(kHistTrackEffGTU, part->Pt(), part->Eta(), part->Phi()); } } else { AliWarning(Form("GTU track at %p with no label", trdTrack)); const Int_t nLayers = 6; for (Int_t iLayer = 0; iLayer < nLayers; ++iLayer) { AliVTrdTracklet *trkl = trdTrack->GetTracklet(iLayer); if (trkl) AliWarning(Form("tracklet in layer %i has label %i\n", iLayer, trkl->GetLabel())); } } } } else { AliError(Form("Invalid sector or stack: %i %i", sec, stack)); } } // loop over jets from AOD event if (aodEvent) { TClonesArray *jetArray = dynamic_cast (aodEvent->FindListObject(fJetBranchName)); if (jetArray) { Int_t nJets = jetArray->GetEntriesFast(); FillH1(kHistNoJets, nJets); AliAODJet *leadJet = 0x0; // AliAODJet *subleadJet = 0x0; for (Int_t iJet = 0; iJet < nJets; ++iJet) { AliAODJet *jet = (AliAODJet*) (*jetArray)[iJet]; if (TMath::Abs(jet->Eta()) < 0.5) { // check contributing tracks Int_t nJetTracks = jet->GetRefTracks()->GetEntriesFast(); Int_t nJetTracks3 = 0; AliAODTrack *leadingTrack = 0x0; for (Int_t iTrack=0; iTrack < nJetTracks; ++iTrack) { AliAODTrack *track = (AliAODTrack*) jet->GetRefTracks()->At(iTrack); // count constituents above 3 GeV/c if (track->Pt() > 3.) ++nJetTracks3; // find the leading track if (!leadingTrack || (track->Pt() > leadingTrack->Pt())) leadingTrack = track; } // find leading jet if (TMath::Abs(jet->Pt()) > leadingJetPtRec) { leadingJetPtRec = TMath::Abs(jet->Pt()); // subleadJet = leadJet; leadJet = jet; } // jet pt spectrum and // fragmentation function for (Int_t iTrigger = 1; iTrigger < fNoTriggers; ++iTrigger) if (IsTrigger(Trigger_t (iTrigger))) { Float_t jetPt = jet->Pt(); FillH1(kHistJetPt, jetPt, iTrigger); FillH3(kHistJetPtEta, jet->Pt(), iTrigger, jet->Eta()); FillH3(kHistJetPtPhi, jet->Pt(), iTrigger, jet->Phi()); if (jet->Pt() > 50.) FillH3(kHistJetEtaPhi, jet->Eta(), iTrigger, jet->Phi()); Int_t nTracks = jet->GetRefTracks()->GetEntriesFast(); for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack) { AliVParticle *track = dynamic_cast(jet->GetRefTracks()->At(iTrack)); if (track) { Float_t trackPt = track->Pt(); Float_t z = trackPt / jetPt; Float_t xi = - TMath::Log(z); FillH3(kHistJetPtTrackPt, trackPt, iTrigger, jet->Pt()); FillH3(kHistJetPtZ, z, iTrigger, jet->Pt()); FillH3(kHistJetPtXi, xi, iTrigger, jet->Pt()); } } } // integrated over all triggers FillH2(kHistJetPtNoTracks3, jet->Pt(), nJetTracks3); // limit to jets with leading track having an ITS contribution // with a hit in any SPD layer if (leadingTrack && (leadingTrack->GetFlags() & AliVTrack::kITSrefit) && (leadingTrack->HasPointOnITSLayer(0) || leadingTrack->HasPointOnITSLayer(1))) for (Int_t iTrigger = 1; iTrigger < fNoTriggers; ++iTrigger) if (IsTrigger(Trigger_t (iTrigger))) FillH1(kHistJetPtITS, jet->Pt(), iTrigger); // limit to jets having 3 tracks above 3 GeV/c if (nJetTracks3 > 2) for (Int_t iTrigger = 1; iTrigger < fNoTriggers; ++iTrigger) if (IsTrigger(Trigger_t (iTrigger))) FillH1(kHistJetPt3x3, jet->Pt(), iTrigger); } } // fill leading jet information for (Int_t iTrigger = 1; iTrigger < fNoTriggers; ++iTrigger) if (IsTrigger(Trigger_t (iTrigger))) { FillH2(kHistLeadJetPt, leadJet ? leadJet->Pt() : 0., iTrigger); if (leadJet) { Float_t leadJetPt = leadJet->Pt(); FillH3(kHistLeadJetPtEta, leadJet->Pt(), iTrigger, leadJet->Eta()); FillH3(kHistLeadJetPtPhi, leadJet->Pt(), iTrigger, leadJet->Phi()); if (leadJet->Pt() > 50.) FillH3(kHistLeadJetEtaPhi, leadJet->Eta(), iTrigger, leadJet->Phi()); Int_t nTracks = leadJet->GetRefTracks()->GetEntriesFast(); for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack) { AliVParticle *track = dynamic_cast(leadJet->GetRefTracks()->At(iTrack)); if (track) { Float_t trackPt = track->Pt(); Float_t z = trackPt / leadJetPt; Float_t xi = - TMath::Log(z); FillH3(kHistLeadJetPtTrackPt, trackPt, iTrigger, leadJet->Pt()); FillH3(kHistLeadJetPtZ, z, iTrigger, leadJet->Pt()); FillH3(kHistLeadJetPtXi, xi, iTrigger, leadJet->Pt()); } } } } } else { printf("no jet array found as branch %s\n", fJetBranchName); aodEvent->Print(); } } else { printf("no AOD event found\n"); } PostData(1, fOutputList); } void AliAnalysisTaskJetsTriggerTRD::Terminate(const Option_t * /* option */) { // actions at task termination printf("total trials: %d\n", fNTrials); } Bool_t AliAnalysisTaskJetsTriggerTRD::DetectTriggers() { fTriggerMask = 0; AliInputEventHandler *inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler(); AliVEvent::EOfflineTriggerTypes physSel = (AliVEvent::EOfflineTriggerTypes) inputHandler->IsEventSelected(); TString trgClasses = InputEvent()->GetFiredTriggerClasses(); fTrdTrg.CalcTriggers(InputEvent()); UInt_t inputMaskL1 = 0; if (AliESDEvent *esdEvent = dynamic_cast (InputEvent())) inputMaskL1 = esdEvent->GetHeader()->GetL1TriggerInputs(); else if (AliAODEvent *aodEvent = dynamic_cast (InputEvent())) inputMaskL1 = aodEvent->GetHeader()->GetL1TriggerInputs(); if (fDebug > 2) printf("trg: %8s %8s %8s %8s %8s %8s %8s (%s)\n", (physSel & AliVEvent::kAnyINT) ? "kAnyINT" : " ", (physSel & AliVEvent::kCINT5) ? "kCINT5" : " ", (physSel & AliVEvent::kINT7) ? "kINT7" : " ", (physSel & AliVEvent::kINT8) ? "kINT8" : " ", (physSel & AliVEvent::kMB) ? "kMB" : " ", (physSel & AliVEvent::kEMC7) ? "kEMC7" : " ", (physSel & AliVEvent::kTRD) ? "kTRD" : " ", trgClasses.Data() ); if (fDebug > 2) printf("trg: %8s %8s %8s %8s (%s)\n", fTrdTrg.IsFired(AliTRDTriggerAnalysis::kHJT) ? "HJT" : " ", fTrdTrg.IsFired(AliTRDTriggerAnalysis::kHSE) ? "HSE" : " ", fTrdTrg.IsFired(AliTRDTriggerAnalysis::kHQU) ? "HQU" : " ", fTrdTrg.IsFired(AliTRDTriggerAnalysis::kHEE) ? "HEE" : " ", trgClasses.Data() ); // physics selection if ((physSel & (AliVEvent::kMB))) MarkTrigger(kTrgMinBias); if ((physSel & (AliVEvent::kAnyINT | AliVEvent::kCentral | AliVEvent::kSemiCentral))) MarkTrigger(kTrgPbPb); if ((physSel & (AliVEvent::kINT7))) { MarkTrigger(kTrgInt7); if (trgClasses.Contains("WU")) { MarkTrigger(kTrgInt7_WU); if (inputMaskL1 & (1 << 9)) MarkTrigger(kTrgInt7_WUHJT); } } if ((physSel & (AliVEvent::kINT8))) MarkTrigger(kTrgInt8); if ((physSel & (AliVEvent::kEMC7)) && trgClasses.Contains("CEMC7")) MarkTrigger(kTrgEMC7); if ((physSel & (AliVEvent::kEMC8)) && trgClasses.Contains("CEMC8")) MarkTrigger(kTrgEMC8); if ((physSel & (AliVEvent::kEMCEJE))) { MarkTrigger(kTrgEMCEJE); if (trgClasses.Contains("WU")) { MarkTrigger(kTrgEMCEJE_WU); if (inputMaskL1 & (1 << 9)) MarkTrigger(kTrgEMCEJE_WUHJT); } } if ((physSel & (AliVEvent::kEMCEGA))) MarkTrigger(kTrgEMCEGA); // for the TRD-triggered events we use the classes if (trgClasses.Contains("CINT7WUHJT-")) MarkTrigger(kTrgInt7WUHJT); if (trgClasses.Contains("CINT8WUHJT-")) MarkTrigger(kTrgInt8WUHJT); if (trgClasses.Contains("CEMC7WUHJT-")) MarkTrigger(kTrgEMC7WUHJT); if (trgClasses.Contains("CEMC8WUHJT-")) MarkTrigger(kTrgEMC8WUHJT); if (HasMC()) DetectMCTriggers(); return kTRUE; } Bool_t AliAnalysisTaskJetsTriggerTRD::DetectMCTriggers() { AliMCEvent *mcEvent = MCEvent(); Int_t nTracksMC = mcEvent ? mcEvent->GetNumberOfTracks() : 0; // no. of MC tracks TList partMC; // list of particles // fill tracks passing cuts to list for (Int_t iTrackMC = 0; iTrackMC < nTracksMC; ++iTrackMC) { AliVParticle *part = mcEvent->GetTrack(iTrackMC); // only consider primaries if (!mcEvent->IsPhysicalPrimary(iTrackMC)) continue; // only consider charged particles if (part->Charge() == 0) continue; // only look at particles in eta-acceptance of the central barrel if (TMath::Abs(part->Eta()) >= 0.9) continue; partMC.Add(part); } // sort, starting from highest pt partMC.Sort(kSortAscending); Int_t nTracksInStack[4][90] = { { 0 } }; // iterate over accepted tracks TIter nextMcPart(&partMC); while (AliVParticle *part = (AliMCParticle*) (nextMcPart())) { Float_t pt = part->Pt(); Float_t eta = part->Eta(); Float_t phi = part->Phi(); Int_t phiBin = (Int_t) (phi / (2 / 18. *TMath::Pi())); // ??? use actual geometry Int_t etaBin = (Int_t) ((eta + .9) / (1.8 / 5.)); Int_t pdgCode = ((AliMCParticle*) part)->Particle()->GetPdgCode(); if (fDebug > 2) printf("pt = %f, eta = %f, phi = %f, phiBin = %d, etaBin = %d, pdgCode = %d\n", pt, eta, phi, phiBin, etaBin, pdgCode); // increment number of tracks in stack Int_t bin = phiBin * 5 + etaBin; ++nTracksInStack[0][bin]; // if minimum number reached and // the track has pt above threshold // trigger fired if ((nTracksInStack[0][bin] >= 3) && (pt >= 3.)) MarkTrigger(kTrgMC3x3Vtx); // only propagate particles that do not decay before ??? // propagate to TRD // why not use track references? Float_t rTrack = pt / .15; Float_t rTrd = 3.; Float_t phiTrd = phi + TMath::ASin(rTrd / 2 / rTrack); if (phiTrd < 0) phiTrd += 2*TMath::Pi(); else if (phiTrd > 2*TMath::Pi()) phiTrd -= 2*TMath::Pi(); Float_t etaTrd = eta; Int_t secTrd = (Int_t) (phiTrd / (2.*TMath::Pi()/18.)); Int_t stackTrd = (Int_t) ((etaTrd+0.9) / (1.8/5.)); Int_t binTrd = secTrd * 5 + stackTrd; // increment number of tracks in stack ++nTracksInStack[1][binTrd]; if (gRandom->Uniform() < fGlobalEfficiencyGTU) ++nTracksInStack[2][binTrd]; if (gRandom->Uniform() < GetEfficiencyTRD(pt, eta, phi)) ++nTracksInStack[3][binTrd]; // if minimum number reached and // the track has pt above threshold // trigger fired if (pt >= 3.) { if (nTracksInStack[1][binTrd] >= 3) MarkTrigger(kTrgMC3x3TRD); if (nTracksInStack[2][binTrd] >= 3) MarkTrigger(kTrgMC3x3TRDeff); if (nTracksInStack[3][binTrd] >= 3) MarkTrigger(kTrgMC3x3TRDeffmap); } } return kTRUE; } Bool_t AliAnalysisTaskJetsTriggerTRD::AcceptTrackMC(Int_t track) const { // only consider primaries if (!MCEvent()->IsPhysicalPrimary(track)) return kFALSE; const AliVParticle *part = MCEvent()->GetTrack(track); // only consider charged particles if (part->Charge() == 0) return kFALSE; // only look at particles in eta-acceptance of the central barrel if (TMath::Abs(part->Eta()) >= 0.9) return kFALSE; return kTRUE; } // ----- histogram management ----- TH1* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title, Int_t xbins, Float_t xmin, Float_t xmax, Int_t binType) { TString hName; hName.Form("%s_%s", fShortTaskId, hid); hName.ToLower(); TH1 *h = 0x0; if (binType == 0) h = new TH1I(hName.Data(), title, xbins, xmin, xmax); else h = new TH1F(hName.Data(), title, xbins, xmin, xmax); GetHistogram(hist) = h; fOutputList->Add(h); return h; } TH2* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title, Int_t xbins, Float_t xmin, Float_t xmax, Int_t ybins, Float_t ymin, Float_t ymax, Int_t binType) { TString hName; hName.Form("%s_%s", fShortTaskId, hid); hName.ToLower(); TH1 *h = 0x0; if (binType == 0) h = GetHistogram(hist) = new TH2I(hName.Data(), title, xbins, xmin, xmax, ybins, ymin, ymax); else h = GetHistogram(hist) = new TH2F(hName.Data(), title, xbins, xmin, xmax, ybins, ymin, ymax); fOutputList->Add(h); return (TH2*) h; } TH3* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title, Int_t xbins, Float_t xmin, Float_t xmax, Int_t ybins, Float_t ymin, Float_t ymax, Int_t zbins, Float_t zmin, Float_t zmax, Int_t binType) { TString hName; hName.Form("%s_%s", fShortTaskId, hid); hName.ToLower(); TH1 *h = 0x0; if (binType == 0) h = GetHistogram(hist) = new TH3I(hName.Data(), title, xbins, xmin, xmax, ybins, ymin, ymax, zbins, zmin, zmax); else h = GetHistogram(hist) = new TH3F(hName.Data(), title, xbins, xmin, xmax, ybins, ymin, ymax, zbins, zmin, zmax); fOutputList->Add(h); return (TH3*) h; }