From Jochen:
authormvl <mvl@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 3 Oct 2013 13:17:43 +0000 (13:17 +0000)
committermvl <mvl@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 3 Oct 2013 13:17:43 +0000 (13:17 +0000)
- distinguish TRD-relevant trigger classes
- generalise to VEvent
- prepare for MC analysis

PWGJE/UserTasks/AliAnalysisTaskJetsTriggerTRD.cxx
PWGJE/UserTasks/AliAnalysisTaskJetsTriggerTRD.h

index b6af76a..08b8d6e 100644 (file)
@@ -41,7 +41,10 @@ AliAnalysisTaskJetsTriggerTRD::AliAnalysisTaskJetsTriggerTRD(const char *name) :
   fHist(),
   fShortTaskId("jets_trg_trd"),
   fNoJetPtBins(80),
-  fJetPtBinMax(400)
+  fJetPtBinMax(400),
+  fXsection(0.),
+  fAvgTrials(0.),
+  fPtHard(0.)
 {
   // default ctor
 
@@ -67,30 +70,59 @@ void AliAnalysisTaskJetsTriggerTRD::UserCreateOutputObjects()
   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(kStatMB));
+  histStat->GetXaxis()->SetBinLabel(ID(kStatEvCuts));
+
+  AddHistogram(ID(kHistJetPtMC), "leading jet spectrum (MC, |#eta| < 0.5);p_{T} (GeV/c);counts",
+              fNoJetPtBins, 0., fJetPtBinMax);
 
   AddHistogram(ID(kHistNoJets), "number of jets;N^{jet};counts",
-               100, -.5, 99.5);
+               400, -.5, 399.5);
 
   AddHistogram(ID(kHistTrackGTU), "GTU track p_{T};p_{T} (GeV/c);counts",
                100, 0., 25.);
 
-  AddHistogram(ID(kHistJetPt), "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);counts",
-              fNoJetPtBins, 0., fJetPtBinMax);
-  AddHistogram(ID(kHistJetPtITS), "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);counts",
-              fNoJetPtBins, 0., fJetPtBinMax);
-  AddHistogram(ID(kHistJetPt3x3), "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);counts",
-              fNoJetPtBins, 0., fJetPtBinMax);
-
-  AddHistogram(ID(kHistJetPtEMC), "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);counts",
-              fNoJetPtBins, 0., fJetPtBinMax);
-  AddHistogram(ID(kHistJetPtHJT), "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);counts",
-              fNoJetPtBins, 0., fJetPtBinMax);
+  AddHistogram(ID(kHistNPtMin),
+              "rejection;p_{T}^{min};N_{trk};trigger",
+              100, 0., 10.,
+              20, 0., 20.,
+              kTrgLast - 1, .5, kTrgLast - .5);
+
+  AddHistogram(ID(kHistLeadJetPt),
+              "leading jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);counts",
+              fNoJetPtBins, 0., fJetPtBinMax,
+              kTrgLast - 1, .5, kTrgLast - .5);
+  AddHistogram(ID(kHistJetPt),
+              "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);trigger",
+              fNoJetPtBins, 0., fJetPtBinMax,
+              kTrgLast - 1, .5, kTrgLast - .5);
+  AddHistogram(ID(kHistJetPtITS),
+              "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);trigger",
+              fNoJetPtBins, 0., fJetPtBinMax,
+              kTrgLast - 1, .5, kTrgLast - .5);
+  AddHistogram(ID(kHistJetPt3x3),
+              "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);trigger",
+              fNoJetPtBins, 0., fJetPtBinMax,
+              kTrgLast - 1, .5, kTrgLast - .5);
+
+  for (Int_t iHist = kHistLeadJetPt; iHist <= kHistJetPt3x3; ++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(kTrgInt7WUHJT));
+    h->GetYaxis()->SetBinLabel(ID(kTrgInt8WUHJT));
+    h->GetYaxis()->SetBinLabel(ID(kTrgEMC7WUHJT));
+    h->GetYaxis()->SetBinLabel(ID(kTrgEMC8WUHJT));
+  }
 
-  AddHistogram(ID(kHistJetPtNoTracks3), "number of tracks above 3 GeV;p_{T}^{jet};no. of tracks",
+  AddHistogram(ID(kHistJetPtNoTracks3),
+              "number of tracks above 3 GeV;p_{T}^{jet};no. of tracks",
                fNoJetPtBins, 0., fJetPtBinMax,
-               20, -.5, 19.5);
+               40, -.5, 39.5);
 
   PostData(1, fOutputList);
 }
@@ -99,6 +131,30 @@ Bool_t AliAnalysisTaskJetsTriggerTRD::Notify()
 {
   // actions to be taken upon notification about input file change
 
+  // ??? check ???
+
+  fXsection = 0.;
+  fAvgTrials = 1.;
+
+  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
+
+  if (tree) {
+    TFile *curfile = tree->GetCurrentFile();
+    if (!curfile) {
+      Error("Notify","No current file");
+      return kFALSE;
+    }
+
+    Float_t nEntries = (Float_t) tree->GetTree()->GetEntries();
+    if (!AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(), fXsection, fAvgTrials)) {
+      AliError("retrieval of cross section failed");
+      // ??? what to set as cross section?
+    }
+
+    if (nEntries > 0.)
+      fAvgTrials /= nEntries;
+  }
+
   return AliAnalysisTaskSE::Notify();
 }
 
@@ -111,7 +167,7 @@ void AliAnalysisTaskJetsTriggerTRD::UserExec(Option_t * /* option */)
   // esdEvent: ESD input
   // outEvent: AOD output
   // aodEvent: AOD input if available, otherwise AOD output
-  // AliMCEvent  *mcEvent   = this->MCEvent();
+  AliMCEvent  *mcEvent   = this->MCEvent();
   AliESDEvent *esdEvent  = dynamic_cast<AliESDEvent*>(this->InputEvent()); // could also be AOD input
   AliAODEvent* outEvent  = this->AODEvent();
   AliAODEvent *aodEvent  = outEvent;
@@ -121,19 +177,15 @@ void AliAnalysisTaskJetsTriggerTRD::UserExec(Option_t * /* option */)
   if ((fDebug > 0) && esdEvent)
     printf("event: %s-%06i\n", CurrentFileName(), esdEvent->GetEventNumberInFile());
 
-  // Int_t nTracksMC  = 0; // no. of MC tracks
-  // Int_t nTracksESD = 0; // no. of global ESD tracks
-  Int_t nTracksGTU = esdEvent ? esdEvent->GetNumberOfTrdTracks() : 0; // no. of GTU tracks
-
-  // Int_t nTracks[6][90]; // tracks above lower pt threshold, counted stack-wise
-  // memset(nTracks, 0, sizeof(Int_t)*6*90);
-  // Int_t nMax[6] = { 0 };
+  // 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 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
@@ -143,78 +195,165 @@ void AliAnalysisTaskJetsTriggerTRD::UserExec(Option_t * /* option */)
     return;
   }
 
-  if (!IsTrigger(kTrgInt))
+  // 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);
 
-  for (Int_t iTrack = 0; iTrack < nTracksGTU; ++iTrack) {
-    AliVTrdTrack *trk = esdEvent->GetTrdTrack(iTrack);
-    FillH1(kHistTrackGTU, trk->Pt());
+  // 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<AliGenPythiaEventHeader*> (mcEvent->GenEventHeader());
+    if (!pythiaHeader) {
+      AliWarning("MC event without PYTHIA event header!\n");
+    }
+    else {
+      fPtHard  = pythiaHeader->GetPtHard();
+      // Int_t nTrials = pythiaHeader->Trials();
+
+      // loop over jets from PYTHIA
+      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;
+      }
+      // fill histogram for leading jet pt spectrum
+      FillH1(kHistJetPtMC, leadingJetPtMC, fXsection);
+    }
   }
 
-  AliAODJet *leadJet = 0x0;
-  AliAODJet *subleadJet = 0x0;
+  // loop over GTU tracks
+  for (Int_t iTrack = 0; iTrack < nTracksTRD; ++iTrack) {
+    AliVTrdTrack *trk = InputEvent()->GetTrdTrack(iTrack);
+    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 < kTrgLast; ++iTrigger)
+         if (IsTrigger(Trigger_t (iTrigger)))
+             FillH3(kHistNPtMin,
+                    TMath::Abs(trdTrack->Pt()), nTracksPerStackMax, iTrigger);
+      }
+    }
+    else {
+      AliError(Form("Invalid sector or stack: %i %i",
+                    sec, stack));
+    }
+    
+  }
 
+  // loop over jets from AOD event
   if (aodEvent) {
-    TClonesArray *jetArray = dynamic_cast<TClonesArray*> (aodEvent->FindListObject(fJetBranchName));
+    TClonesArray *jetArray =
+      dynamic_cast<TClonesArray*> (aodEvent->FindListObject(fJetBranchName));
     if (jetArray) {
       Int_t nJets = jetArray->GetEntriesFast();
       FillH1(kHistNoJets, nJets);
-      for (Int_t iJet = 0; iJet < nJets; ++iJet) {
-        AliAODJet *jet = (AliAODJet*) (*jetArray)[iJet];
 
-        // check contributing tracks
-        Int_t nJetTracks = jet->GetRefTracks()->GetEntriesFast();
-        Int_t nJetTracks3 = 0;
-        Int_t iLeadingTrack = -1;
-        Float_t ptLeadingTrack = 0.;
-        for (Int_t iTrack=0; iTrack < nJetTracks; ++iTrack) {
-          AliAODTrack *track = (AliAODTrack*) jet->GetRefTracks()->At(iTrack);
-          // count tracks above 3 GeV/c
-          if (track->Pt() > 3.)
-            ++nJetTracks3;
-          // find the leading track
-          if (track->Pt() > ptLeadingTrack) {
-            ptLeadingTrack = track->Pt();
-            iLeadingTrack = iTrack;
-          }
-        }
-
-        // retrieve the leading track
-        AliAODTrack *leadingTrack = (AliAODTrack*) jet->GetRefTracks()->At(iLeadingTrack);
+      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;
+            // subleadJet = leadJet;
             leadJet = jet;
           }
 
-         FillH1(kHistJetPt, jet->Pt());
+         // jet pt spectrum
+         for (Int_t iTrigger = 1; iTrigger < kTrgLast; ++iTrigger)
+           if (IsTrigger(Trigger_t (iTrigger)))
+             FillH1(kHistJetPt, jet->Pt(), iTrigger);
+
+         // integrated over all triggers
+          FillH2(kHistJetPtNoTracks3, jet->Pt(), nJetTracks3);
 
-          // discard the jet if the leading track has no ITS contribution with a hit in any SPD layer
+          // 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)))
-            FillH1(kHistJetPtITS, jet->Pt());
+           for (Int_t iTrigger = 1; iTrigger < kTrgLast; ++iTrigger)
+             if (IsTrigger(Trigger_t (iTrigger)))
+               FillH1(kHistJetPtITS, jet->Pt(), iTrigger);
 
-          // discard the jet if it does not have 3 tracks above 3 GeV/c
-          FillH2(kHistJetPtNoTracks3, jet->Pt(), nJetTracks3);
+          // limit to jets having 3 tracks above 3 GeV/c
           if (nJetTracks3 > 2)
-            FillH1(kHistJetPt3x3, jet->Pt());
-
-          if (IsTrigger(kTrgEMC))
-            FillH1(kHistJetPtEMC, jet->Pt());
-          if (IsTrigger(kTrgHJT))
-            FillH1(kHistJetPtHJT, jet->Pt());
+           for (Int_t iTrigger = 1; iTrigger < kTrgLast; ++iTrigger)
+             if (IsTrigger(Trigger_t (iTrigger)))
+               FillH1(kHistJetPt3x3, jet->Pt(), iTrigger);
         }
       }
+      // fill leading jet information
+      for (Int_t iTrigger = 1; iTrigger < kTrgLast; ++iTrigger)
+       if (IsTrigger(Trigger_t (iTrigger)))
+         FillH1(kHistLeadJetPt, leadJet ? leadJet->Pt() : 0., iTrigger);
     }
     else {
-      printf("no jet array found\n");
+      printf("no jet array found as branch %s\n", fJetBranchName);
+      aodEvent->Print();
     }
+  } else {
+    printf("no AOD event found\n");
   }
 
   PostData(1, fOutputList);
@@ -230,7 +369,10 @@ Bool_t AliAnalysisTaskJetsTriggerTRD::DetectTriggers()
 {
   fTriggerMask = 0;
 
-  AliVEvent::EOfflineTriggerTypes physSel = (AliVEvent::EOfflineTriggerTypes) ((AliInputEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
+  AliInputEventHandler *inputHandler =
+    (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
+
+  AliVEvent::EOfflineTriggerTypes physSel = (AliVEvent::EOfflineTriggerTypes) inputHandler->IsEventSelected();
   TString trgClasses = InputEvent()->GetFiredTriggerClasses();
 
   if (fDebug > 2)
@@ -246,14 +388,33 @@ Bool_t AliAnalysisTaskJetsTriggerTRD::DetectTriggers()
            );
 
   // physics selection
-  if (physSel & (AliVEvent::kAnyINT))
-    MarkTrigger(kTrgInt);
-
-  // trigger classes
-  if (trgClasses.Contains("CINT7WUHJT") || trgClasses.Contains("CINT8WUHJT"))
-    MarkTrigger(kTrgHJT);
-  if (trgClasses.Contains("CEMC7WU-") || trgClasses.Contains("CEMC8WU-"))
-    MarkTrigger(kTrgEMC);
+  if ((physSel & (AliVEvent::kMB)))
+    MarkTrigger(kTrgMinBias);
+
+  if ((physSel & (AliVEvent::kINT7)))
+    MarkTrigger(kTrgInt7);
+
+  if ((physSel & (AliVEvent::kINT8)))
+    MarkTrigger(kTrgInt8);
+
+  if ((physSel & (AliVEvent::kEMC7)))
+    MarkTrigger(kTrgEMC7);
+
+  if ((physSel & (AliVEvent::kEMC8)))
+    MarkTrigger(kTrgEMC8);
+
+  // for the 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);
 
   return kTRUE;
 }
index 7e1549f..ab154d1 100644 (file)
@@ -36,10 +36,14 @@ public:
   // histograms
   enum Hist_t {
       kHistStat = 0,
+      kHistJetPtMC,
       kHistNoJets,
       kHistTrackGTU,
-      kHistJetPt, kHistJetPtITS, kHistJetPt3x3,
-      kHistJetPtEMC, kHistJetPtHJT,
+      kHistNPtMin,
+      kHistLeadJetPt,
+      kHistJetPt,
+      kHistJetPtITS,
+      kHistJetPt3x3,
       kHistJetPtNoTracks3,
       kHistLast
   };
@@ -47,19 +51,27 @@ public:
   // statistics
   enum Stat_t {
       kStatSeen = 1,
+      kStatTrg,
       kStatUsed,
-      kStatMB,
+      kStatEvCuts,
       kStatLast
   };
 
   // trigger conditions
-  enum Trigger_t { 
-      kTrgMB = 0,
-      kTrgInt,
-      kTrgInt78,
-      kTrgHJT,
-      kTrgEMC,
-      kTrgLast
+  enum Trigger_t {
+    // untriggered
+    kTrgMinBias = 1, // CINT1
+    kTrgInt7,
+    kTrgInt8,
+    kTrgEMC7,
+    kTrgEMC8,
+    // TRD jet trigger (HJT)
+    kTrgInt7WUHJT,
+    kTrgInt8WUHJT,
+    kTrgEMC7WUHJT,
+    kTrgEMC8WUHJT,
+    //
+    kTrgLast
   };      
 
 protected:
@@ -99,6 +111,11 @@ protected:
   Int_t    fNoJetPtBins;                // number of bins for jet pt
   Float_t  fJetPtBinMax;                // max jet pt (GeV) in histograms
 
+  Float_t  fXsection;                  // x-section from PYTHIA
+  Float_t  fAvgTrials;                 // ratio of PYTHIA events
+                                       // over accepted events
+  Float_t  fPtHard;                    // pt hard
+
   static const Int_t fgkStringLength = 100; // max length for the jet branch name
   char fJetBranchName[fgkStringLength];     // jet branch name