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   fXsection(0.),
  46   fAvgTrials(0.),
  47   fPtHard(0.)
  48 {
  49   // default ctor
  50
  51   DefineOutput(1, TList::Class());
  52 }
  53
  54 AliAnalysisTaskJetsTriggerTRD::~AliAnalysisTaskJetsTriggerTRD()
  55 {
  56   // dtor
  57
  58 }
  59
  60 void AliAnalysisTaskJetsTriggerTRD::UserCreateOutputObjects()
  61 {
  62   // create user output objects
  63
  64   // setup list
  65   OpenFile(1);
  66   fOutputList = new TList();
  67   fOutputList->SetOwner();
  68
  69   // setup histograms
  70   TH1 *histStat = AddHistogram(ID(kHistStat), "event statistics;;counts",
  71                                kStatLast-1, .5, kStatLast-.5);
  72   histStat->GetXaxis()->SetBinLabel(ID(kStatSeen));
  73   histStat->GetXaxis()->SetBinLabel(ID(kStatTrg));
  74   histStat->GetXaxis()->SetBinLabel(ID(kStatUsed));
  75   histStat->GetXaxis()->SetBinLabel(ID(kStatEvCuts));
  76
  77   AddHistogram(ID(kHistJetPtMC), "leading jet spectrum (MC, |#eta| < 0.5);p_{T} (GeV/c);counts",
  78                fNoJetPtBins, 0., fJetPtBinMax);
  79
  80   AddHistogram(ID(kHistNoJets), "number of jets;N^{jet};counts",
  81                400, -.5, 399.5);
  82
  83   AddHistogram(ID(kHistTrackGTU), "GTU track p_{T};p_{T} (GeV/c);counts",
  84                100, 0., 25.);
  85
  86   AddHistogram(ID(kHistNPtMin),
  87                "rejection;p_{T}^{min};N_{trk};trigger",
  88                100, 0., 10.,
  89                20, 0., 20.,
  90                kTrgLast - 1, .5, kTrgLast - .5);
  91
  92   AddHistogram(ID(kHistLeadJetPt),
  93                "leading jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);counts",
  94                fNoJetPtBins, 0., fJetPtBinMax,
  95                kTrgLast - 1, .5, kTrgLast - .5);
  96   AddHistogram(ID(kHistJetPt),
  97                "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);trigger",
  98                fNoJetPtBins, 0., fJetPtBinMax,
  99                kTrgLast - 1, .5, kTrgLast - .5);
 100   AddHistogram(ID(kHistJetPtITS),
 101                "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);trigger",
 102                fNoJetPtBins, 0., fJetPtBinMax,
 103                kTrgLast - 1, .5, kTrgLast - .5);
 104   AddHistogram(ID(kHistJetPt3x3),
 105                "jet spectrum (|#eta| < 0.5);p_{T} (GeV/c);trigger",
 106                fNoJetPtBins, 0., fJetPtBinMax,
 107                kTrgLast - 1, .5, kTrgLast - .5);
 108
 109   for (Int_t iHist = kHistLeadJetPt; iHist <= kHistJetPt3x3; ++iHist) {
 110     TH1 *h = GetHistogram(Hist_t (iHist));
 111     h->GetYaxis()->SetBinLabel(ID(kTrgMinBias));
 112     h->GetYaxis()->SetBinLabel(ID(kTrgInt7));
 113     h->GetYaxis()->SetBinLabel(ID(kTrgInt8));
 114     h->GetYaxis()->SetBinLabel(ID(kTrgEMC7));
 115     h->GetYaxis()->SetBinLabel(ID(kTrgEMC8));
 116     h->GetYaxis()->SetBinLabel(ID(kTrgInt7WUHJT));
 117     h->GetYaxis()->SetBinLabel(ID(kTrgInt8WUHJT));
 118     h->GetYaxis()->SetBinLabel(ID(kTrgEMC7WUHJT));
 119     h->GetYaxis()->SetBinLabel(ID(kTrgEMC8WUHJT));
 120     h->GetYaxis()->SetBinLabel(ID(kTrgEMCEJE));
 121     h->GetYaxis()->SetBinLabel(ID(kTrgEMCEGA));
 122   }
 123
 124   AddHistogram(ID(kHistJetPtNoTracks3),
 125                "number of tracks above 3 GeV;p_{T}^{jet};no. of tracks",
 126                fNoJetPtBins, 0., fJetPtBinMax,
 127                40, -.5, 39.5);
 128
 129   PostData(1, fOutputList);
 130 }
 131
 132 Bool_t AliAnalysisTaskJetsTriggerTRD::Notify()
 133 {
 134   // actions to be taken upon notification about input file change
 135
 136   // ??? check ???
 137
 138   fXsection = 0.;
 139   fAvgTrials = 1.;
 140
 141   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 142
 143   if (tree) {
 144     TFile *curfile = tree->GetCurrentFile();
 145     if (!curfile) {
 146       Error("Notify","No current file");
 147       return kFALSE;
 148     }
 149
 150     Float_t nEntries = (Float_t) tree->GetTree()->GetEntries();
 151     if (!AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(), fXsection, fAvgTrials)) {
 152       AliError("retrieval of cross section failed");
 153       // ??? what to set as cross section?
 154     }
 155
 156     if (nEntries > 0.)
 157       fAvgTrials /= nEntries;
 158   }
 159
 160   return AliAnalysisTaskSE::Notify();
 161 }
 162
 163 void AliAnalysisTaskJetsTriggerTRD::UserExec(Option_t * /* option */)
 164 {
 165   // actual work
 166
 167   // setup pointers to input data (null if unavailable)
 168   // mcEvent:  MC input
 169   // esdEvent: ESD input
 170   // outEvent: AOD output
 171   // aodEvent: AOD input if available, otherwise AOD output
 172   AliMCEvent  *mcEvent   = this->MCEvent();
 173   AliESDEvent *esdEvent  = dynamic_cast<AliESDEvent*>(this->InputEvent()); // could also be AOD input
 174   AliAODEvent* outEvent  = this->AODEvent();
 175   AliAODEvent *aodEvent  = outEvent;
 176   if (dynamic_cast<AliAODEvent*>(this->InputEvent()))
 177     aodEvent = (AliAODEvent*) (this->InputEvent());
 178
 179   if ((fDebug > 0) && esdEvent)
 180     printf("event: %s-%06i\n", CurrentFileName(), esdEvent->GetEventNumberInFile());
 181
 182   // Int_t nTracksMC  = mcEvent ? mcEvent->GetNumberOfTracks() : 0; // no. of MC tracks
 183   // Int_t nTracks    = InputEvent()->GetNumberOfTracks(); // no. of global tracks
 184   Int_t nTracksTRD = InputEvent()->GetNumberOfTrdTracks(); // no. of GTU tracks
 185
 186   TList partMC;
 187   TList partEsd;
 188   TList partGtu;
 189
 190   Float_t leadingJetPtMC = 0.; // leading jet energy from MC information
 191   Float_t leadingJetPtRec = 0.; // leading jet energy from AOD information
 192
 193   // record number of sampled events and detect trigger contributions
 194   FillH1(kHistStat, kStatSeen);
 195   if (!DetectTriggers()) {
 196     AliError("Failed to detect the triggers");
 197     return;
 198   }
 199
 200   // only continue for events from interesting triggers
 201   if (fTriggerMask == 0)
 202     return;
 203   FillH1(kHistStat, kStatTrg);
 204
 205   // no further technical requirements for the event at the moment
 206   FillH1(kHistStat, kStatUsed);
 207
 208   // apply event cuts
 209   const AliVVertex *vtx = InputEvent()->GetPrimaryVertex();
 210   if (!vtx ||
 211       (vtx->GetNContributors() < 3.) ||
 212       (vtx->GetZ() > 10.))
 213     return;
 214   FillH1(kHistStat, kStatEvCuts);
 215
 216   // extract MC information
 217   if (mcEvent) {
 218     // check for PYTHIA event header
 219     AliGenPythiaEventHeader *pythiaHeader =
 220       dynamic_cast<AliGenPythiaEventHeader*> (mcEvent->GenEventHeader());
 221     if (!pythiaHeader) {
 222       AliWarning("MC event without PYTHIA event header!\n");
 223     }
 224     else {
 225       fPtHard  = pythiaHeader->GetPtHard();
 226       // Int_t nTrials = pythiaHeader->Trials();
 227
 228       // loop over jets from PYTHIA
 229       for (Int_t iJet = 0; iJet < pythiaHeader->NTriggerJets(); iJet++) {
 230         Float_t p[4];
 231         pythiaHeader->TriggerJet(iJet, p);
 232         TLorentzVector pJet(p);
 233         Float_t pt  = pJet.Pt();
 234         // only consider jets with |eta| < 0.5
 235         Float_t eta = pJet.Eta();
 236         if (TMath::Abs(eta) > 0.5)
 237           continue;
 238         if (pt > leadingJetPtMC)
 239           leadingJetPtMC = pt;
 240       }
 241       // fill histogram for leading jet pt spectrum
 242       FillH1(kHistJetPtMC, leadingJetPtMC, fXsection);
 243     }
 244   }
 245
 246   // loop over GTU tracks
 247   for (Int_t iTrack = 0; iTrack < nTracksTRD; ++iTrack) {
 248     AliVTrdTrack *trk = InputEvent()->GetTrdTrack(iTrack);
 249     FillH1(kHistTrackGTU, TMath::Abs(trk->Pt()));
 250     partGtu.Add(trk);
 251   }
 252   partGtu.Sort(kSortAscending);
 253
 254   Int_t nTracksPerStack[90] = { 0 };
 255   Int_t nTracksPerStackMax = 0;
 256
 257   TIter nextPartGtu(&partGtu);
 258   while (AliVTrdTrack *trdTrack = (AliVTrdTrack*) nextPartGtu()) {
 259     // count no. of tracks in stack,
 260     // check whether this number was reached before,
 261     // if not store pt^min(n),
 262     // i.e. pt of current track because of sorting
 263
 264     Int_t sec    = trdTrack->GetSector();
 265     Int_t stack  = trdTrack->GetStack();
 266
 267     if ((sec > -1) && (sec < 18) &&
 268         (stack > -1) && (stack < 5)) {
 269       ++nTracksPerStack[5*sec + stack];
 270       if (nTracksPerStack[5*sec + stack] > nTracksPerStackMax) {
 271         ++nTracksPerStackMax;
 272
 273         for (Int_t iTrigger = 1; iTrigger < kTrgLast; ++iTrigger)
 274           if (IsTrigger(Trigger_t (iTrigger)))
 275               FillH3(kHistNPtMin,
 276                      TMath::Abs(trdTrack->Pt()), nTracksPerStackMax, iTrigger);
 277       }
 278     }
 279     else {
 280       AliError(Form("Invalid sector or stack: %i %i",
 281                     sec, stack));
 282     }
 283
 284   }
 285
 286   // loop over jets from AOD event
 287   if (aodEvent) {
 288     TClonesArray *jetArray =
 289       dynamic_cast<TClonesArray*> (aodEvent->FindListObject(fJetBranchName));
 290     if (jetArray) {
 291       Int_t nJets = jetArray->GetEntriesFast();
 292       FillH1(kHistNoJets, nJets);
 293
 294       AliAODJet *leadJet = 0x0;
 295       // AliAODJet *subleadJet = 0x0;
 296
 297       for (Int_t iJet = 0; iJet < nJets; ++iJet) {
 298         AliAODJet *jet = (AliAODJet*) (*jetArray)[iJet];
 299         if (TMath::Abs(jet->Eta()) < 0.5) {
 300           // check contributing tracks
 301           Int_t nJetTracks = jet->GetRefTracks()->GetEntriesFast();
 302           Int_t nJetTracks3 = 0;
 303           AliAODTrack *leadingTrack = 0x0;
 304           for (Int_t iTrack=0; iTrack < nJetTracks; ++iTrack) {
 305             AliAODTrack *track = (AliAODTrack*) jet->GetRefTracks()->At(iTrack);
 306
 307             // count constituents above 3 GeV/c
 308             if (track->Pt() > 3.)
 309               ++nJetTracks3;
 310
 311             // find the leading track
 312             if (!leadingTrack ||
 313                 (track->Pt() > leadingTrack->Pt()))
 314               leadingTrack = track;
 315           }
 316
 317           // find leading jet
 318           if (TMath::Abs(jet->Pt()) > leadingJetPtRec) {
 319             leadingJetPtRec = TMath::Abs(jet->Pt());
 320             // subleadJet = leadJet;
 321             leadJet = jet;
 322           }
 323
 324           // jet pt spectrum
 325           for (Int_t iTrigger = 1; iTrigger < kTrgLast; ++iTrigger)
 326             if (IsTrigger(Trigger_t (iTrigger)))
 327               FillH1(kHistJetPt, jet->Pt(), iTrigger);
 328
 329           // integrated over all triggers
 330           FillH2(kHistJetPtNoTracks3, jet->Pt(), nJetTracks3);
 331
 332           // limit to jets with leading track having an ITS contribution
 333           // with a hit in any SPD layer
 334           if (leadingTrack &&
 335               (leadingTrack->GetFlags() & AliVTrack::kITSrefit) &&
 336               (leadingTrack->HasPointOnITSLayer(0) || leadingTrack->HasPointOnITSLayer(1)))
 337             for (Int_t iTrigger = 1; iTrigger < kTrgLast; ++iTrigger)
 338               if (IsTrigger(Trigger_t (iTrigger)))
 339                 FillH1(kHistJetPtITS, jet->Pt(), iTrigger);
 340
 341           // limit to jets having 3 tracks above 3 GeV/c
 342           if (nJetTracks3 > 2)
 343             for (Int_t iTrigger = 1; iTrigger < kTrgLast; ++iTrigger)
 344               if (IsTrigger(Trigger_t (iTrigger)))
 345                 FillH1(kHistJetPt3x3, jet->Pt(), iTrigger);
 346         }
 347       }
 348       // fill leading jet information
 349       for (Int_t iTrigger = 1; iTrigger < kTrgLast; ++iTrigger)
 350         if (IsTrigger(Trigger_t (iTrigger)))
 351           FillH1(kHistLeadJetPt, leadJet ? leadJet->Pt() : 0., iTrigger);
 352     }
 353     else {
 354       printf("no jet array found as branch %s\n", fJetBranchName);
 355       aodEvent->Print();
 356     }
 357   } else {
 358     printf("no AOD event found\n");
 359   }
 360
 361   PostData(1, fOutputList);
 362 }
 363
 364 void AliAnalysisTaskJetsTriggerTRD::Terminate(const Option_t * /* option */)
 365 {
 366   // actions at task termination
 367
 368 }
 369
 370 Bool_t AliAnalysisTaskJetsTriggerTRD::DetectTriggers()
 371 {
 372   fTriggerMask = 0;
 373
 374   AliInputEventHandler *inputHandler =
 375     (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 376
 377   AliVEvent::EOfflineTriggerTypes physSel = (AliVEvent::EOfflineTriggerTypes) inputHandler->IsEventSelected();
 378   TString trgClasses = InputEvent()->GetFiredTriggerClasses();
 379
 380   if (fDebug > 2)
 381     printf("trg: %8s %8s %8s %8s %8s %8s %8s (%s)\n",
 382            (physSel & AliVEvent::kAnyINT)  ? "kAnyINT"  : " ",
 383            (physSel & AliVEvent::kCINT5)   ? "kCINT5" : " ",
 384            (physSel & AliVEvent::kINT7)    ? "kINT7" : " ",
 385            (physSel & AliVEvent::kINT8)    ? "kINT8" : " ",
 386            (physSel & AliVEvent::kMB)      ? "kMB"   : " ",
 387            (physSel & AliVEvent::kEMC7)    ? "kEMC7" : " ",
 388            (physSel & AliVEvent::kTRD)     ? "kTRD" : " ",
 389            trgClasses.Data()
 390            );
 391
 392   // physics selection
 393   if ((physSel & (AliVEvent::kMB)))
 394     MarkTrigger(kTrgMinBias);
 395
 396   if ((physSel & (AliVEvent::kINT7)))
 397     MarkTrigger(kTrgInt7);
 398
 399   if ((physSel & (AliVEvent::kINT8)))
 400     MarkTrigger(kTrgInt8);
 401
 402   if ((physSel & (AliVEvent::kEMC7)) &&
 403       trgClasses.Contains("CEMC7"))
 404     MarkTrigger(kTrgEMC7);
 405
 406   if ((physSel & (AliVEvent::kEMC8)) &&
 407       trgClasses.Contains("CEMC8"))
 408     MarkTrigger(kTrgEMC8);
 409
 410   if ((physSel & (AliVEvent::kEMCEJE)))
 411     MarkTrigger(kTrgEMCEJE);
 412
 413   if ((physSel & (AliVEvent::kEMCEGA)))
 414     MarkTrigger(kTrgEMCEGA);
 415
 416   // for the TRD-triggered events we use the classes
 417   if (trgClasses.Contains("CINT7WUHJT-"))
 418     MarkTrigger(kTrgInt7WUHJT);
 419
 420   if (trgClasses.Contains("CINT8WUHJT-"))
 421     MarkTrigger(kTrgInt8WUHJT);
 422
 423   if (trgClasses.Contains("CEMC7WUHJT-"))
 424     MarkTrigger(kTrgEMC7WUHJT);
 425
 426   if (trgClasses.Contains("CEMC8WUHJT-"))
 427     MarkTrigger(kTrgEMC8WUHJT);
 428
 429   return kTRUE;
 430 }
 431
 432 // ----- histogram management -----
 433 TH1* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title,
 434                                                  Int_t xbins, Float_t xmin, Float_t xmax,
 435                                                  Int_t binType)
 436 {
 437   TString hName;
 438   hName.Form("%s_%s", fShortTaskId, hid);
 439   hName.ToLower();
 440   TH1 *h = 0x0;
 441   if (binType == 0)
 442     h = new TH1I(hName.Data(), title,
 443                  xbins, xmin, xmax);
 444   else
 445     h = new TH1F(hName.Data(), title,
 446                  xbins, xmin, xmax);
 447   GetHistogram(hist) = h;
 448   fOutputList->Add(h);
 449   return h;
 450 }
 451
 452 TH2* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title,
 453                                                  Int_t xbins, Float_t xmin, Float_t xmax,
 454                                                  Int_t ybins, Float_t ymin, Float_t ymax,
 455                                                  Int_t binType)
 456 {
 457   TString hName;
 458   hName.Form("%s_%s", fShortTaskId, hid);
 459   hName.ToLower();
 460   TH1 *h = 0x0;
 461   if (binType == 0)
 462     h = GetHistogram(hist) = new TH2I(hName.Data(), title,
 463                                      xbins, xmin, xmax,
 464                                      ybins, ymin, ymax);
 465   else
 466     h = GetHistogram(hist) = new TH2F(hName.Data(), title,
 467                                      xbins, xmin, xmax,
 468                                      ybins, ymin, ymax);
 469   fOutputList->Add(h);
 470   return (TH2*) h;
 471 }
 472
 473 TH3* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title,
 474                                                  Int_t xbins, Float_t xmin, Float_t xmax,
 475                                                  Int_t ybins, Float_t ymin, Float_t ymax,
 476                                                  Int_t zbins, Float_t zmin, Float_t zmax,
 477                                                  Int_t binType)
 478 {
 479   TString hName;
 480   hName.Form("%s_%s", fShortTaskId, hid);
 481   hName.ToLower();
 482   TH1 *h = 0x0;
 483   if (binType == 0)
 484     h = GetHistogram(hist) = new TH3I(hName.Data(), title,
 485                                      xbins, xmin, xmax,
 486                                      ybins, ymin, ymax,
 487                                      zbins, zmin, zmax);
 488   else
 489     h = GetHistogram(hist) = new TH3F(hName.Data(), title,
 490                                      xbins, xmin, xmax,
 491                                      ybins, ymin, ymax,
 492                                      zbins, zmin, zmax);
 493   fOutputList->Add(h);
 494   return (TH3*) h;
 495 }