PWGJE/UserTasks/AliAnalysisTaskJetsTriggerTRD.cxx

   1 // ROOT
   2 #include "TFile.h"
   3 #include "TList.h"
   4 #include "TH1.h"
   5 #include "TH2.h"
   6 #include "TH3.h"
   7
   8 // analysis framework
   9 #include "AliAnalysisManager.h"
  10 #include "AliInputEventHandler.h"
  11 #include "AliVEvent.h"
  12 #include "AliVTrdTrack.h"
  13
  14 // MC stuff
  15 #include "AliMCEvent.h"
  16 #include "AliGenPythiaEventHeader.h"
  17
  18 // ESD stuff
  19 #include "AliESDEvent.h"
  20 #include "AliESDInputHandler.h"
  21 #include "AliESDtrack.h"
  22 #include "AliESDTrdTrack.h"
  23 #include "AliESDTrdTracklet.h"
  24 #include "AliESDTrdTrigger.h"
  25
  26 // AOD stuff
  27 #include "AliAODEvent.h"
  28 #include "AliAODJet.h"
  29 #include "AliAODTrack.h"
  30
  31 // jet tasks
  32 #include "AliAnalysisTaskJetServices.h"
  33 #include "AliAnalysisHelperJetTasks.h"
  34
  35 #include "AliAnalysisTaskJetsTriggerTRD.h"
  36
  37 AliAnalysisTaskJetsTriggerTRD::AliAnalysisTaskJetsTriggerTRD(const char *name) :
  38   AliAnalysisTaskSE(name),
  39   fTriggerMask(0),
  40   fOutputList(),
  41   fHist(),
  42   fShortTaskId("jets_trg_trd"),
  43   fNoJetPtBins(80),
  44   fJetPtBinMax(400)
  45 {
  46   // default ctor
  47
  48   DefineOutput(1, TList::Class());
  49 }
  50
  51 AliAnalysisTaskJetsTriggerTRD::~AliAnalysisTaskJetsTriggerTRD()
  52 {
  53   // dtor
  54
  55 }
  56
  57 void AliAnalysisTaskJetsTriggerTRD::UserCreateOutputObjects()
  58 {
  59   // create user output objects
  60
  61   // setup list
  62   OpenFile(1);
  63   fOutputList = new TList();
  64   fOutputList->SetOwner();
  65
  66   // setup histograms
  67   TH1 *histStat = AddHistogram(ID(kHistStat), "event statistics;;counts",
  68                                kStatLast-1, .5, kStatLast-.5);
  69   histStat->GetXaxis()->SetBinLabel(ID(kStatSeen));
  70   histStat->GetXaxis()->SetBinLabel(ID(kStatUsed));
  71   histStat->GetXaxis()->SetBinLabel(ID(kStatMB));
  72
  73   AddHistogram(ID(kHistNoJets), "number of jets;N^{jet};counts",
  74                100, -.5, 99.5);
  75
  76   AddHistogram(ID(kHistTrackGTU), "GTU track p_{T};p_{T} (GeV/c);counts",
  77                100, 0., 25.);
  78
  79   AddHistogram(ID(kHistJetPt), "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);counts",
  80                fNoJetPtBins, 0., fJetPtBinMax);
  81   AddHistogram(ID(kHistJetPtITS), "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);counts",
  82                fNoJetPtBins, 0., fJetPtBinMax);
  83   AddHistogram(ID(kHistJetPt3x3), "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);counts",
  84                fNoJetPtBins, 0., fJetPtBinMax);
  85
  86   AddHistogram(ID(kHistJetPtEMC), "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);counts",
  87                fNoJetPtBins, 0., fJetPtBinMax);
  88   AddHistogram(ID(kHistJetPtHJT), "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);counts",
  89                fNoJetPtBins, 0., fJetPtBinMax);
  90
  91   AddHistogram(ID(kHistJetPtNoTracks3), "number of tracks above 3 GeV;p_{T}^{jet};no. of tracks",
  92                fNoJetPtBins, 0., fJetPtBinMax,
  93                20, -.5, 19.5);
  94
  95   PostData(1, fOutputList);
  96 }
  97
  98 Bool_t AliAnalysisTaskJetsTriggerTRD::Notify()
  99 {
 100   // actions to be taken upon notification about input file change
 101
 102   return AliAnalysisTaskSE::Notify();
 103 }
 104
 105 void AliAnalysisTaskJetsTriggerTRD::UserExec(Option_t * /* option */)
 106 {
 107   // actual work
 108
 109   // setup pointers to input data (null if unavailable)
 110   // mcEvent:  MC input
 111   // esdEvent: ESD input
 112   // outEvent: AOD output
 113   // aodEvent: AOD input if available, otherwise AOD output
 114   // AliMCEvent  *mcEvent   = this->MCEvent();
 115   AliESDEvent *esdEvent  = dynamic_cast<AliESDEvent*>(this->InputEvent()); // could also be AOD input
 116   AliAODEvent* outEvent  = this->AODEvent();
 117   AliAODEvent *aodEvent  = outEvent;
 118   if (dynamic_cast<AliAODEvent*>(this->InputEvent()))
 119     aodEvent = (AliAODEvent*) (this->InputEvent());
 120
 121   if ((fDebug > 0) && esdEvent)
 122     printf("event: %s-%06i\n", CurrentFileName(), esdEvent->GetEventNumberInFile());
 123
 124   // Int_t nTracksMC  = 0; // no. of MC tracks
 125   // Int_t nTracksESD = 0; // no. of global ESD tracks
 126   Int_t nTracksGTU = 0; // no. of GTU tracks
 127
 128   // Int_t nTracks[6][90]; // tracks above lower pt threshold, counted stack-wise
 129   // memset(nTracks, 0, sizeof(Int_t)*6*90);
 130   // Int_t nMax[6] = { 0 };
 131
 132   TList partMC;
 133   TList partEsd;
 134   TList partGtu;
 135
 136   // Float_t leadingJetPtMC = 0.; // leading jet energy from MC information
 137   Float_t leadingJetPtRec = 0.; // leading jet energy from AOD information
 138
 139   // record number of sampled events and detect trigger contributions
 140   FillH1(kHistStat, kStatSeen);
 141   if (!DetectTriggers()) {
 142     AliError("Failed to detect the triggers");
 143     return;
 144   }
 145
 146   if (!IsTrigger(kTrgInt))
 147     return;
 148   FillH1(kHistStat, kStatUsed);
 149
 150   for (Int_t iTrack = 0; iTrack < nTracksGTU; ++iTrack) {
 151     AliVTrdTrack *trk = 0x0;
 152     FillH1(kHistTrackGTU, trk->Pt());
 153   }
 154
 155   AliAODJet *leadJet = 0x0;
 156   AliAODJet *subleadJet = 0x0;
 157
 158   if (aodEvent) {
 159     TClonesArray *jetArray = dynamic_cast<TClonesArray*> (aodEvent->FindListObject(fJetBranchName));
 160     if (jetArray) {
 161       Int_t nJets = jetArray->GetEntriesFast();
 162       FillH1(kHistNoJets, nJets);
 163       for (Int_t iJet = 0; iJet < nJets; ++iJet) {
 164         AliAODJet *jet = (AliAODJet*) (*jetArray)[iJet];
 165
 166         // check contributing tracks
 167         Int_t nJetTracks = jet->GetRefTracks()->GetEntriesFast();
 168         Int_t nJetTracks3 = 0;
 169         Int_t iLeadingTrack = -1;
 170         Float_t ptLeadingTrack = 0.;
 171         for (Int_t iTrack=0; iTrack < nJetTracks; ++iTrack) {
 172           AliAODTrack *track = (AliAODTrack*) jet->GetRefTracks()->At(iTrack);
 173           // count tracks above 3 GeV/c
 174           if (track->Pt() > 3.)
 175             ++nJetTracks3;
 176           // find the leading track
 177           if (track->Pt() > ptLeadingTrack) {
 178             ptLeadingTrack = track->Pt();
 179             iLeadingTrack = iTrack;
 180           }
 181         }
 182
 183         // retrieve the leading track
 184         AliAODTrack *leadingTrack = (AliAODTrack*) jet->GetRefTracks()->At(iLeadingTrack);
 185
 186         if (TMath::Abs(jet->Eta()) < 0.5) {
 187
 188           // find leading jet
 189           if (TMath::Abs(jet->Pt()) > leadingJetPtRec) {
 190             leadingJetPtRec = TMath::Abs(jet->Pt());
 191             subleadJet = leadJet;
 192             leadJet = jet;
 193           }
 194
 195           FillH1(kHistJetPt, jet->Pt());
 196
 197           // discard the jet if the leading track has no ITS contribution with a hit in any SPD layer
 198           if (leadingTrack &&
 199               (leadingTrack->GetFlags() & AliVTrack::kITSrefit) &&
 200               (leadingTrack->HasPointOnITSLayer(0) || leadingTrack->HasPointOnITSLayer(1)))
 201             FillH1(kHistJetPtITS, jet->Pt());
 202
 203           // discard the jet if it does not have 3 tracks above 3 GeV/c
 204           FillH2(kHistJetPtNoTracks3, jet->Pt(), nJetTracks3);
 205           if (nJetTracks3 > 2)
 206             FillH1(kHistJetPt3x3, jet->Pt());
 207
 208           if (IsTrigger(kTrgEMC))
 209             FillH1(kHistJetPtEMC, jet->Pt());
 210           if (IsTrigger(kTrgHJT))
 211             FillH1(kHistJetPtHJT, jet->Pt());
 212         }
 213       }
 214     }
 215     else {
 216       printf("no jet array found\n");
 217     }
 218   }
 219
 220   PostData(1, fOutputList);
 221 }
 222
 223 void AliAnalysisTaskJetsTriggerTRD::Terminate(const Option_t * /* option */)
 224 {
 225   // actions at task termination
 226
 227 }
 228
 229 Bool_t AliAnalysisTaskJetsTriggerTRD::DetectTriggers()
 230 {
 231   fTriggerMask = 0;
 232
 233   AliVEvent::EOfflineTriggerTypes physSel = (AliVEvent::EOfflineTriggerTypes) ((AliInputEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 234   TString trgClasses = InputEvent()->GetFiredTriggerClasses();
 235
 236   if (fDebug > 2)
 237     printf("trg: %8s %8s %8s %8s %8s %8s %8s (%s)\n",
 238            (physSel & AliVEvent::kAnyINT)  ? "kAnyINT"  : " ",
 239            (physSel & AliVEvent::kCINT5)   ? "kCINT5" : " ",
 240            (physSel & AliVEvent::kINT7)    ? "kINT7" : " ",
 241            (physSel & AliVEvent::kINT8)    ? "kINT8" : " ",
 242            (physSel & AliVEvent::kMB)      ? "kMB"   : " ",
 243            (physSel & AliVEvent::kEMC7)    ? "kEMC7" : " ",
 244            (physSel & AliVEvent::kTRD)     ? "kTRD" : " ",
 245            trgClasses.Data()
 246            );
 247
 248   // physics selection
 249   if (physSel & (AliVEvent::kAnyINT))
 250     MarkTrigger(kTrgInt);
 251
 252   // trigger classes
 253   if (trgClasses.Contains("CINT7WUHJT") || trgClasses.Contains("CINT8WUHJT"))
 254     MarkTrigger(kTrgHJT);
 255   if (trgClasses.Contains("CEMC7WU-") || trgClasses.Contains("CEMC8WU-"))
 256     MarkTrigger(kTrgEMC);
 257
 258   return kTRUE;
 259 }
 260
 261 // ----- histogram management -----
 262 TH1* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title,
 263                                                  Int_t xbins, Float_t xmin, Float_t xmax,
 264                                                  Int_t binType)
 265 {
 266   TString hName;
 267   hName.Form("%s_%s", fShortTaskId, hid);
 268   hName.ToLower();
 269   TH1 *h = 0x0;
 270   if (binType == 0)
 271     h = new TH1I(hName.Data(), title,
 272                  xbins, xmin, xmax);
 273   else
 274     h = new TH1F(hName.Data(), title,
 275                  xbins, xmin, xmax);
 276   GetHistogram(hist) = h;
 277   fOutputList->Add(h);
 278   return h;
 279 }
 280
 281 TH2* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title,
 282                                                  Int_t xbins, Float_t xmin, Float_t xmax,
 283                                                  Int_t ybins, Float_t ymin, Float_t ymax,
 284                                                  Int_t binType)
 285 {
 286   TString hName;
 287   hName.Form("%s_%s", fShortTaskId, hid);
 288   hName.ToLower();
 289   TH1 *h = 0x0;
 290   if (binType == 0)
 291     h = GetHistogram(hist) = new TH2I(hName.Data(), title,
 292                                      xbins, xmin, xmax,
 293                                      ybins, ymin, ymax);
 294   else
 295     h = GetHistogram(hist) = new TH2F(hName.Data(), title,
 296                                      xbins, xmin, xmax,
 297                                      ybins, ymin, ymax);
 298   fOutputList->Add(h);
 299   return (TH2*) h;
 300 }
 301
 302 TH3* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title,
 303                                                  Int_t xbins, Float_t xmin, Float_t xmax,
 304                                                  Int_t ybins, Float_t ymin, Float_t ymax,
 305                                                  Int_t zbins, Float_t zmin, Float_t zmax,
 306                                                  Int_t binType)
 307 {
 308   TString hName;
 309   hName.Form("%s_%s", fShortTaskId, hid);
 310   hName.ToLower();
 311   TH1 *h = 0x0;
 312   if (binType == 0)
 313     h = GetHistogram(hist) = new TH3I(hName.Data(), title,
 314                                      xbins, xmin, xmax,
 315                                      ybins, ymin, ymax,
 316                                      zbins, zmin, zmax);
 317   else
 318     h = GetHistogram(hist) = new TH3F(hName.Data(), title,
 319                                      xbins, xmin, xmax,
 320                                      ybins, ymin, ymax,
 321                                      zbins, zmin, zmax);
 322   fOutputList->Add(h);
 323   return (TH3*) h;
 324 }