PWGGA/EMCALTasks/AliAnalysisTaskSAJF.cxx

   1 // $Id$
   2 //
   3 // Jet finder analysis task (S.Aiola).
   4 //
   5 //
   6
   7 #include <TChain.h>
   8 #include <TClonesArray.h>
   9 #include <TH1F.h>
  10 #include <TH2F.h>
  11 #include <TList.h>
  12 #include <TLorentzVector.h>
  13 #include <TParticle.h>
  14 #include <TRandom3.h>
  15
  16 #include "AliAnalysisManager.h"
  17 #include "AliCentrality.h"
  18 #include "AliVCluster.h"
  19 #include "AliVTrack.h"
  20 #include "AliEmcalJet.h"
  21 #include "AliVEventHandler.h"
  22 #include "AliLog.h"
  23
  24 #include "AliAnalysisTaskSAJF.h"
  25
  26 ClassImp(AliAnalysisTaskSAJF)
  27
  28 //________________________________________________________________________
  29 AliAnalysisTaskSAJF::AliAnalysisTaskSAJF() :
  30   AliAnalysisTaskSE("AliAnalysisTaskSAJF"),
  31   fAnaType(kEMCAL),
  32   fMinEta(-1),
  33   fMaxEta(1),
  34   fMinPhi(0),
  35   fMaxPhi(2 * TMath::Pi()),
  36   fJetRadius(0.4),
  37   fOutput(0),
  38   fTracksName("Tracks"),
  39   fCaloName("CaloClusters"),
  40   fJetsName("Jets"),
  41   fKtJetsName("KtJets"),
  42   fEmbJetsName("EmbJets"),
  43   fTrgClusName("ClustersL1GAMMAFEE"),
  44   fTracks(0),
  45   fCaloClusters(0),
  46   fJets(0),
  47   fKtJets(0),
  48   fEmbJets(0),
  49   fTrgClusters(0),
  50   fCent(0),
  51   fCentBin(-1),
  52   fHistCentrality(0),
  53   fHistJetPhiEta(0),
  54   fHistRhoPartVSleadJetPt(0),
  55   fHistMedKtVSRhoPart(0),
  56   fHistRCPtVSRhoPart(0),
  57   fNbins(500),
  58   fMinPt(0),
  59   fMaxPt(250)
  60 {
  61   // Default constructor.
  62
  63   for (Int_t i = 0; i < 4; i++) {
  64     fHistJetsPt[i] = 0;
  65     fHistJetsNEF[i] = 0;
  66     fHistJetsZ[i] = 0;
  67     fHistLeadingJetPt[i] = 0;
  68     fHist2LeadingJetPt[i] = 0;
  69     fHistTracksPtLJ[i] = 0;
  70     fHistClusEtLJ[i] = 0;
  71     fHistTracksPtBkg[i] = 0;
  72     fHistClusEtBkg[i] = 0;
  73     fHistBkgClusPhiCorr[i] = 0;
  74     fHistBkgTracksPhiCorr[i] = 0;
  75     fHistBkgClusMeanRho[i] = 0;
  76     fHistBkgTracksMeanRho[i] = 0;
  77     fHistBkgLJetPhiCorr[i] = 0;
  78     fHistMedianPtKtJet[i] = 0;
  79     fHistDeltaPtRC[i] = 0;
  80     fHistDeltaPtRCExLJ[i] = 0;
  81     fHistRCPt[i] = 0;
  82     fHistRCPtExLJ[i] = 0;
  83     fHistEmbJets[i] = 0;
  84     fHistEmbPart[i] = 0;
  85     fHistDeltaPtMedKtEmb[i] = 0;
  86   }
  87
  88   // Output slot #1 writes into a TH1 container
  89   DefineOutput(1, TList::Class());
  90 }
  91
  92 //________________________________________________________________________
  93 AliAnalysisTaskSAJF::AliAnalysisTaskSAJF(const char *name) :
  94   AliAnalysisTaskSE(name),
  95   fAnaType(kEMCAL),
  96   fMinEta(-1),
  97   fMaxEta(1),
  98   fMinPhi(0),
  99   fMaxPhi(2 * TMath::Pi()),
 100   fJetRadius(0.4),
 101   fOutput(0),
 102   fTracksName("Tracks"),
 103   fCaloName("CaloClusters"),
 104   fJetsName("Jets"),
 105   fKtJetsName("KtJets"),
 106   fEmbJetsName("EmbJets"),
 107   fTrgClusName("ClustersL1GAMMAFEE"),
 108   fTracks(0),
 109   fCaloClusters(0),
 110   fJets(0),
 111   fKtJets(0),
 112   fEmbJets(0),
 113   fTrgClusters(0),
 114   fCent(0),
 115   fCentBin(-1),
 116   fHistCentrality(0),
 117   fHistJetPhiEta(0),
 118   fHistRhoPartVSleadJetPt(0),
 119   fHistMedKtVSRhoPart(0),
 120   fHistRCPtVSRhoPart(0),
 121   fNbins(500),
 122   fMinPt(0),
 123   fMaxPt(250)
 124 {
 125   // Standard constructor.
 126
 127   for (Int_t i = 0; i < 4; i++) {
 128     fHistJetsPt[i] = 0;
 129     fHistJetsNEF[i] = 0;
 130     fHistJetsZ[i] = 0;
 131     fHistLeadingJetPt[i] = 0;
 132     fHist2LeadingJetPt[i] = 0;
 133     fHistTracksPtLJ[i] = 0;
 134     fHistClusEtLJ[i] = 0;
 135     fHistTracksPtBkg[i] = 0;
 136     fHistClusEtBkg[i] = 0;
 137     fHistBkgClusPhiCorr[i] = 0;
 138     fHistBkgTracksPhiCorr[i] = 0;
 139     fHistBkgClusMeanRho[i] = 0;
 140     fHistBkgTracksMeanRho[i] = 0;
 141     fHistBkgLJetPhiCorr[i] = 0;
 142     fHistMedianPtKtJet[i] = 0;
 143     fHistDeltaPtRC[i] = 0;
 144     fHistDeltaPtRCExLJ[i] = 0;
 145     fHistRCPt[i] = 0;
 146     fHistRCPtExLJ[i] = 0;
 147     fHistEmbJets[i] = 0;
 148     fHistEmbPart[i] = 0;
 149     fHistDeltaPtMedKtEmb[i] = 0;
 150   }
 151
 152   // Output slot #1 writes into a TH1 container
 153   DefineOutput(1, TList::Class());
 154 }
 155
 156 //________________________________________________________________________
 157 AliAnalysisTaskSAJF::~AliAnalysisTaskSAJF()
 158 {
 159   // Destructor
 160 }
 161
 162 //________________________________________________________________________
 163 void AliAnalysisTaskSAJF::UserCreateOutputObjects()
 164 {
 165   // Create histograms
 166
 167   fOutput = new TList();
 168   fOutput->SetOwner();  // IMPORTANT!
 169
 170   fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", fNbins, 0, 100);
 171   fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
 172   fHistCentrality->GetYaxis()->SetTitle("counts");
 173   fOutput->Add(fHistCentrality);
 174
 175   fHistJetPhiEta = new TH2F("fHistJetPhiEta","Phi-Eta distribution of jets", 20, -2, 2, 32, 0, 6.4);
 176   fHistJetPhiEta->GetXaxis()->SetTitle("Eta");
 177   fHistJetPhiEta->GetYaxis()->SetTitle("Phi");
 178   fOutput->Add(fHistJetPhiEta);
 179
 180   fHistRhoPartVSleadJetPt = new TH2F("fHistRhoPartVSleadJetPt","fHistRhoPartVSleadJetPt", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
 181   fHistRhoPartVSleadJetPt->GetXaxis()->SetTitle("#rho [GeV]");
 182   fHistRhoPartVSleadJetPt->GetYaxis()->SetTitle("Leading jet energy [GeV]");
 183   fOutput->Add(fHistRhoPartVSleadJetPt);
 184
 185   fHistMedKtVSRhoPart = new TH2F("fHistMedKtVSRhoPart","fHistMedKtVSRhoPart", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
 186   fHistMedKtVSRhoPart->GetXaxis()->SetTitle("median kt P_{T} [GeV]");
 187   fHistMedKtVSRhoPart->GetYaxis()->SetTitle("#rho [GeV]");
 188   fOutput->Add(fHistMedKtVSRhoPart);
 189
 190   fHistRCPtVSRhoPart = new TH2F("fHistRCPtVSRhoPart","fHistRCPtVSRhoPart", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
 191   fHistRCPtVSRhoPart->GetXaxis()->SetTitle("rigid cone P_{T} [GeV]");
 192   fHistRCPtVSRhoPart->GetYaxis()->SetTitle("#rho [GeV]");
 193   fOutput->Add(fHistRCPtVSRhoPart);
 194
 195   TString histname;
 196
 197   for (Int_t i = 0; i < 4; i++) {
 198     histname = "fHistJetsPt_";
 199     histname += i;
 200     fHistJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 201     fHistJetsPt[i]->GetXaxis()->SetTitle("E [GeV]");
 202     fHistJetsPt[i]->GetYaxis()->SetTitle("counts");
 203     fOutput->Add(fHistJetsPt[i]);
 204
 205     histname = "fHistJetsNEF_";
 206     histname += i;
 207     fHistJetsNEF[i] = new TH1F(histname.Data(), histname.Data(), fNbins, 0, 1.2);
 208     fHistJetsNEF[i]->GetXaxis()->SetTitle("NEF");
 209     fHistJetsNEF[i]->GetYaxis()->SetTitle("counts");
 210     fOutput->Add(fHistJetsNEF[i]);
 211
 212     histname = "fHistJetsZ_";
 213     histname += i;
 214     fHistJetsZ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, 0, 1.2);
 215     fHistJetsZ[i]->GetXaxis()->SetTitle("Z");
 216     fHistJetsZ[i]->GetYaxis()->SetTitle("counts");
 217     fOutput->Add(fHistJetsZ[i]);
 218
 219     histname = "fHistLeadingJetPt_";
 220     histname += i;
 221     fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 222     fHistLeadingJetPt[i]->GetXaxis()->SetTitle("P_{T} [GeV]");
 223     fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
 224     fOutput->Add(fHistLeadingJetPt[i]);
 225
 226     histname = "fHistTracksPtLJ_";
 227     histname += i;
 228     fHistTracksPtLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 229     fHistTracksPtLJ[i]->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 230     fHistTracksPtLJ[i]->GetYaxis()->SetTitle("counts");
 231     fOutput->Add(fHistTracksPtLJ[i]);
 232
 233     histname = "fHistClusEtLJ_";
 234     histname += i;
 235     fHistClusEtLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 236     fHistClusEtLJ[i]->GetXaxis()->SetTitle("E_{T} [GeV]");
 237     fHistClusEtLJ[i]->GetYaxis()->SetTitle("counts");
 238     fOutput->Add(fHistClusEtLJ[i]);
 239
 240     histname = "fHistTracksPtBkg_";
 241     histname += i;
 242     fHistTracksPtBkg[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 243     fHistTracksPtBkg[i]->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 244     fHistTracksPtBkg[i]->GetYaxis()->SetTitle("counts");
 245     fOutput->Add(fHistTracksPtBkg[i]);
 246
 247     histname = "fHistClusEtBkg_";
 248     histname += i;
 249     fHistClusEtBkg[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 250     fHistClusEtBkg[i]->GetXaxis()->SetTitle("E_{T} [GeV]");
 251     fHistClusEtBkg[i]->GetYaxis()->SetTitle("counts");
 252     fOutput->Add(fHistClusEtBkg[i]);
 253
 254     histname = "fHistBkgClusPhiCorr_";
 255     histname += i;
 256     fHistBkgClusPhiCorr[i] = new TH1F(histname.Data(), histname.Data(), 128, -1.6, 4.8);
 257     fHistBkgClusPhiCorr[i]->GetXaxis()->SetTitle("#Delta#phi");
 258     fHistBkgClusPhiCorr[i]->GetYaxis()->SetTitle("counts");
 259     fOutput->Add(fHistBkgClusPhiCorr[i]);
 260
 261     histname = "fHistBkgTracksPhiCorr_";
 262     histname += i;
 263     fHistBkgTracksPhiCorr[i] = new TH1F(histname.Data(), histname.Data(), 128, -1.6, 4.8);
 264     fHistBkgTracksPhiCorr[i]->GetXaxis()->SetTitle("#Delta#phi");
 265     fHistBkgTracksPhiCorr[i]->GetYaxis()->SetTitle("counts");
 266     fOutput->Add(fHistBkgTracksPhiCorr[i]);
 267
 268     histname = "fHistBkgClusMeanRho_";
 269     histname += i;
 270     fHistBkgClusMeanRho[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 271     fHistBkgClusMeanRho[i]->GetXaxis()->SetTitle("#rho [GeV]");
 272     fHistBkgClusMeanRho[i]->GetYaxis()->SetTitle("counts");
 273     fOutput->Add(fHistBkgClusMeanRho[i]);
 274
 275     histname = "fHistBkgTracksMeanRho_";
 276     histname += i;
 277     fHistBkgTracksMeanRho[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 278     fHistBkgTracksMeanRho[i]->GetXaxis()->SetTitle("#rho [GeV/c]");
 279     fHistBkgTracksMeanRho[i]->GetYaxis()->SetTitle("counts");
 280     fOutput->Add(fHistBkgTracksMeanRho[i]);
 281
 282     histname = "fHistBkgLJetPhiCorr_";
 283     histname += i;
 284     fHistBkgLJetPhiCorr[i] = new TH1F(histname.Data(), histname.Data(), 128, -1.6, 4.8);
 285     fHistBkgLJetPhiCorr[i]->GetXaxis()->SetTitle("#Delta#phi");
 286     fHistBkgLJetPhiCorr[i]->GetYaxis()->SetTitle("counts");
 287     fOutput->Add(fHistBkgLJetPhiCorr[i]);
 288
 289     histname = "fHistMedianPtKtJet_";
 290     histname += i;
 291     fHistMedianPtKtJet[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 292     fHistMedianPtKtJet[i]->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 293     fHistMedianPtKtJet[i]->GetYaxis()->SetTitle("counts");
 294     fOutput->Add(fHistMedianPtKtJet[i]);
 295
 296     histname = "fHistDeltaPtRC_";
 297     histname += i;
 298     fHistDeltaPtRC[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2, fMaxPt / 2);
 299     fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#deltaP_{T} [GeV/c]");
 300     fHistDeltaPtRC[i]->GetYaxis()->SetTitle("counts");
 301     fOutput->Add(fHistDeltaPtRC[i]);
 302
 303     histname = "fHistDeltaPtRCExLJ_";
 304     histname += i;
 305     fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2, fMaxPt / 2);
 306     fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#deltaP_{T} [GeV/c]");
 307     fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
 308     fOutput->Add(fHistDeltaPtRCExLJ[i]);
 309
 310     histname = "fHistRCPt_";
 311     histname += i;
 312     fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 313     fHistRCPt[i]->GetXaxis()->SetTitle("rigid cone P_{T} [GeV/c]");
 314     fHistRCPt[i]->GetYaxis()->SetTitle("counts");
 315     fOutput->Add(fHistRCPt[i]);
 316
 317     histname = "fHistRCPtExLJ_";
 318     histname += i;
 319     fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 320     fHistRCPtExLJ[i]->GetXaxis()->SetTitle("rigid cone P_{T} [GeV/c]");
 321     fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
 322     fOutput->Add(fHistRCPtExLJ[i]);
 323
 324     histname = "fHistEmbJets_";
 325     histname += i;
 326     fHistEmbJets[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 327     fHistEmbJets[i]->GetXaxis()->SetTitle("embedded jet P_{T} [GeV/c]");
 328     fHistEmbJets[i]->GetYaxis()->SetTitle("counts");
 329     fOutput->Add(fHistEmbJets[i]);
 330
 331     histname = "fHistEmbPart_";
 332     histname += i;
 333     fHistEmbPart[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 334     fHistEmbPart[i]->GetXaxis()->SetTitle("embedded particle P_{T} [GeV/c]");
 335     fHistEmbPart[i]->GetYaxis()->SetTitle("counts");
 336     fOutput->Add(fHistEmbPart[i]);
 337
 338     histname = "fHistDeltaPtMedKtEmb_";
 339     histname += i;
 340     fHistDeltaPtMedKtEmb[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2, fMaxPt / 2);
 341     fHistDeltaPtMedKtEmb[i]->GetXaxis()->SetTitle("#deltaP_{T} [GeV/c]");
 342     fHistDeltaPtMedKtEmb[i]->GetYaxis()->SetTitle("counts");
 343     fOutput->Add(fHistDeltaPtMedKtEmb[i]);
 344   }
 345
 346   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 347 }
 348
 349 //________________________________________________________________________
 350 void AliAnalysisTaskSAJF::RetrieveEventObjects()
 351 {
 352   fCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
 353   if (!fCaloClusters) {
 354     AliWarning(Form("Could not retrieve clusters!"));
 355   }
 356
 357   fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
 358   if (!fTracks) {
 359     AliWarning(Form("Could not retrieve tracks!"));
 360   }
 361
 362   fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
 363   if (!fJets) {
 364     AliWarning(Form("Could not retrieve jets!"));
 365   }
 366
 367   if (strcmp(fKtJetsName,"")) {
 368     fKtJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fKtJetsName));
 369     if (!fKtJets) {
 370       AliWarning(Form("Could not retrieve Kt jets!"));
 371     }
 372   }
 373
 374   if (strcmp(fEmbJetsName,"")) {
 375     fEmbJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
 376     if (!fEmbJets) {
 377       AliWarning(Form("Could not retrieve emb jets!"));
 378     }
 379   }
 380
 381   if (strcmp(fTrgClusName,"")) {
 382     fTrgClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrgClusName));
 383     if (!fTrgClusters) {
 384       AliWarning(Form("Could not retrieve trigger clusters!"));
 385     }
 386   }
 387
 388   fCent = InputEvent()->GetCentrality();
 389   if (fCent) {
 390     Float_t cent = fCent->GetCentralityPercentile("V0M");
 391     if      (cent >=  0 && cent <   10) fCentBin = 0;
 392     else if (cent >= 10 && cent <   30) fCentBin = 1;
 393     else if (cent >= 30 && cent <   50) fCentBin = 2;
 394     else if (cent >= 50 && cent <= 100) fCentBin = 3;
 395     else {
 396       AliWarning(Form("Negative centrality: %f. Assuming 99", cent));
 397       fCentBin = 3;
 398     }
 399   }
 400   else {
 401     AliWarning(Form("Could not retrieve centrality information! Assuming 99"));
 402     fCentBin = 3;
 403   }
 404 }
 405
 406 //________________________________________________________________________
 407 AliVTrack* AliAnalysisTaskSAJF::GetTrack(const Int_t i) const
 408 {
 409   if (fTracks)
 410     return dynamic_cast<AliVTrack*>(fTracks->At(i));
 411   else
 412     return 0;
 413 }
 414
 415 //________________________________________________________________________
 416 Int_t AliAnalysisTaskSAJF::GetNumberOfTracks() const
 417 {
 418   if (fTracks)
 419     return fTracks->GetEntriesFast();
 420   else
 421     return 0;
 422 }
 423
 424 //________________________________________________________________________
 425 AliVCluster* AliAnalysisTaskSAJF::GetCaloCluster(const Int_t i) const
 426 {
 427   if (fCaloClusters)
 428     return dynamic_cast<AliVCluster*>(fCaloClusters->At(i));
 429   else
 430     return 0;
 431 }
 432
 433 //________________________________________________________________________
 434 Int_t AliAnalysisTaskSAJF::GetNumberOfCaloClusters() const
 435 {
 436   if (fCaloClusters)
 437     return fCaloClusters->GetEntriesFast();
 438   else
 439     return 0;
 440 }
 441
 442 //________________________________________________________________________
 443 AliEmcalJet* AliAnalysisTaskSAJF::GetJet(const Int_t i) const
 444 {
 445   if (fJets)
 446     return dynamic_cast<AliEmcalJet*>(fJets->At(i));
 447   else
 448     return 0;
 449 }
 450
 451 //________________________________________________________________________
 452 Int_t AliAnalysisTaskSAJF::GetNumberOfJets() const
 453 {
 454   if (fJets)
 455     return fJets->GetEntriesFast();
 456   else
 457     return 0;
 458 }
 459
 460 //________________________________________________________________________
 461 AliEmcalJet* AliAnalysisTaskSAJF::GetKtJet(const Int_t i) const
 462 {
 463   if (fKtJets)
 464     return dynamic_cast<AliEmcalJet*>(fKtJets->At(i));
 465   else
 466     return 0;
 467 }
 468
 469 //________________________________________________________________________
 470 Int_t AliAnalysisTaskSAJF::GetNumberOfKtJets() const
 471 {
 472   if (fKtJets)
 473     return fKtJets->GetEntriesFast();
 474   else
 475     return 0;
 476 }
 477
 478 //________________________________________________________________________
 479 AliEmcalJet* AliAnalysisTaskSAJF::GetEmbJet(const Int_t i) const
 480 {
 481   if (fEmbJets)
 482     return dynamic_cast<AliEmcalJet*>(fEmbJets->At(i));
 483   else
 484     return 0;
 485 }
 486
 487 //________________________________________________________________________
 488 Int_t AliAnalysisTaskSAJF::GetNumberOfEmbJets() const
 489 {
 490   if (fEmbJets)
 491     return fEmbJets->GetEntriesFast();
 492   else
 493     return 0;
 494 }
 495
 496 //________________________________________________________________________
 497 AliVCluster* AliAnalysisTaskSAJF::GetTrgCluster(const Int_t i) const
 498 {
 499   if (fTrgClusters)
 500     return dynamic_cast<AliVCluster*>(fTrgClusters->At(i));
 501   else
 502     return 0;
 503 }
 504
 505 //________________________________________________________________________
 506 Int_t AliAnalysisTaskSAJF::GetNumberOfTrgClusters() const
 507 {
 508   if (fTrgClusters)
 509     return fTrgClusters->GetEntriesFast();
 510   else
 511     return 0;
 512 }
 513
 514 //________________________________________________________________________
 515 void AliAnalysisTaskSAJF::FillHistograms()
 516 {
 517   Float_t A = fJetRadius * fJetRadius * TMath::Pi();
 518
 519   Float_t cent = 100;
 520
 521   if (fCent)
 522     cent = fCent->GetCentralityPercentile("V0M");
 523
 524   fHistCentrality->Fill(cent);
 525
 526   Int_t maxJetIndex  = -1;
 527   Int_t max2JetIndex = -1;
 528
 529   DoJetLoop(maxJetIndex, max2JetIndex);
 530
 531   if (maxJetIndex < 0)
 532     return;
 533
 534   AliEmcalJet* jet = GetJet(maxJetIndex);
 535   if (!jet)
 536     return;
 537
 538   fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
 539   jet->SortConstituents();
 540
 541   AliEmcalJet* jet2 = 0;
 542   if (max2JetIndex >= 0)
 543     jet2 = GetJet(max2JetIndex);
 544
 545   if (jet2) {
 546     //fHistLeadingJetPt[fCentBin]->Fill(jet2->Pt());
 547     jet2->SortConstituents();
 548   }
 549
 550   Float_t rhoKt = DoKtJetLoop();
 551   fHistMedianPtKtJet[fCentBin]->Fill(rhoKt);
 552
 553   Float_t rhoTracks = DoTrackLoop(maxJetIndex, max2JetIndex);
 554   fHistBkgTracksMeanRho[fCentBin]->Fill(rhoTracks);
 555
 556   Float_t rhoClus = 0;
 557   if (fAnaType == kEMCAL || fAnaType == kEMCALFiducial)
 558     rhoClus = DoClusterLoop(maxJetIndex, max2JetIndex);
 559
 560   fHistBkgClusMeanRho[fCentBin]->Fill(rhoClus);
 561
 562   fHistRhoPartVSleadJetPt->Fill(A * (rhoClus + rhoTracks), jet->Pt());
 563
 564   fHistMedKtVSRhoPart->Fill(rhoKt, rhoClus + rhoTracks);
 565
 566   Int_t nRCs = GetArea() / A - 1;
 567
 568   for (Int_t i = 0; i < nRCs; i++) {
 569     Float_t RCpt = GetRigidConePt(0);
 570     Float_t RCptExLJ = GetRigidConePt(jet);
 571
 572     fHistDeltaPtRC[fCentBin]->Fill(RCpt - A * rhoKt);
 573     fHistDeltaPtRCExLJ[fCentBin]->Fill(RCptExLJ - A * rhoKt);
 574
 575     fHistRCPt[fCentBin]->Fill(RCpt / A);
 576     fHistRCPtExLJ[fCentBin]->Fill(RCptExLJ / A);
 577     fHistRCPtVSRhoPart->Fill(RCptExLJ / A, rhoClus + rhoTracks);
 578   }
 579
 580   Float_t maxEmbJetPt = 0;
 581   Float_t maxEmbPartPt = 0;
 582
 583   Bool_t embOK = DoEmbJetLoop(maxEmbJetPt, maxEmbPartPt);
 584
 585   if (embOK) {
 586     fHistEmbJets[fCentBin]->Fill(maxEmbJetPt);
 587     fHistEmbPart[fCentBin]->Fill(maxEmbPartPt);
 588     fHistDeltaPtMedKtEmb[fCentBin]->Fill(maxEmbJetPt - A * rhoKt - maxEmbPartPt);
 589   }
 590   else {
 591     if (maxEmbPartPt != 0)
 592       AliWarning("Embedded particle is not the leading particle of the leading jet!");
 593   }
 594 }
 595
 596 //________________________________________________________________________
 597 void AliAnalysisTaskSAJF::DoJetLoop(Int_t &maxJetIndex, Int_t &max2JetIndex)
 598 {
 599   Double_t vertex[3] = {0, 0, 0};
 600   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 601
 602   Int_t njets = GetNumberOfJets();
 603
 604   Float_t maxJetPt = 0;
 605   Float_t max2JetPt = 0;
 606   for (Int_t ij = 0; ij < njets; ij++) {
 607
 608     AliEmcalJet* jet = GetJet(ij);
 609
 610     if (!jet) {
 611       AliError(Form("Could not receive jet %d", ij));
 612       continue;
 613     }
 614
 615     if (jet->Pt() <= 0)
 616       continue;
 617
 618     if (!AcceptJet(jet))
 619       continue;
 620
 621     fHistJetPhiEta->Fill(jet->Eta(), jet->Phi());
 622     fHistJetsPt[fCentBin]->Fill(jet->Pt());
 623     fHistJetsNEF[fCentBin]->Fill(jet->NEF());
 624
 625     for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
 626       AliVTrack *track = jet->TrackAt(it, fTracks);
 627       if (track)
 628         fHistJetsZ[fCentBin]->Fill(track->Pt() / jet->Pt());
 629     }
 630
 631     for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
 632       AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
 633
 634       if (cluster) {
 635         TLorentzVector nPart;
 636         cluster->GetMomentum(nPart, vertex);
 637         fHistJetsZ[fCentBin]->Fill(nPart.Et() / jet->Pt());
 638       }
 639     }
 640
 641     if (maxJetPt < jet->Pt()) {
 642       max2JetPt = maxJetPt;
 643       max2JetIndex = maxJetIndex;
 644       maxJetPt = jet->Pt();
 645       maxJetIndex = ij;
 646     }
 647     else if (max2JetPt < jet->Pt()) {
 648       max2JetPt = jet->Pt();
 649       max2JetIndex = ij;
 650     }
 651   } //jet loop
 652 }
 653
 654 //________________________________________________________________________
 655 Float_t AliAnalysisTaskSAJF::DoKtJetLoop(Int_t nLJs)
 656 {
 657   Float_t ktJetsMedian = 0;
 658   Int_t nktjets =  GetNumberOfKtJets();
 659
 660   Int_t NoOfZeroJets = 0;
 661   if (nktjets > 0) {
 662
 663     TArrayF ktJets(nktjets);
 664     for (Int_t ij = 0; ij < nktjets; ij++) {
 665
 666       AliEmcalJet* jet = GetKtJet(ij);
 667
 668       if (!jet) {
 669         AliError(Form("Could not receive jet %d", ij));
 670         continue;
 671       }
 672
 673       if (jet->Pt() <= 0) {
 674         NoOfZeroJets++;
 675         continue;
 676       }
 677
 678       if (!AcceptJet(jet)) {
 679         NoOfZeroJets++;
 680         continue;
 681       }
 682
 683       Float_t rho = jet->Pt() / jet->Area();
 684       Int_t i = nktjets - 1;
 685       while (rho < ktJets[i] && i > 0)
 686         i--;
 687       memmove(ktJets.GetArray() + nktjets - ij - 1, ktJets.GetArray() + nktjets - ij, (ij + i - nktjets + 1) * sizeof(Float_t));
 688       ktJets[i] = rho;
 689     } //kt jet loop
 690
 691     nktjets -= NoOfZeroJets;
 692
 693     if (nktjets < 1) return 0;
 694
 695     memmove(ktJets.GetArray(), ktJets.GetArray() + NoOfZeroJets, nktjets * sizeof(Float_t));
 696
 697     nktjets -= nLJs;
 698
 699     if (nktjets < 1) return 0;
 700
 701     if (nktjets % 2)
 702       ktJetsMedian = ktJets[nktjets / 2];
 703     else
 704       ktJetsMedian = (ktJets[nktjets / 2] + ktJets[nktjets / 2 - 1]) / 2;
 705   }
 706
 707   return ktJetsMedian;
 708 }
 709
 710 //________________________________________________________________________
 711 Bool_t AliAnalysisTaskSAJF::DoEmbJetLoop(Float_t &maxJetPt, Float_t &maxPartPt)
 712 {
 713   Double_t vertex[3] = {0, 0, 0};
 714   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 715
 716   Int_t nembjets = GetNumberOfEmbJets();
 717
 718   Int_t maxJetIdx = -1;
 719   Int_t maxTrackIdx = -1;
 720   Int_t maxClusIdx = -1;
 721
 722   Float_t maxTrackPt = 0;
 723   Float_t maxClusEt = 0;
 724
 725   for (Int_t ij = 0; ij < nembjets; ij++) {
 726
 727     AliEmcalJet* jet = GetEmbJet(ij);
 728
 729     if (!jet) {
 730       AliError(Form("Could not receive jet %d", ij));
 731       continue;
 732     }
 733
 734     if (jet->Pt() <= 0)
 735         continue;
 736
 737     if (!AcceptJet(jet))
 738       continue;
 739
 740     if (jet->Pt() > maxJetPt) {
 741       maxJetPt = jet->Pt();
 742       maxJetIdx = ij;
 743     }
 744   }
 745
 746   if (maxJetPt <= 0)
 747     return kFALSE;
 748
 749   AliEmcalJet* jet = GetEmbJet(maxJetIdx);
 750
 751   for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
 752     AliVTrack *track = jet->TrackAt(it, fTracks);
 753
 754     if (!track) continue;
 755
 756     if (track->Pt() > maxTrackPt) {
 757       maxTrackPt = track->Pt();
 758       maxTrackIdx = it;
 759     }
 760   }
 761
 762   for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
 763     AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
 764
 765     if (!cluster) continue;
 766
 767     TLorentzVector nPart;
 768     cluster->GetMomentum(nPart, vertex);
 769
 770     if (nPart.Et() > maxClusEt) {
 771       maxClusEt = nPart.Et();
 772       maxClusIdx = ic;
 773     }
 774   }
 775
 776   if (maxClusEt > maxTrackPt) {
 777     AliVCluster *cluster = jet->ClusterAt(maxClusIdx, fCaloClusters);
 778     maxPartPt = maxClusEt;
 779     return (Bool_t)(cluster->Chi2() == 100);
 780   }
 781   else {
 782     AliVTrack *track = jet->TrackAt(maxTrackIdx, fTracks);
 783     maxPartPt = maxTrackPt;
 784     return (Bool_t)(track->GetLabel() == 100);
 785   }
 786 }
 787
 788 //________________________________________________________________________
 789 Float_t AliAnalysisTaskSAJF::DoTrackLoop(Int_t maxJetIndex, Int_t max2JetIndex)
 790 {
 791   AliEmcalJet* jet = 0;
 792   if (max2JetIndex >= 0)
 793     jet = GetJet(maxJetIndex);
 794
 795   AliEmcalJet* jet2 = 0;
 796   if (max2JetIndex >= 0)
 797     jet2 = GetJet(max2JetIndex);
 798
 799   Float_t rhoTracks = 0;
 800   Int_t ntracks = GetNumberOfTracks();
 801
 802   for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 803     AliVTrack* track = GetTrack(iTracks);
 804     if(!track) {
 805       AliError(Form("Could not retrieve track %d",iTracks));
 806       continue;
 807     }
 808
 809     if (!AcceptTrack(track)) continue;
 810
 811     if (jet && IsJetTrack(jet, iTracks)) {
 812       fHistTracksPtLJ[fCentBin]->Fill(track->Pt());
 813     }
 814     else if (!jet2 || !IsJetTrack(jet2, iTracks)) {
 815       fHistTracksPtBkg[fCentBin]->Fill(track->Pt());
 816       rhoTracks += track->Pt();
 817
 818       if (jet) {
 819         Float_t dphijet = jet->Phi() - track->Phi();
 820         if (dphijet < -1.6) dphijet += TMath::Pi() * 2;
 821         if (dphijet > 4.8) dphijet -= TMath::Pi() * 2;
 822         fHistBkgLJetPhiCorr[fCentBin]->Fill(dphijet);
 823       }
 824
 825       for(Int_t it2 = iTracks+1; it2 < ntracks; it2++) {
 826         AliVTrack* track2 = GetTrack(it2);
 827         if(!track2) {
 828           AliError(Form("Could not retrieve track %d", it2));
 829           continue;
 830         }
 831
 832         if (!AcceptTrack(track2)) continue;
 833
 834         if (jet && IsJetTrack(jet, it2)) continue;
 835
 836         if (jet2 && IsJetTrack(jet2, it2)) continue;
 837
 838         Float_t dphi = track->Phi() - track2->Phi();
 839         if (dphi < -1.6) dphi += TMath::Pi() * 2;
 840         if (dphi > 4.8) dphi -= TMath::Pi() * 2;
 841         fHistBkgTracksPhiCorr[fCentBin]->Fill(dphi);
 842       } // second track loop
 843     }
 844   }
 845   rhoTracks /= GetArea();
 846   return rhoTracks;
 847 }
 848
 849 //________________________________________________________________________
 850 Float_t AliAnalysisTaskSAJF::DoClusterLoop(Int_t maxJetIndex, Int_t max2JetIndex)
 851 {
 852   Double_t vertex[3] = {0, 0, 0};
 853   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 854
 855   AliEmcalJet* jet = 0;
 856   if (max2JetIndex >= 0)
 857     jet = GetJet(maxJetIndex);
 858
 859   AliEmcalJet* jet2 = 0;
 860   if (max2JetIndex >= 0)
 861     jet2 = GetJet(max2JetIndex);
 862
 863   Float_t rhoClus = 0;
 864   Int_t nclusters =  GetNumberOfCaloClusters();
 865   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 866     AliVCluster* cluster = GetCaloCluster(iClusters);
 867     if (!cluster) {
 868       AliError(Form("Could not receive cluster %d", iClusters));
 869       continue;
 870     }
 871
 872     if (!(cluster->IsEMCAL())) continue;
 873
 874     if (!AcceptCluster(cluster)) continue;
 875
 876     TLorentzVector nPart;
 877     cluster->GetMomentum(nPart, vertex);
 878
 879     if (jet && IsJetCluster(jet, iClusters)) {
 880       fHistClusEtLJ[fCentBin]->Fill(nPart.Et());
 881     }
 882     else if (!jet2 || !IsJetCluster(jet2, iClusters)) {
 883       fHistClusEtBkg[fCentBin]->Fill(nPart.Et());
 884       rhoClus += nPart.Et();
 885
 886       Float_t pos1[3];
 887       cluster->GetPosition(pos1);
 888       TVector3 clusVec1(pos1);
 889
 890       if (jet) {
 891         Float_t dphijet = jet->Phi() - clusVec1.Phi();
 892         if (dphijet < -1.6) dphijet += TMath::Pi() * 2;
 893         if (dphijet > 4.8) dphijet -= TMath::Pi() * 2;
 894         fHistBkgLJetPhiCorr[fCentBin]->Fill(dphijet);
 895       }
 896
 897       for(Int_t ic2 = iClusters+1; ic2 < nclusters; ic2++) {
 898         AliVCluster* cluster2 = GetCaloCluster(ic2);
 899         if (!cluster2) {
 900           AliError(Form("Could not receive cluster %d", ic2));
 901           continue;
 902         }
 903
 904         if (!(cluster2->IsEMCAL())) continue;
 905
 906         if (!AcceptCluster(cluster)) continue;
 907
 908         if (jet && IsJetCluster(jet, ic2)) continue;
 909
 910         if (jet2 && IsJetCluster(jet2, ic2)) continue;
 911
 912         Float_t pos2[3];
 913         cluster2->GetPosition(pos2);
 914         TVector3 clusVec2(pos2);
 915
 916         Float_t dphi = clusVec1.Phi() - clusVec2.Phi();
 917         if (dphi < -1.6) dphi += TMath::Pi() * 2;
 918         if (dphi > 4.8) dphi -= TMath::Pi() * 2;
 919         fHistBkgClusPhiCorr[fCentBin]->Fill(dphi);
 920       }
 921     }
 922   } //cluster loop
 923   rhoClus /= GetArea();
 924   return rhoClus;
 925 }
 926
 927 //________________________________________________________________________
 928 void AliAnalysisTaskSAJF::Init()
 929 {
 930   if (fAnaType == kFullAcceptance) {
 931     fMinEta = -1;
 932     fMaxEta = 1;
 933     fMinPhi = 0;
 934     fMaxPhi = 2 * TMath::Pi();
 935   }
 936   else if (fAnaType == kEMCAL || fAnaType == kEMCALFiducial) {
 937     fMinEta = -0.7;
 938     fMaxEta = 0.7;
 939     fMinPhi = 80 * TMath::DegToRad();
 940     fMaxPhi = 180 * TMath::DegToRad();
 941   }
 942   else {
 943     AliWarning("Analysis type not recognized! Assuming kFullAcceptance...");
 944     fMinEta = -1;
 945     fMaxEta = 1;
 946     fMinPhi = 0;
 947     fMaxPhi = 2 * TMath::Pi();
 948   }
 949 }
 950
 951 //________________________________________________________________________
 952 Float_t AliAnalysisTaskSAJF::GetRigidConePt(AliEmcalJet *jet, Float_t minD)
 953 {
 954   static TRandom3 random;
 955
 956   Double_t vertex[3] = {0, 0, 0};
 957   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 958
 959   Float_t eta = 0;
 960   Float_t phi = 0;
 961
 962   Float_t LJeta = 999;
 963   Float_t LJphi = 999;
 964
 965   if (jet) {
 966     LJeta = jet->Eta();
 967     LJphi = jet->Phi();
 968   }
 969
 970   Float_t dLJ = 0;
 971   do {
 972     eta = random.Rndm() * (fMaxEta - fMinEta) + fMinEta;
 973     phi = random.Rndm() * (fMaxPhi - fMinPhi) + fMinPhi;
 974     dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
 975   } while (dLJ < minD && !AcceptJet(eta, phi));
 976
 977   Float_t pt = 0;
 978
 979   Int_t nclusters =  GetNumberOfCaloClusters();
 980   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 981     AliVCluster* cluster = GetCaloCluster(iClusters);
 982     if (!cluster) {
 983       AliError(Form("Could not receive cluster %d", iClusters));
 984       continue;
 985     }
 986
 987     if (!(cluster->IsEMCAL())) continue;
 988
 989     if (!AcceptCluster(cluster)) continue;
 990
 991     TLorentzVector nPart;
 992     cluster->GetMomentum(nPart, vertex);
 993
 994     Float_t pos[3];
 995     cluster->GetPosition(pos);
 996     TVector3 clusVec(pos);
 997
 998     Float_t d = TMath::Sqrt((clusVec.Eta() - eta) * (clusVec.Eta() - eta) + (clusVec.Phi() - phi) * (clusVec.Phi() - phi));
 999
1000     if (d <= fJetRadius)
1001       pt += nPart.Pt();
1002   }
1003
1004   Int_t ntracks = GetNumberOfTracks();
1005   for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1006     AliVTrack* track = GetTrack(iTracks);
1007     if(!track) {
1008       AliError(Form("Could not retrieve track %d",iTracks));
1009       continue;
1010     }
1011
1012     if (!AcceptTrack(track)) continue;
1013
1014     Float_t tracketa = track->Eta();
1015     Float_t trackphi = track->Phi();
1016
1017     if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
1018       trackphi += 2 * TMath::Pi();
1019     if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
1020       trackphi -= 2 * TMath::Pi();
1021
1022     Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
1023
1024     if (d <= fJetRadius)
1025       pt += track->Pt();
1026   }
1027
1028   return pt;
1029 }
1030
1031 //________________________________________________________________________
1032 Float_t AliAnalysisTaskSAJF::GetArea() const
1033 {
1034   return (fMaxEta - fMinEta) * (fMaxPhi - fMinPhi);
1035 }
1036
1037 //________________________________________________________________________
1038 Bool_t AliAnalysisTaskSAJF::IsJetTrack(AliEmcalJet* jet, Int_t itrack, Bool_t sorted) const
1039 {
1040   for (Int_t i = 0; i < jet->GetNumberOfTracks(); i++) {
1041     Int_t ijettrack = jet->TrackAt(i);
1042     if (sorted && ijettrack > itrack)
1043       return kFALSE;
1044     if (ijettrack == itrack)
1045       return kTRUE;
1046   }
1047   return kFALSE;
1048 }
1049
1050 //________________________________________________________________________
1051 Bool_t AliAnalysisTaskSAJF::IsJetCluster(AliEmcalJet* jet, Int_t iclus, Bool_t sorted) const
1052 {
1053   for (Int_t i = 0; i < jet->GetNumberOfClusters(); i++) {
1054     Int_t ijetclus = jet->ClusterAt(i);
1055     if (sorted && ijetclus > iclus)
1056       return kFALSE;
1057     if (ijetclus == iclus)
1058       return kTRUE;
1059   }
1060   return kFALSE;
1061 }
1062
1063 //________________________________________________________________________
1064 Bool_t AliAnalysisTaskSAJF::AcceptJet(Float_t eta, Float_t phi) const
1065 {
1066   if (fAnaType == kFullAcceptance) {
1067     return kTRUE;
1068   }
1069   else if (fAnaType == kEMCAL) {
1070     return (Bool_t)(eta > fMinEta && eta < fMaxEta && phi > fMinPhi && phi < fMaxPhi);
1071   }
1072   else if (fAnaType == kEMCALFiducial) {
1073     return (Bool_t)(eta > fMinEta + fJetRadius && eta < fMaxEta - fJetRadius && phi > fMinPhi + fJetRadius && phi < fMaxPhi - fJetRadius);
1074   }
1075   else {
1076     AliWarning("Analysis type not recognized! Assuming kFullAcceptance...");
1077     return kTRUE;
1078   }
1079 }
1080
1081 //________________________________________________________________________
1082 Bool_t AliAnalysisTaskSAJF::AcceptJet(AliEmcalJet* jet) const
1083 {
1084   return AcceptJet(jet->Eta(), jet->Phi());
1085 }
1086
1087 //________________________________________________________________________
1088 Bool_t AliAnalysisTaskSAJF::AcceptCluster(AliVCluster* clus, Bool_t acceptMC) const
1089 {
1090   if (!acceptMC && clus->Chi2() == 100)
1091     return kFALSE;
1092
1093   return kTRUE;
1094 }
1095
1096
1097 //________________________________________________________________________
1098 Bool_t AliAnalysisTaskSAJF::AcceptTrack(AliVTrack* track, Bool_t acceptMC) const
1099 {
1100   if (!acceptMC && track->GetLabel() == 100)
1101     return kFALSE;
1102
1103   if (fAnaType == kFullAcceptance) {
1104     return kTRUE;
1105   }
1106   else if (fAnaType == kEMCAL || fAnaType == kEMCALFiducial) {
1107     return (Bool_t)(track->Eta() > fMinEta && track->Eta() < fMaxEta && track->Phi() > fMinPhi && track->Phi() < fMaxPhi);
1108   }
1109   else {
1110     AliWarning("Analysis type not recognized! Assuming kFullAcceptance...");
1111     return kTRUE;
1112   }
1113 }
1114
1115 //________________________________________________________________________
1116 void AliAnalysisTaskSAJF::UserExec(Option_t *)
1117 {
1118   Init();
1119
1120   RetrieveEventObjects();
1121
1122   FillHistograms();
1123
1124   // information for this iteration of the UserExec in the container
1125   PostData(1, fOutput);
1126 }
1127
1128 //________________________________________________________________________
1129 void AliAnalysisTaskSAJF::Terminate(Option_t *)
1130 {
1131   // Called once at the end of the analysis.
1132 }