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 #include "TRandom.h"
   8
   9 // analysis framework
  10 #include "AliAnalysisManager.h"
  11 #include "AliInputEventHandler.h"
  12 #include "AliVEvent.h"
  13 #include "AliVTrdTrack.h"
  14
  15 // MC stuff
  16 #include "AliMCEvent.h"
  17 #include "AliGenPythiaEventHeader.h"
  18
  19 // ESD stuff
  20 #include "AliESDEvent.h"
  21 #include "AliESDInputHandler.h"
  22 #include "AliESDtrack.h"
  23 #include "AliESDTrdTrack.h"
  24 #include "AliESDTrdTracklet.h"
  25 #include "AliESDTrdTrigger.h"
  26
  27 // AOD stuff
  28 #include "AliAODEvent.h"
  29 #include "AliAODJet.h"
  30 #include "AliAODTrack.h"
  31
  32 // jet tasks
  33 #include "AliAnalysisTaskJetServices.h"
  34 #include "AliAnalysisHelperJetTasks.h"
  35
  36 #include "AliAnalysisTaskJetsTriggerTRD.h"
  37
  38 AliAnalysisTaskJetsTriggerTRD::AliAnalysisTaskJetsTriggerTRD(const char *name) :
  39   AliAnalysisTaskSE(name),
  40   fTriggerMask(0),
  41   fTrdTrg(),
  42   fOutputList(),
  43   fHist(),
  44   fShortTaskId("jets_trg_trd"),
  45   fNoTriggers(kTrgLast),
  46   fNoJetPtBins(80),
  47   fJetPtBinMax(400),
  48   fAvgXsection(0.),
  49   fAvgTrials(0.),
  50   fPtHard(0.),
  51   fNTrials(0),
  52   fGlobalEfficiencyGTU(.8),
  53   fGtuLabel(-1) // -3 for hw, -1821 for re-simulation
  54 {
  55   // default ctor
  56
  57   DefineOutput(1, TList::Class());
  58 }
  59
  60 AliAnalysisTaskJetsTriggerTRD::~AliAnalysisTaskJetsTriggerTRD()
  61 {
  62   // dtor
  63
  64 }
  65
  66 void AliAnalysisTaskJetsTriggerTRD::UserCreateOutputObjects()
  67 {
  68   // create user output objects
  69
  70   // setup list
  71   OpenFile(1);
  72   fOutputList = new TList();
  73   fOutputList->SetOwner();
  74
  75   // setup histograms
  76   TH1 *histStat = AddHistogram(ID(kHistStat), "event statistics;;counts",
  77                                kStatLast-1, .5, kStatLast-.5);
  78   histStat->GetXaxis()->SetBinLabel(ID(kStatSeen));
  79   histStat->GetXaxis()->SetBinLabel(ID(kStatTrg));
  80   histStat->GetXaxis()->SetBinLabel(ID(kStatUsed));
  81   histStat->GetXaxis()->SetBinLabel(ID(kStatEvCuts));
  82
  83   if (HasMC()) {
  84     fNoTriggers = kTrgMCLast;
  85     AddHistogram(ID(kHistXsection), "xsection stats",
  86                  1, 0., 1.);
  87     AddHistogram(ID(kHistPtHard), "pt hard;#hat{p}_{T} (GeV/c);counts",
  88                  fNoJetPtBins, 0., fJetPtBinMax);
  89     AddHistogram(ID(kHistJetPtMC), "leading jet spectrum (MC, |#eta| < 0.5);p_{T}^{jet} (GeV/c);counts",
  90                  fNoJetPtBins, 0., fJetPtBinMax);
  91   }
  92
  93   AddHistogram(ID(kHistJetEtaAvg), "average eta",
  94                100, -1., 1.);
  95
  96   AddHistogram(ID(kHistNoJets), "number of jets;N^{jet};counts",
  97                400, -.5, 399.5);
  98
  99   AddHistogram(ID(kHistTrackGTU), "GTU track p_{T};p_{T} (GeV/c);counts",
 100                100, 0., 25.);
 101
 102   AddHistogram(ID(kHistTrackEffGTU), "found GTU tracks;p_{T} (GeV/c);#eta;#varphi",
 103                100, 0., 25.,
 104                100, -1., 1.,
 105                100, 0., 2.*TMath::Pi());
 106   AddHistogram(ID(kHistTrackEffMC), "wanted GTU tracks;p_{T} (GeV/c);#eta;#varphi",
 107                100, 0., 25.,
 108                100, -1., 1.,
 109                100, 0., 2.*TMath::Pi());
 110
 111   AddHistogram(ID(kHistNPtMin),
 112                "rejection;p_{T}^{min};N_{trk};trigger",
 113                100, 0., 10.,
 114                20, 0., 20.,
 115                fNoTriggers - 1, .5, fNoTriggers - .5);
 116
 117   // leading jet
 118   AddHistogram(ID(kHistLeadJetPt),
 119                "leading jet spectrum (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);counts",
 120                fNoJetPtBins, 0., fJetPtBinMax,
 121                fNoTriggers - 1, .5, fNoTriggers - .5);
 122   AddHistogram(ID(kHistLeadJetPtEta),
 123                "leading jet #eta (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger;#eta",
 124                fNoJetPtBins, 0., fJetPtBinMax,
 125                fNoTriggers - 1, .5, fNoTriggers - .5,
 126                40, -.5, .5);
 127   AddHistogram(ID(kHistLeadJetPtPhi),
 128                "leading jet #phi (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger;#varphi",
 129                fNoJetPtBins, 0., fJetPtBinMax,
 130                fNoTriggers - 1, .5, fNoTriggers - .5,
 131                40, 0., 2. * TMath::Pi());
 132   AddHistogram(ID(kHistLeadJetEtaPhi),
 133                "leading jet #eta - #varphi (|#eta| < 0.5);#eta;trigger;#varphi",
 134                40, -.5, .5,
 135                fNoTriggers - 1, .5, fNoTriggers - .5,
 136                40, 0., 2. * TMath::Pi());
 137
 138   AddHistogram(ID(kHistLeadJetPtTrackPt),
 139                "leading jet pt vs track pt;p_{T} (GeV/c);trigger;p_{T}^{jet,ch} (GeV/c)",
 140                100., 0., 20.,
 141                fNoTriggers - 1, .5, fNoTriggers - .5,
 142                fNoJetPtBins, 0., fJetPtBinMax);
 143   AddHistogram(ID(kHistLeadJetPtZ),
 144                "leading jet pt vs z;z;trigger;p_{T}^{jet,ch} (GeV/c)",
 145                100., 0., 1.,
 146                fNoTriggers - 1, .5, fNoTriggers - .5,
 147                fNoJetPtBins, 0., fJetPtBinMax);
 148   AddHistogram(ID(kHistLeadJetPtXi),
 149                "leading jet pt vs #xi;#xi;trigger;p_{T}^{jet,ch} (GeV/c)",
 150                100., 0., 10.,
 151                fNoTriggers - 1, .5, fNoTriggers - .5,
 152                fNoJetPtBins, 0., fJetPtBinMax);
 153
 154   // inclusive jets
 155   AddHistogram(ID(kHistJetPt),
 156                "jet spectrum (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger",
 157                fNoJetPtBins, 0., fJetPtBinMax,
 158                fNoTriggers - 1, .5, fNoTriggers - .5);
 159   AddHistogram(ID(kHistJetPtEta),
 160                "jet #eta (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger;#eta",
 161                fNoJetPtBins, 0., fJetPtBinMax,
 162                fNoTriggers - 1, .5, fNoTriggers - .5,
 163                40, -.5, .5);
 164   AddHistogram(ID(kHistJetPtPhi),
 165                "jet #phi (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger;#varphi",
 166                fNoJetPtBins, 0., fJetPtBinMax,
 167                fNoTriggers - 1, .5, fNoTriggers - .5,
 168                40, 0., 2. * TMath::Pi());
 169   AddHistogram(ID(kHistJetEtaPhi),
 170                "jet #eta - #varphi (|#eta| < 0.5);#eta;trigger;#varphi",
 171                40, -.5, .5,
 172                fNoTriggers - 1, .5, fNoTriggers - .5,
 173                40, 0., 2. * TMath::Pi());
 174
 175   AddHistogram(ID(kHistJetPtITS),
 176                "jet spectrum (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger",
 177                fNoJetPtBins, 0., fJetPtBinMax,
 178                fNoTriggers - 1, .5, fNoTriggers - .5);
 179   AddHistogram(ID(kHistJetPt3x3),
 180                "jet spectrum (|#eta| < 0.5);p_{T}^{jet,ch} (GeV/c);trigger",
 181                fNoJetPtBins, 0., fJetPtBinMax,
 182                fNoTriggers - 1, .5, fNoTriggers - .5);
 183
 184   AddHistogram(ID(kHistJetPtTrackPt),
 185                "jet pt vs track pt;p_{T} (GeV/c);trigger;p_{T}^{jet,ch} (GeV/c)",
 186                100., 0., 20.,
 187                fNoTriggers - 1, .5, fNoTriggers - .5,
 188                fNoJetPtBins, 0., fJetPtBinMax);
 189   AddHistogram(ID(kHistJetPtZ),
 190                "jet pt vs z;z;trigger;p_{T}^{jet,ch} (GeV/c)",
 191                100., 0., 1.,
 192                fNoTriggers - 1, .5, fNoTriggers - .5,
 193                fNoJetPtBins, 0., fJetPtBinMax);
 194   AddHistogram(ID(kHistJetPtXi),
 195                "jet pt vs #xi;#xi;trigger;p_{T}^{jet,ch} (GeV/c)",
 196                100., 0., 10.,
 197                fNoTriggers - 1, .5, fNoTriggers - .5,
 198                fNoJetPtBins, 0., fJetPtBinMax);
 199
 200   for (Int_t iHist = kHistLeadJetPt; iHist <= kHistJetPtXi; ++iHist) {
 201     TH1 *h = GetHistogram(Hist_t (iHist));
 202     h->GetYaxis()->SetBinLabel(ID(kTrgMinBias));
 203     h->GetYaxis()->SetBinLabel(ID(kTrgInt7));
 204     h->GetYaxis()->SetBinLabel(ID(kTrgInt8));
 205     h->GetYaxis()->SetBinLabel(ID(kTrgEMC7));
 206     h->GetYaxis()->SetBinLabel(ID(kTrgEMC8));
 207     h->GetYaxis()->SetBinLabel(ID(kTrgPbPb));
 208     h->GetYaxis()->SetBinLabel(ID(kTrgCentral));
 209     h->GetYaxis()->SetBinLabel(ID(kTrgSemiCentral));
 210     h->GetYaxis()->SetBinLabel(ID(kTrgInt7WUHJT));
 211     h->GetYaxis()->SetBinLabel(ID(kTrgInt8WUHJT));
 212     h->GetYaxis()->SetBinLabel(ID(kTrgEMC7WUHJT));
 213     h->GetYaxis()->SetBinLabel(ID(kTrgEMC8WUHJT));
 214     h->GetYaxis()->SetBinLabel(ID(kTrgEMCEJE));
 215     h->GetYaxis()->SetBinLabel(ID(kTrgEMCEGA));
 216     h->GetYaxis()->SetBinLabel(ID(kTrgInt7_WU));
 217     h->GetYaxis()->SetBinLabel(ID(kTrgInt7_WUHJT));
 218     h->GetYaxis()->SetBinLabel(ID(kTrgEMCEJE_WU));
 219     h->GetYaxis()->SetBinLabel(ID(kTrgEMCEJE_WUHJT));
 220     if (HasMC()) {
 221       h->GetYaxis()->SetBinLabel(ID(kTrgMC3x3Vtx));
 222       h->GetYaxis()->SetBinLabel(ID(kTrgMC3x3TRD));
 223       h->GetYaxis()->SetBinLabel(ID(kTrgMC3x3TRDeff));
 224       h->GetYaxis()->SetBinLabel(ID(kTrgMC3x3TRDeffmap));
 225     }
 226   }
 227
 228   AddHistogram(ID(kHistJetPtNoTracks3),
 229                "number of tracks above 3 GeV;p_{T}^{jet,ch};no. of tracks",
 230                fNoJetPtBins, 0., fJetPtBinMax,
 231                40, -.5, 39.5);
 232
 233   PostData(1, fOutputList);
 234 }
 235
 236 Bool_t AliAnalysisTaskJetsTriggerTRD::Notify()
 237 {
 238   // actions to be taken upon notification about input file change
 239
 240   // ??? check ???
 241   // we should only count the cross section we see in the analysis,
 242   // i.e. fXsection / nTrials
 243
 244   if (HasMC()) {
 245     fAvgXsection = 0.;
 246     Float_t xSection = 0.;
 247     Float_t nTrials = 0.;
 248     Float_t nEntries = 0.;
 249
 250     TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 251
 252     if (tree) {
 253       TFile *curfile = tree->GetCurrentFile();
 254       if (!curfile) {
 255         AliError("No current file");
 256         return kFALSE;
 257       }
 258
 259       if (!AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(), xSection, nTrials)) {
 260         AliError("retrieval of cross section failed");
 261       }
 262
 263       nEntries = (Float_t) tree->GetTree()->GetEntries();
 264       if (nEntries > 0.) {
 265         fAvgTrials = nTrials / nEntries;
 266         if (nTrials > 0.)
 267           fAvgXsection = xSection / (nTrials / nEntries);
 268       }
 269     }
 270     printf("n_trials = %f\n", nTrials);
 271   }
 272
 273   return AliAnalysisTaskSE::Notify();
 274 }
 275
 276 void AliAnalysisTaskJetsTriggerTRD::UserExec(Option_t * /* option */)
 277 {
 278   // actual work
 279
 280   // setup pointers to input data (null if unavailable)
 281   // mcEvent:  MC input
 282   // esdEvent: ESD input
 283   // outEvent: AOD output
 284   // aodEvent: AOD input if available, otherwise AOD output
 285   AliMCEvent  *mcEvent   = this->MCEvent();
 286   AliESDEvent *esdEvent  = dynamic_cast<AliESDEvent*>(this->InputEvent()); // could also be AOD input
 287   AliAODEvent* outEvent  = this->AODEvent();
 288   AliAODEvent *aodEvent  = outEvent;
 289   if (dynamic_cast<AliAODEvent*>(this->InputEvent()))
 290     aodEvent = (AliAODEvent*) (this->InputEvent());
 291
 292   if ((fDebug > 0) && esdEvent)
 293     printf("event: %s-%06i\n", CurrentFileName(), esdEvent->GetEventNumberInFile());
 294
 295   Int_t nTracksMC  = mcEvent ? mcEvent->GetNumberOfTracks() : 0; // no. of MC tracks
 296   // Int_t nTracks    = InputEvent()->GetNumberOfTracks(); // no. of global tracks
 297   Int_t nTracksTRD = InputEvent()->GetNumberOfTrdTracks(); // no. of GTU tracks
 298
 299   TList partMC;
 300   TList partEsd;
 301   TList partGtu;
 302
 303   Float_t leadingJetPtMC = 0.; // leading jet energy from MC information
 304   Float_t leadingJetPtRec = 0.; // leading jet energy from AOD information
 305
 306   // record number of sampled events and detect trigger contributions
 307   FillH1(kHistStat, kStatSeen);
 308   if (!DetectTriggers()) {
 309     AliError("Failed to detect the triggers");
 310     return;
 311   }
 312
 313   // only continue for events from interesting triggers
 314   if (fTriggerMask == 0)
 315     return;
 316   FillH1(kHistStat, kStatTrg);
 317
 318   // no further technical requirements for the event at the moment
 319   FillH1(kHistStat, kStatUsed);
 320   if (HasMC()) {
 321     if (fAvgXsection == 0.) {
 322       AliError("zero cross section");
 323       return;
 324     }
 325     FillH1(kHistXsection, .5, fAvgXsection);
 326   }
 327
 328   // apply event cuts
 329   const AliVVertex *vtx = InputEvent()->GetPrimaryVertex();
 330   if (!vtx ||
 331       (vtx->GetNContributors() < 3.) ||
 332       (vtx->GetZ() > 10.))
 333     return;
 334   FillH1(kHistStat, kStatEvCuts);
 335
 336   // extract MC information
 337   if (mcEvent) {
 338     // check for PYTHIA event header
 339     AliGenPythiaEventHeader *pythiaHeader =
 340       dynamic_cast<AliGenPythiaEventHeader*> (mcEvent->GenEventHeader());
 341     if (!pythiaHeader) {
 342       AliWarning("MC event without PYTHIA event header!\n");
 343     }
 344     else {
 345       fPtHard  = pythiaHeader->GetPtHard();
 346       Int_t nTrials = pythiaHeader->Trials();
 347       fNTrials += nTrials;
 348
 349       FillH1(kHistPtHard, fPtHard);
 350
 351       // loop over jets from PYTHIA
 352       Float_t eta1 = -10.;
 353       Float_t eta2 = -10.;
 354       for (Int_t iJet = 0; iJet < pythiaHeader->NTriggerJets(); iJet++) {
 355         Float_t p[4];
 356         pythiaHeader->TriggerJet(iJet, p);
 357         TLorentzVector pJet(p);
 358         Float_t pt  = pJet.Pt();
 359         // only consider jets with |eta| < 0.5
 360         Float_t eta = pJet.Eta();
 361         if (TMath::Abs(eta) > 0.5)
 362           continue;
 363         if (pt > leadingJetPtMC)
 364           leadingJetPtMC = pt;
 365
 366         // for eta averge determination consider only jets above 60 GeV
 367         if (pt < 60.)
 368           continue;
 369         if (eta1 < -1.)
 370           eta1 = eta;
 371         else if (eta2 < -1.)
 372           eta2 = eta;
 373       }
 374       // fill histogram for leading jet pt spectrum
 375       FillH1(kHistJetPtMC, leadingJetPtMC);
 376       // fill histogram for eta average
 377       if ((eta1 > -1.) && (eta2 > -1.))
 378         FillH1(kHistJetEtaAvg, (eta1 + eta2)/2.);
 379     }
 380
 381     for (Int_t iTrack = 0; iTrack < nTracksMC; ++iTrack) {
 382       AliVParticle *part = mcEvent->GetTrack(iTrack);
 383       if (AcceptTrackMC(iTrack)) {
 384         FillH3(kHistTrackEffMC, part->Pt(), part->Eta(), part->Phi());
 385       }
 386     }
 387   }
 388
 389   // loop over GTU tracks
 390   for (Int_t iTrack = 0; iTrack < nTracksTRD; ++iTrack) {
 391     AliVTrdTrack *trk = InputEvent()->GetTrdTrack(iTrack);
 392     // printf("trk %p has pt %5.2f and label %i\n",
 393     //     trk, trk->Pt(), trk->GetLabel());
 394     if ((fGtuLabel != -1) && (trk->GetLabel() != fGtuLabel))
 395       continue;
 396     FillH1(kHistTrackGTU, TMath::Abs(trk->Pt()));
 397     partGtu.Add(trk);
 398   }
 399   partGtu.Sort(kSortAscending);
 400
 401   Int_t nTracksPerStack[90] = { 0 };
 402   Int_t nTracksPerStackMax = 0;
 403
 404   TIter nextPartGtu(&partGtu);
 405   while (AliVTrdTrack *trdTrack = (AliVTrdTrack*) nextPartGtu()) {
 406     // count no. of tracks in stack,
 407     // check whether this number was reached before,
 408     // if not store pt^min(n),
 409     // i.e. pt of current track because of sorting
 410
 411     Int_t sec    = trdTrack->GetSector();
 412     Int_t stack  = trdTrack->GetStack();
 413
 414     if ((sec > -1) && (sec < 18) &&
 415         (stack > -1) && (stack < 5)) {
 416       ++nTracksPerStack[5*sec + stack];
 417       if (nTracksPerStack[5*sec + stack] > nTracksPerStackMax) {
 418         ++nTracksPerStackMax;
 419
 420         for (Int_t iTrigger = 1; iTrigger < fNoTriggers; ++iTrigger)
 421           if (IsTrigger(Trigger_t (iTrigger)))
 422               FillH3(kHistNPtMin,
 423                      TMath::Abs(trdTrack->Pt()), nTracksPerStackMax, iTrigger);
 424       }
 425       if (HasMC()) {
 426         Int_t label = trdTrack->GetLabel();
 427         if (label > -1) {
 428           if (AcceptTrackMC(label)) {
 429             AliVParticle *part = MCEvent()->GetTrack(label);
 430             FillH3(kHistTrackEffGTU, part->Pt(), part->Eta(), part->Phi());
 431           }
 432         }
 433         else {
 434           AliWarning(Form("GTU track at %p with no label", trdTrack));
 435           const Int_t nLayers = 6;
 436           for (Int_t iLayer = 0; iLayer < nLayers; ++iLayer) {
 437             AliVTrdTracklet *trkl = trdTrack->GetTracklet(iLayer);
 438             if (trkl)
 439               AliWarning(Form("tracklet in layer %i has label %i\n",
 440                               iLayer, trkl->GetLabel()));
 441           }
 442         }
 443       }
 444     }
 445     else {
 446       AliError(Form("Invalid sector or stack: %i %i",
 447                     sec, stack));
 448     }
 449
 450   }
 451
 452   // loop over jets from AOD event
 453   if (aodEvent) {
 454     TClonesArray *jetArray =
 455       dynamic_cast<TClonesArray*> (aodEvent->FindListObject(fJetBranchName));
 456     if (jetArray) {
 457       Int_t nJets = jetArray->GetEntriesFast();
 458       FillH1(kHistNoJets, nJets);
 459
 460       AliAODJet *leadJet = 0x0;
 461       // AliAODJet *subleadJet = 0x0;
 462
 463       for (Int_t iJet = 0; iJet < nJets; ++iJet) {
 464         AliAODJet *jet = (AliAODJet*) (*jetArray)[iJet];
 465         if (TMath::Abs(jet->Eta()) < 0.5) {
 466           // check contributing tracks
 467           Int_t nJetTracks = jet->GetRefTracks()->GetEntriesFast();
 468           Int_t nJetTracks3 = 0;
 469           AliAODTrack *leadingTrack = 0x0;
 470           for (Int_t iTrack=0; iTrack < nJetTracks; ++iTrack) {
 471             AliAODTrack *track = (AliAODTrack*) jet->GetRefTracks()->At(iTrack);
 472
 473             // count constituents above 3 GeV/c
 474             if (track->Pt() > 3.)
 475               ++nJetTracks3;
 476
 477             // find the leading track
 478             if (!leadingTrack ||
 479                 (track->Pt() > leadingTrack->Pt()))
 480               leadingTrack = track;
 481           }
 482
 483           // find leading jet
 484           if (TMath::Abs(jet->Pt()) > leadingJetPtRec) {
 485             leadingJetPtRec = TMath::Abs(jet->Pt());
 486             // subleadJet = leadJet;
 487             leadJet = jet;
 488           }
 489
 490           // jet pt spectrum and
 491           // fragmentation function
 492           for (Int_t iTrigger = 1; iTrigger < fNoTriggers; ++iTrigger)
 493             if (IsTrigger(Trigger_t (iTrigger))) {
 494               Float_t jetPt = jet->Pt();
 495
 496               FillH1(kHistJetPt, jetPt, iTrigger);
 497
 498               FillH3(kHistJetPtEta, jet->Pt(), iTrigger, jet->Eta());
 499               FillH3(kHistJetPtPhi, jet->Pt(), iTrigger, jet->Phi());
 500               if (jet->Pt() > 50.)
 501                 FillH3(kHistJetEtaPhi, jet->Eta(), iTrigger, jet->Phi());
 502
 503               Int_t nTracks = jet->GetRefTracks()->GetEntriesFast();
 504               for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack) {
 505                 AliVParticle *track = dynamic_cast<AliVParticle*>(jet->GetRefTracks()->At(iTrack));
 506                 if (track) {
 507                   Float_t trackPt = track->Pt();
 508                   Float_t z = trackPt / jetPt;
 509                   Float_t xi = - TMath::Log(z);
 510                   FillH3(kHistJetPtTrackPt, trackPt, iTrigger, jet->Pt());
 511                   FillH3(kHistJetPtZ, z, iTrigger, jet->Pt());
 512                   FillH3(kHistJetPtXi, xi, iTrigger, jet->Pt());
 513                 }
 514               }
 515             }
 516
 517           // integrated over all triggers
 518           FillH2(kHistJetPtNoTracks3, jet->Pt(), nJetTracks3);
 519
 520           // limit to jets with leading track having an ITS contribution
 521           // with a hit in any SPD layer
 522           if (leadingTrack &&
 523               (leadingTrack->GetFlags() & AliVTrack::kITSrefit) &&
 524               (leadingTrack->HasPointOnITSLayer(0) || leadingTrack->HasPointOnITSLayer(1)))
 525             for (Int_t iTrigger = 1; iTrigger < fNoTriggers; ++iTrigger)
 526               if (IsTrigger(Trigger_t (iTrigger)))
 527                 FillH1(kHistJetPtITS, jet->Pt(), iTrigger);
 528
 529           // limit to jets having 3 tracks above 3 GeV/c
 530           if (nJetTracks3 > 2)
 531             for (Int_t iTrigger = 1; iTrigger < fNoTriggers; ++iTrigger)
 532               if (IsTrigger(Trigger_t (iTrigger)))
 533                 FillH1(kHistJetPt3x3, jet->Pt(), iTrigger);
 534         }
 535       }
 536
 537       // fill leading jet information
 538       for (Int_t iTrigger = 1; iTrigger < fNoTriggers; ++iTrigger)
 539         if (IsTrigger(Trigger_t (iTrigger))) {
 540           FillH2(kHistLeadJetPt, leadJet ? leadJet->Pt() : 0., iTrigger);
 541           if (leadJet) {
 542             Float_t leadJetPt = leadJet->Pt();
 543
 544             FillH3(kHistLeadJetPtEta, leadJet->Pt(), iTrigger, leadJet->Eta());
 545             FillH3(kHistLeadJetPtPhi, leadJet->Pt(), iTrigger, leadJet->Phi());
 546             if (leadJet->Pt() > 50.)
 547               FillH3(kHistLeadJetEtaPhi, leadJet->Eta(), iTrigger, leadJet->Phi());
 548
 549             Int_t nTracks = leadJet->GetRefTracks()->GetEntriesFast();
 550             for (Int_t iTrack = 0; iTrack < nTracks; ++iTrack) {
 551               AliVParticle *track = dynamic_cast<AliVParticle*>(leadJet->GetRefTracks()->At(iTrack));
 552               if (track) {
 553                 Float_t trackPt = track->Pt();
 554                 Float_t z = trackPt / leadJetPt;
 555                 Float_t xi = - TMath::Log(z);
 556                 FillH3(kHistLeadJetPtTrackPt, trackPt, iTrigger, leadJet->Pt());
 557                 FillH3(kHistLeadJetPtZ, z, iTrigger, leadJet->Pt());
 558                 FillH3(kHistLeadJetPtXi, xi, iTrigger, leadJet->Pt());
 559               }
 560             }
 561           }
 562         }
 563     }
 564     else {
 565       printf("no jet array found as branch %s\n", fJetBranchName);
 566       aodEvent->Print();
 567     }
 568   } else {
 569     printf("no AOD event found\n");
 570   }
 571
 572   PostData(1, fOutputList);
 573 }
 574
 575 void AliAnalysisTaskJetsTriggerTRD::Terminate(const Option_t * /* option */)
 576 {
 577   // actions at task termination
 578
 579   printf("total trials: %d\n", fNTrials);
 580 }
 581
 582 Bool_t AliAnalysisTaskJetsTriggerTRD::DetectTriggers()
 583 {
 584   fTriggerMask = 0;
 585
 586   AliInputEventHandler *inputHandler =
 587     (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 588
 589   AliVEvent::EOfflineTriggerTypes physSel = (AliVEvent::EOfflineTriggerTypes) inputHandler->IsEventSelected();
 590   TString trgClasses = InputEvent()->GetFiredTriggerClasses();
 591
 592   fTrdTrg.CalcTriggers(InputEvent());
 593
 594   UInt_t inputMaskL1 = 0;
 595   if (AliESDEvent *esdEvent = dynamic_cast<AliESDEvent*> (InputEvent()))
 596     inputMaskL1 = esdEvent->GetHeader()->GetL1TriggerInputs();
 597   else if (AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*> (InputEvent()))
 598     inputMaskL1 = aodEvent->GetHeader()->GetL1TriggerInputs();
 599
 600   if (fDebug > 2)
 601     printf("trg: %8s %8s %8s %8s %8s %8s %8s (%s)\n",
 602            (physSel & AliVEvent::kAnyINT)  ? "kAnyINT"  : " ",
 603            (physSel & AliVEvent::kCINT5)   ? "kCINT5" : " ",
 604            (physSel & AliVEvent::kINT7)    ? "kINT7" : " ",
 605            (physSel & AliVEvent::kINT8)    ? "kINT8" : " ",
 606            (physSel & AliVEvent::kMB)      ? "kMB"   : " ",
 607            (physSel & AliVEvent::kEMC7)    ? "kEMC7" : " ",
 608            (physSel & AliVEvent::kTRD)     ? "kTRD" : " ",
 609            trgClasses.Data()
 610            );
 611
 612   if (fDebug > 2)
 613     printf("trg: %8s %8s %8s %8s (%s)\n",
 614            fTrdTrg.IsFired(AliTRDTriggerAnalysis::kHJT) ? "HJT"  : " ",
 615            fTrdTrg.IsFired(AliTRDTriggerAnalysis::kHSE) ? "HSE"  : " ",
 616            fTrdTrg.IsFired(AliTRDTriggerAnalysis::kHQU) ? "HQU"  : " ",
 617            fTrdTrg.IsFired(AliTRDTriggerAnalysis::kHEE) ? "HEE"  : " ",
 618            trgClasses.Data()
 619            );
 620
 621   // physics selection
 622   if ((physSel & (AliVEvent::kMB)))
 623     MarkTrigger(kTrgMinBias);
 624
 625   if ((physSel & (AliVEvent::kAnyINT | AliVEvent::kCentral | AliVEvent::kSemiCentral)))
 626     MarkTrigger(kTrgPbPb);
 627
 628   if ((physSel & (AliVEvent::kINT7))) {
 629     MarkTrigger(kTrgInt7);
 630     if (trgClasses.Contains("WU")) {
 631       MarkTrigger(kTrgInt7_WU);
 632       if (inputMaskL1 & (1 << 9))
 633         MarkTrigger(kTrgInt7_WUHJT);
 634     }
 635   }
 636
 637   if ((physSel & (AliVEvent::kINT8)))
 638     MarkTrigger(kTrgInt8);
 639
 640   if ((physSel & (AliVEvent::kEMC7)) &&
 641       trgClasses.Contains("CEMC7"))
 642     MarkTrigger(kTrgEMC7);
 643
 644   if ((physSel & (AliVEvent::kEMC8)) &&
 645       trgClasses.Contains("CEMC8"))
 646     MarkTrigger(kTrgEMC8);
 647
 648   if ((physSel & (AliVEvent::kEMCEJE))) {
 649     MarkTrigger(kTrgEMCEJE);
 650     if (trgClasses.Contains("WU")) {
 651       MarkTrigger(kTrgEMCEJE_WU);
 652       if (inputMaskL1 & (1 << 9))
 653         MarkTrigger(kTrgEMCEJE_WUHJT);
 654     }
 655   }
 656
 657   if ((physSel & (AliVEvent::kEMCEGA)))
 658     MarkTrigger(kTrgEMCEGA);
 659
 660   // for the TRD-triggered events we use the classes
 661   if (trgClasses.Contains("CINT7WUHJT-"))
 662     MarkTrigger(kTrgInt7WUHJT);
 663
 664   if (trgClasses.Contains("CINT8WUHJT-"))
 665     MarkTrigger(kTrgInt8WUHJT);
 666
 667   if (trgClasses.Contains("CEMC7WUHJT-"))
 668     MarkTrigger(kTrgEMC7WUHJT);
 669
 670   if (trgClasses.Contains("CEMC8WUHJT-"))
 671     MarkTrigger(kTrgEMC8WUHJT);
 672
 673   if (HasMC())
 674     DetectMCTriggers();
 675
 676   return kTRUE;
 677 }
 678
 679 Bool_t AliAnalysisTaskJetsTriggerTRD::DetectMCTriggers()
 680 {
 681   AliMCEvent  *mcEvent   = MCEvent();
 682
 683   Int_t nTracksMC  = mcEvent ? mcEvent->GetNumberOfTracks() : 0; // no. of MC tracks
 684   TList partMC; // list of particles
 685
 686   // fill tracks passing cuts to list
 687   for (Int_t iTrackMC = 0; iTrackMC < nTracksMC; ++iTrackMC) {
 688     AliVParticle *part = mcEvent->GetTrack(iTrackMC);
 689
 690     // only consider primaries
 691     if (!mcEvent->IsPhysicalPrimary(iTrackMC))
 692       continue;
 693
 694     // only consider charged particles
 695     if (part->Charge() == 0)
 696       continue;
 697
 698     // only look at particles in eta-acceptance of the central barrel
 699     if (TMath::Abs(part->Eta()) >= 0.9)
 700       continue;
 701
 702     partMC.Add(part);
 703   }
 704
 705   // sort, starting from highest pt
 706   partMC.Sort(kSortAscending);
 707
 708   Int_t nTracksInStack[4][90] = { { 0 } };
 709
 710   // iterate over accepted tracks
 711   TIter nextMcPart(&partMC);
 712   while (AliVParticle *part = (AliMCParticle*) (nextMcPart())) {
 713     Float_t pt  = part->Pt();
 714     Float_t eta = part->Eta();
 715     Float_t phi = part->Phi();
 716
 717     Int_t phiBin = (Int_t) (phi / (2 / 18. *TMath::Pi()));
 718     // ??? use actual geometry
 719     Int_t etaBin = (Int_t) ((eta + .9) / (1.8 / 5.));
 720
 721     Int_t pdgCode = ((AliMCParticle*) part)->Particle()->GetPdgCode();
 722     if (fDebug > 2)
 723       printf("pt = %f, eta = %f, phi = %f, phiBin = %d, etaBin = %d, pdgCode = %d\n",
 724              pt, eta, phi, phiBin, etaBin, pdgCode);
 725
 726     // increment number of tracks in stack
 727     Int_t bin = phiBin * 5 + etaBin;
 728     ++nTracksInStack[0][bin];
 729
 730     // if minimum number reached and
 731     // the track has pt above threshold
 732     // trigger fired
 733     if ((nTracksInStack[0][bin] >= 3) &&
 734         (pt >= 3.))
 735       MarkTrigger(kTrgMC3x3Vtx);
 736
 737     // only propagate particles that do not decay before ???
 738
 739     // propagate to TRD
 740     // why not use track references?
 741     Float_t rTrack = pt / .15;
 742     Float_t rTrd   = 3.;
 743     Float_t phiTrd = phi + TMath::ASin(rTrd / 2 / rTrack);
 744     if (phiTrd < 0)
 745       phiTrd += 2*TMath::Pi();
 746     else if (phiTrd > 2*TMath::Pi())
 747       phiTrd -= 2*TMath::Pi();
 748     Float_t etaTrd = eta;
 749
 750     Int_t secTrd   = (Int_t) (phiTrd / (2.*TMath::Pi()/18.));
 751     Int_t stackTrd = (Int_t) ((etaTrd+0.9) / (1.8/5.));
 752     Int_t binTrd   = secTrd * 5 + stackTrd;
 753
 754     // increment number of tracks in stack
 755     ++nTracksInStack[1][binTrd];
 756     if (gRandom->Uniform() < fGlobalEfficiencyGTU)
 757       ++nTracksInStack[2][binTrd];
 758     if (gRandom->Uniform() < GetEfficiencyTRD(pt, eta, phi))
 759       ++nTracksInStack[3][binTrd];
 760
 761     // if minimum number reached and
 762     // the track has pt above threshold
 763     // trigger fired
 764     if (pt >= 3.) {
 765       if (nTracksInStack[1][binTrd] >= 3)
 766         MarkTrigger(kTrgMC3x3TRD);
 767       if (nTracksInStack[2][binTrd] >= 3)
 768         MarkTrigger(kTrgMC3x3TRDeff);
 769       if (nTracksInStack[3][binTrd] >= 3)
 770         MarkTrigger(kTrgMC3x3TRDeffmap);
 771     }
 772   }
 773
 774   return kTRUE;
 775 }
 776
 777 Bool_t AliAnalysisTaskJetsTriggerTRD::AcceptTrackMC(Int_t track) const
 778 {
 779   // only consider primaries
 780   if (!MCEvent()->IsPhysicalPrimary(track))
 781     return kFALSE;
 782
 783   const AliVParticle *part = MCEvent()->GetTrack(track);
 784
 785   // only consider charged particles
 786   if (part->Charge() == 0)
 787     return kFALSE;
 788
 789   // only look at particles in eta-acceptance of the central barrel
 790   if (TMath::Abs(part->Eta()) >= 0.9)
 791     return kFALSE;
 792
 793   return kTRUE;
 794 }
 795
 796
 797 // ----- histogram management -----
 798 TH1* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title,
 799                                                  Int_t xbins, Float_t xmin, Float_t xmax,
 800                                                  Int_t binType)
 801 {
 802   TString hName;
 803   hName.Form("%s_%s", fShortTaskId, hid);
 804   hName.ToLower();
 805   TH1 *h = 0x0;
 806   if (binType == 0)
 807     h = new TH1I(hName.Data(), title,
 808                  xbins, xmin, xmax);
 809   else
 810     h = new TH1F(hName.Data(), title,
 811                  xbins, xmin, xmax);
 812   GetHistogram(hist) = h;
 813   fOutputList->Add(h);
 814   return h;
 815 }
 816
 817 TH2* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title,
 818                                                  Int_t xbins, Float_t xmin, Float_t xmax,
 819                                                  Int_t ybins, Float_t ymin, Float_t ymax,
 820                                                  Int_t binType)
 821 {
 822   TString hName;
 823   hName.Form("%s_%s", fShortTaskId, hid);
 824   hName.ToLower();
 825   TH1 *h = 0x0;
 826   if (binType == 0)
 827     h = GetHistogram(hist) = new TH2I(hName.Data(), title,
 828                                      xbins, xmin, xmax,
 829                                      ybins, ymin, ymax);
 830   else
 831     h = GetHistogram(hist) = new TH2F(hName.Data(), title,
 832                                      xbins, xmin, xmax,
 833                                      ybins, ymin, ymax);
 834   fOutputList->Add(h);
 835   return (TH2*) h;
 836 }
 837
 838 TH3* AliAnalysisTaskJetsTriggerTRD::AddHistogram(Hist_t hist, const char *hid, TString title,
 839                                                  Int_t xbins, Float_t xmin, Float_t xmax,
 840                                                  Int_t ybins, Float_t ymin, Float_t ymax,
 841                                                  Int_t zbins, Float_t zmin, Float_t zmax,
 842                                                  Int_t binType)
 843 {
 844   TString hName;
 845   hName.Form("%s_%s", fShortTaskId, hid);
 846   hName.ToLower();
 847   TH1 *h = 0x0;
 848   if (binType == 0)
 849     h = GetHistogram(hist) = new TH3I(hName.Data(), title,
 850                                      xbins, xmin, xmax,
 851                                      ybins, ymin, ymax,
 852                                      zbins, zmin, zmax);
 853   else
 854     h = GetHistogram(hist) = new TH3F(hName.Data(), title,
 855                                      xbins, xmin, xmax,
 856                                      ybins, ymin, ymax,
 857                                      zbins, zmin, zmax);
 858   fOutputList->Add(h);
 859   return (TH3*) h;
 860 }