PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskQualityAssurancePA.cxx

   1 #ifndef ALIANALYSISTASKSE_H
   2 #include <Riostream.h>
   3 #include <TROOT.h>
   4 #include <TString.h>
   5 #include <TFile.h>
   6 #include <TCint.h>
   7 #include <TChain.h>
   8 #include <TTree.h>
   9 #include <TKey.h>
  10 #include <TProfile.h>
  11 #include <TH1F.h>
  12 #include <TH2F.h>
  13 #include <TCanvas.h>
  14 #include <TList.h>
  15 #include <TClonesArray.h>
  16 #include <TObject.h>
  17 #include <TMath.h>
  18 #include <TSystem.h>
  19 #include <TInterpreter.h>
  20 #include <TH1.h>
  21 #include "AliAnalysisTask.h"
  22 #include "AliCentrality.h"
  23 #include "AliStack.h"
  24 #include "AliESDEvent.h"
  25 #include "AliESDInputHandler.h"
  26 #include "AliAODEvent.h"
  27 #include "AliAODHandler.h"
  28 #include "AliAnalysisManager.h"
  29 #include "AliAnalysisTaskSE.h"
  30 #endif
  31
  32 #include <TRandom3.h>
  33 #include "AliGenPythiaEventHeader.h"
  34 #include "AliMCEvent.h"
  35 #include "AliLog.h"
  36 #include <AliEmcalJet.h>
  37 #include <AliRhoParameter.h>
  38 #include "AliVEventHandler.h"
  39 #include "AliVParticle.h"
  40 #include "AliAnalysisUtils.h"
  41
  42 #include "AliAnalysisTaskQualityAssurancePA.h"
  43
  44
  45 //TODO: FillHistogram can be done better with virtual TH1(?)
  46 ClassImp(AliAnalysisTaskQualityAssurancePA)
  47
  48 // ######################################################################################## DEFINE HISTOGRAMS
  49 void AliAnalysisTaskQualityAssurancePA::Init()
  50 {
  51   #ifdef DEBUGMODE
  52     AliInfo("Creating histograms.");
  53   #endif
  54   *fRunNumbers = *fRunNumbers + " ALL";
  55
  56   TObjArray* runNumberArr = fRunNumbers->Tokenize(" ");
  57
  58   for(Int_t i=0; i<runNumberArr->GetEntries();i++)
  59   {
  60     const char* tmpRunNum = (static_cast<TObjString*>(runNumberArr->At(i)))->String().Data();
  61
  62     // NOTE: Track & Cluster & QA histograms
  63     if (fAnalyzeQA)
  64     {
  65       AddHistogram1D<TH1D>(tmpRunNum, "hNumberEvents", "Number of events (0 = before, 1 = after vertex cuts)", "", 2, 0, 2, "#Delta z(cm)","N^{Events}/cut");
  66       AddHistogram1D<TH1D>(tmpRunNum, "hVertexX", "X distribution of the vertex", "", 500, -1., 1., "#Delta x(cm)","dN^{Events}/dx");
  67       AddHistogram1D<TH1D>(tmpRunNum, "hVertexY", "Y distribution of the vertex", "", 500, -1., 1., "#Delta y(cm)","dN^{Events}/dy");
  68       AddHistogram2D<TH2D>(tmpRunNum, "hVertexXY", "XY distribution of the vertex", "COLZ", 250, -1., 1., 250, -1., 1.,"#Delta x(cm)", "#Delta y(cm)","dN^{Events}/dxdy");
  69       AddHistogram1D<TH1D>(tmpRunNum, "hVertexZ", "Z distribution of the vertex", "", 100, -10., 10., "#Delta z(cm)","dN^{Events}/dz");
  70       AddHistogram1D<TH1D>(tmpRunNum, "hVertexR", "R distribution of the vertex", "", 100, 0., 1., "#Delta r(cm)","dN^{Events}/dr");
  71       AddHistogram1D<TH1D>(tmpRunNum, "hCentralityV0M", "Centrality distribution V0M", "", 100, 0., 100., "Centrality","dN^{Events}");
  72       AddHistogram1D<TH1D>(tmpRunNum, "hCentralityV0A", "Centrality distribution V0A", "", 100, 0., 100., "Centrality","dN^{Events}");
  73       AddHistogram1D<TH1D>(tmpRunNum, "hCentralityV0C", "Centrality distribution V0C", "", 100, 0., 100., "Centrality","dN^{Events}");
  74
  75       AddHistogram2D<TH2D>(tmpRunNum, "hTrackCountAcc", "Number of tracks in acceptance vs. centrality", "LEGO2", 750, 0., 750., 100, 0, 100, "N tracks","Centrality", "dN^{Events}/dN^{Tracks}");
  76       AddHistogram1D<TH1D>(tmpRunNum, "hTrackPt", "Tracks p_{T} distribution", "", 1000, 0., 250., "p_{T} (GeV/c)","dN^{Tracks}/dp_{T}");
  77       AddHistogram1D<TH1D>(tmpRunNum, "hTrackPtNegEta", "Tracks p_{T} distribution (negative #eta)", "", 1000, 0., 250., "p_{T} (GeV/c)","dN^{Tracks}/dp_{T}");
  78       AddHistogram1D<TH1D>(tmpRunNum, "hTrackPtPosEta", "Tracks p_{T} distribution (positive #eta)", "", 1000, 0., 250., "p_{T} (GeV/c)","dN^{Tracks}/dp_{T}");
  79       AddHistogram1D<TH1D>(tmpRunNum, "hTrackCharge", "Charge", "", 11, -5, 5, "Charge (e)","dN^{Tracks}/dq");
  80       AddHistogram1D<TH1D>(tmpRunNum, "hTrackPhi", "Track #phi distribution", "", 360, 0, TMath::TwoPi(), "#phi","dN^{Tracks}/d#phi");
  81       AddHistogram2D<TH2D>(tmpRunNum, "hTrackPhiEta", "Track angular distribution", "LEGO2", 100, 0., 2*TMath::Pi(),100, -2.5, 2.5, "#phi","#eta","dN^{Tracks}/(d#phi d#eta)");
  82
  83       AddHistogram2D<TH2D>(tmpRunNum, "hTrackPhiPtCut", "Track #phi distribution for different pT cuts", "LEGO2", 360, 0, TMath::TwoPi(), 20, 0, 20, "#phi", "p_{T} lower cut", "dN^{Tracks}/d#phi dp_{T}");
  84       AddHistogram2D<TH2D>(tmpRunNum, "hTrackPhiLabel", "Track #phi distribution in different labels", "LEGO2", 360, 0, TMath::TwoPi(), 3, 0, 3, "#phi", "Label", "dN^{Tracks}/d#phi");
  85       AddHistogram1D<TH1D>(tmpRunNum, "hTrackEta", "Track #eta distribution", "", 180, -fTrackEtaWindow, +fTrackEtaWindow, "#eta","dN^{Tracks}/d#eta");
  86     }
  87
  88     // NOTE: Pythia histograms
  89     if (fAnalyzePythia)
  90     {
  91       AddHistogram1D<TH1D>(tmpRunNum, "hPythiaPtHard", "Pythia p_{T} hard distribution", "", 1000, 0, 400, "p_{T} hard","dN^{Events}/dp_{T,hard}");
  92       AddHistogram1D<TProfile>(tmpRunNum, "hPythiaXSection", "Pythia cross section distribution", "", fNumPtHardBins, 0, fNumPtHardBins, "p_{T} hard bin","dN^{Events}/dp_{T,hard}");
  93       AddHistogram1D<TH1D>(tmpRunNum, "hPythiaNTrials", "Pythia trials (no correction for manual cuts)", "", fNumPtHardBins, 0, fNumPtHardBins, "p_{T} hard bin", "Trials");
  94     }
  95
  96     // NOTE: Jet histograms
  97     if (fAnalyzeJets)
  98     {
  99       // ######## Jet spectra
 100       AddHistogram1D<TH1D>(tmpRunNum, "hJetPt", "Jets p_{T} distribution", "", 800, 0., 200., "p_{T} (GeV/c)","dN^{Jets}/dp_{T}");
 101       AddHistogram1D<TH1D>(tmpRunNum, "hJetArea", "Jets area distribution", "", 200, 0., 2., "Area","dN^{Jets}/dA");
 102       AddHistogram2D<TH2D>(tmpRunNum, "hJetAreaVsPt", "Jets area vs. p_{T} distribution", "COLZ", 100, 0., 2., 200, 0., 200., "Area", "p_{T}", "dN^{Jets}/dA dp_{T}");
 103
 104       AddHistogram2D<TH2D>(tmpRunNum, "hJetPhiEta", "Jets angular distribution", "LEGO2", 360, 0., 2*TMath::Pi(),100, -0.6, 0.6, "#phi","#eta","dN^{Jets}/(d#phi d#eta)");
 105       AddHistogram2D<TH2D>(tmpRunNum, "hJetPtVsConstituentCount", "Jets number of constituents vs. jet p_{T}", "COLZ", 400, 0., 200., 100, 0., 100., "p_{T}","N^{Tracks}","dN^{Jets}/(dp_{T} dN^{tracks})");
 106       AddHistogram1D<TH1D>(tmpRunNum, "hJetCountAll", "Number of Jets", "", 100, 0., 100., "N jets","dN^{Events}/dN^{Jets}");
 107       AddHistogram1D<TH1D>(tmpRunNum, "hJetCountAccepted", "Number of accepted Jets", "", 100, 0., 100., "N jets","dN^{Events}/dN^{Jets}");
 108
 109     }
 110   }
 111   // register Histograms
 112   for (Int_t i = 0; i < fHistCount; i++)
 113   {
 114     fOutputList->Add(fHistList->At(i));
 115   }
 116
 117   PostData(1,fOutputList); // important for merging
 118
 119 }
 120
 121 //________________________________________________________________________
 122 AliAnalysisTaskQualityAssurancePA::AliAnalysisTaskQualityAssurancePA(const char *name, const char* trackArrayName, const char* clusterArrayName, const char* jetArrayName) : AliAnalysisTaskSE(name), fOutputList(0), fAnalyzeQA(1), fAnalyzeJets(1), fAnalyzePythia(0), fHasTracks(0), fHasClusters(0), fHasJets(0), fIsMC(0), fJetArray(0), fTrackArray(0), fClusterArray(0), fJetArrayName(0), fTrackArrayName(0), fClusterArrayName(0), fRunNumbers(0), fNumPtHardBins(11), fSignalJetRadius(0.4), fNumberExcludedJets(2), fSignalJetEtaWindow(0.5), fTrackEtaWindow(0.9), fClusterEtaWindow(0.7), fVertexWindow(10.0), fVertexMaxR(1.0), fMinTrackPt(0.150), fMinClusterPt(0.300), fMinJetPt(1.0), fMinJetArea(0.4), fFirstLeadingJet(0), fSecondLeadingJet(0), fNumberSignalJets(0), fCrossSection(0.0), fTrials(0.0), fRandom(0), fHelperClass(0), fInitialized(0), fTaskInstanceCounter(0), fHistList(0), fHistCount(0)
 123 {
 124   #ifdef DEBUGMODE
 125     AliInfo("Calling constructor.");
 126   #endif
 127
 128   // Constructor
 129   // Define input and output slots here (never in the dummy constructor)
 130   // Input slot #0 works with a TChain - it is connected to the default input container
 131   // Output slot #1 writes into a TH1 container
 132   // Constructor
 133
 134   // Every instance of this task gets his own number
 135   static Int_t instance = 0;
 136   fTaskInstanceCounter = instance;
 137   instance++;
 138
 139   fTrackArrayName = new TString(trackArrayName);
 140   fClusterArrayName = new TString(clusterArrayName);
 141   fRunNumbers = new TString("");
 142   if (strcmp(fTrackArrayName->Data(),"") == 0)
 143     fAnalyzeQA = kFALSE;
 144   else
 145   {
 146     fAnalyzeQA = kTRUE;
 147     if (fTrackArrayName->Contains("MCParticles")) //TODO: Hardcoded for now
 148       fIsMC = kTRUE;
 149   }
 150
 151   fJetArrayName = new TString(jetArrayName);
 152   if (strcmp(fJetArrayName->Data(),"") == 0)
 153     fAnalyzeJets = kFALSE;
 154   else
 155     fAnalyzeJets = kTRUE;
 156
 157   DefineOutput(1, TList::Class());
 158
 159   fHistList = new TList();
 160
 161   #ifdef DEBUGMODE
 162     AliInfo("Constructor done.");
 163   #endif
 164
 165 }
 166
 167 //________________________________________________________________________
 168 inline Double_t AliAnalysisTaskQualityAssurancePA::GetPtHard()
 169 {
 170   Double_t tmpPtHard = -1.0;
 171
 172   if (!MCEvent())
 173     AliError("MCEvent not accessible although demanded!");
 174   else
 175   {
 176     AliGenPythiaEventHeader* pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
 177     if (!pythiaHeader)
 178       AliError("Pythia Header not accessible!");
 179     else
 180       tmpPtHard = pythiaHeader->GetPtHard();
 181   }
 182   return tmpPtHard;
 183 }
 184
 185 //________________________________________________________________________
 186 inline Int_t AliAnalysisTaskQualityAssurancePA::GetPtHardBin()
 187 {
 188   // ########## PT HARD BIN EDGES
 189   const Int_t localkPtHardLowerEdges[] = { 0, 5,11,21,36,57, 84,117,152,191,234};
 190   const Int_t localkPtHardHigherEdges[] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 191
 192   Int_t tmpPtHardBin = 0;
 193   Double_t tmpPtHard = GetPtHard();
 194
 195   for (tmpPtHardBin = 0; tmpPtHardBin <= fNumPtHardBins; tmpPtHardBin++)
 196     if (tmpPtHard >= localkPtHardLowerEdges[tmpPtHardBin] && tmpPtHard < localkPtHardHigherEdges[tmpPtHardBin])
 197       break;
 198
 199   return tmpPtHardBin;
 200 }
 201
 202 //________________________________________________________________________
 203 inline Bool_t AliAnalysisTaskQualityAssurancePA::IsTrackInAcceptance(AliVParticle* track)
 204 {
 205   if (track != 0)
 206     if (TMath::Abs(track->Eta()) <= fTrackEtaWindow)
 207       if (track->Pt() >= fMinTrackPt)
 208         return kTRUE;
 209
 210   return kFALSE;
 211 }
 212
 213 //________________________________________________________________________
 214 inline Bool_t AliAnalysisTaskQualityAssurancePA::IsClusterInAcceptance(AliVCluster* cluster)
 215 {
 216   if (cluster != 0)
 217  //   if (TMath::Abs(cluster->Eta()) <= fClusterEtaWindow)
 218 //      if (cluster->Phi() <= 187.0/360.0 * TMath::TwoPi());
 219 //        if (cluster->Phi() >= 80.0/360.0 * TMath::TwoPi());
 220           if (cluster->E() >= fMinClusterPt)
 221             return kTRUE;
 222
 223   return kFALSE;
 224 }
 225
 226 //________________________________________________________________________
 227 inline Bool_t AliAnalysisTaskQualityAssurancePA::IsSignalJetInAcceptance(AliEmcalJet *jet)
 228 {
 229   if (jet != 0)
 230     if (TMath::Abs(jet->Eta()) <= fSignalJetEtaWindow)
 231       if (jet->Pt() >= fMinJetPt)
 232         if (jet->Area() >= fMinJetArea)
 233           return kTRUE;
 234
 235   return kFALSE;
 236 }
 237
 238 //________________________________________________________________________
 239 void AliAnalysisTaskQualityAssurancePA::ExecOnce()
 240 {
 241   #ifdef DEBUGMODE
 242     AliInfo("Starting ExecOnce.");
 243   #endif
 244   fInitialized = kTRUE;
 245
 246   // Check for track array
 247   if (strcmp(fTrackArrayName->Data(), "") != 0)
 248   {
 249     fTrackArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackArrayName->Data()));
 250     fHasTracks = kTRUE;
 251     if (!fTrackArray)
 252     {
 253       AliWarning(Form("%s: Could not retrieve tracks %s! This is OK, if tracks are not demanded.", GetName(), fTrackArrayName->Data()));
 254       fHasTracks = kFALSE;
 255     }
 256     else
 257     {
 258       TClass *cl = fTrackArray->GetClass();
 259       if (!cl->GetBaseClass("AliVParticle"))
 260       {
 261         AliError(Form("%s: Collection %s does not contain AliVParticle objects!", GetName(), fTrackArrayName->Data()));
 262         fTrackArray = 0;
 263         fHasTracks = kFALSE;
 264       }
 265     }
 266   }
 267   // Check for cluster array
 268   if (strcmp(fClusterArrayName->Data(), "") != 0)
 269   {
 270     fClusterArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fClusterArrayName->Data()));
 271     fHasClusters = kTRUE;
 272     if (!fClusterArray)
 273     {
 274       AliWarning(Form("%s: Could not retrieve clusters %s! This is OK, if clusters are not demanded.", GetName(), fClusterArrayName->Data()));
 275       fHasClusters = kFALSE;
 276     }
 277     else
 278     {
 279       TClass *cl = fClusterArray->GetClass();
 280       if (!cl->GetBaseClass("AliVCluster"))
 281       {
 282         AliError(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fClusterArrayName->Data()));
 283         fClusterArray = 0;
 284         fHasClusters = kFALSE;
 285       }
 286     }
 287   }
 288
 289   // Check for jet array
 290   if (strcmp(fJetArrayName->Data(), "") != 0)
 291   {
 292     fJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrayName->Data()));
 293     fHasJets = kTRUE;
 294
 295     if (!fJetArray)
 296     {
 297       AliWarning(Form("%s: Could not retrieve jets %s! This is OK, if jets are not demanded.", GetName(), fJetArrayName->Data()));
 298       fHasJets = kFALSE;
 299     }
 300     else
 301     {
 302       if (!fJetArray->GetClass()->GetBaseClass("AliEmcalJet"))
 303       {
 304         AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJetArrayName->Data()));
 305         fJetArray = 0;
 306         fHasJets = kFALSE;
 307       }
 308     }
 309   }
 310
 311   // Look, if initialization is OK
 312   if (!fHasTracks && fAnalyzeQA)
 313   {
 314     AliError(Form("%s: Tracks NOT successfully casted although demanded! Deactivating QA",GetName()));
 315     fAnalyzeQA = kFALSE;
 316   }
 317   if (!fHasJets && fAnalyzeJets)
 318   {
 319     AliError(Form("%s: Jets NOT successfully casted although demanded!  Deactivating jet",GetName()));
 320     fAnalyzeJets = kFALSE;
 321   }
 322
 323   // Initialize helper class (for vertex selection)
 324   fHelperClass = new AliAnalysisUtils();
 325
 326   // Histogram init
 327   Init();
 328
 329   #ifdef DEBUGMODE
 330     AliInfo("ExecOnce done.");
 331   #endif
 332
 333 }
 334
 335 //________________________________________________________________________
 336 void AliAnalysisTaskQualityAssurancePA::GetSignalJets()
 337 {
 338   // Reset vars
 339   fFirstLeadingJet = NULL;
 340   fSecondLeadingJet = NULL;
 341   fNumberSignalJets = 0;
 342
 343   Float_t maxJetPts[] = {0, 0};
 344   Int_t jetIDArray[]   = {-1, -1};
 345   Int_t jetCount = fJetArray->GetEntries();
 346
 347   // Go through all jets and save signal jets and the two leading ones
 348   for (Int_t i = 0; i < jetCount; i++)
 349   {
 350     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJetArray->At(i));
 351     if (!jet)
 352     {
 353       AliError(Form("%s: Could not receive jet %d", GetName(), i));
 354       continue;
 355     }
 356
 357     if (!IsSignalJetInAcceptance(jet)) continue;
 358
 359     if (jet->Pt() > maxJetPts[0])
 360     {
 361       maxJetPts[1] = maxJetPts[0];
 362       jetIDArray[1] = jetIDArray[0];
 363       maxJetPts[0] = jet->Pt();
 364       jetIDArray[0] = i;
 365     }
 366     else if (jet->Pt() > maxJetPts[1])
 367     {
 368       maxJetPts[1] = jet->Pt();
 369       jetIDArray[1] = i;
 370     }
 371     fSignalJets[fNumberSignalJets] = jet;
 372     fNumberSignalJets++;
 373   }
 374
 375   if (fNumberSignalJets > 0)
 376     fFirstLeadingJet  = static_cast<AliEmcalJet*>(fJetArray->At(jetIDArray[0]));
 377   if (fNumberSignalJets > 1)
 378     fSecondLeadingJet = static_cast<AliEmcalJet*>(fJetArray->At(jetIDArray[1]));
 379
 380 }
 381
 382 //________________________________________________________________________
 383 Int_t AliAnalysisTaskQualityAssurancePA::GetLeadingJets(TClonesArray* jetArray, Int_t* jetIDArray)
 384 {
 385 // Writes first two leading jets into already registered array jetIDArray
 386
 387   if (!jetArray)
 388   {
 389     AliError("Could not get the jet array to get leading jets from it!");
 390     return 0;
 391   }
 392
 393   Float_t maxJetPts[] = {0, 0};
 394   jetIDArray[0] = -1;
 395   jetIDArray[1] = -1;
 396
 397   Int_t jetCount = jetArray->GetEntries();
 398   Int_t jetCountAccepted = 0;
 399
 400   for (Int_t i = 0; i < jetCount; i++)
 401   {
 402     AliEmcalJet* jet = static_cast<AliEmcalJet*>(jetArray->At(i));
 403     if (!jet)
 404     {
 405       AliError(Form("%s: Could not receive jet %d", GetName(), i));
 406       continue;
 407     }
 408
 409     if (!IsSignalJetInAcceptance(jet)) continue;
 410
 411     if (jet->Pt() > maxJetPts[0])
 412     {
 413       maxJetPts[1] = maxJetPts[0];
 414       jetIDArray[1] = jetIDArray[0];
 415       maxJetPts[0] = jet->Pt();
 416       jetIDArray[0] = i;
 417     }
 418     else if (jet->Pt() > maxJetPts[1])
 419     {
 420       maxJetPts[1] = jet->Pt();
 421       jetIDArray[1] = i;
 422     }
 423     jetCountAccepted++;
 424   }
 425   return jetCountAccepted;
 426 }
 427
 428
 429 //________________________________________________________________________
 430 void AliAnalysisTaskQualityAssurancePA::Calculate(AliVEvent* event)
 431 {
 432   #ifdef DEBUGMODE
 433     AliInfo("Starting Calculate().");
 434   #endif
 435   ////////////////////// NOTE: initialization & casting
 436
 437   if (!event) {
 438     AliError("??? Event pointer == 0 ???");
 439     return;
 440   }
 441
 442   if (!fInitialized)
 443     ExecOnce(); // Get tracks, jets from arrays if not already given + Init Histos
 444
 445   // Get run number
 446   TString tmpRunNum("");
 447   tmpRunNum += event->GetRunNumber();
 448   // Additional cuts
 449   FillHistogram(tmpRunNum.Data(), "hNumberEvents", 0.5); // number of events before manual cuts
 450
 451   if(!fHelperClass->IsVertexSelected2013pA(event))
 452     return;
 453
 454   FillHistogram(tmpRunNum.Data(), "hNumberEvents", 1.5); // number of events after manual cuts
 455
 456   #ifdef DEBUGMODE
 457     AliInfo("Calculate()::Init done.");
 458   #endif
 459
 460   ////////////////////// NOTE: Get Centrality, (Leading)Signal jets
 461
 462   // Get centrality (V0A)
 463   AliCentrality* tmpCentrality = NULL;
 464   tmpCentrality = event->GetCentrality();
 465   Double_t centralityPercentileV0A = 0.0;
 466   Double_t centralityPercentileV0C = 0.0;
 467   Double_t centralityPercentileV0M = 0.0;
 468   if (tmpCentrality != NULL)
 469   {
 470     centralityPercentileV0A = tmpCentrality->GetCentralityPercentile("V0A");
 471     centralityPercentileV0C = tmpCentrality->GetCentralityPercentile("V0C");
 472     centralityPercentileV0M = tmpCentrality->GetCentralityPercentile("V0M");
 473   }
 474   // Get jets
 475   if (fAnalyzeJets)
 476     GetSignalJets();
 477
 478   #ifdef DEBUGMODE
 479     AliInfo("Calculate()::Centrality&SignalJets done.");
 480   #endif
 481   ////////////////////// NOTE: Pythia histograms
 482   if(fAnalyzePythia)
 483   {
 484     FillHistogram(tmpRunNum.Data(), "hPythiaPtHard", GetPtHard());
 485     FillHistogram(tmpRunNum.Data(), "hPythiaNTrials", GetPtHardBin()+0.1, fTrials);
 486     FillHistogram(tmpRunNum.Data(), "hPythiaXSection", GetPtHardBin()+0.1, fCrossSection);
 487
 488     #ifdef DEBUGMODE
 489       AliInfo("Calculate()::Pythia done.");
 490     #endif
 491   }
 492
 493   ////////////////////// NOTE: Track & QA histograms
 494
 495   if (fAnalyzeQA)
 496   {
 497     FillHistogram(tmpRunNum.Data(), "hVertexX",event->GetPrimaryVertex()->GetX());
 498     FillHistogram(tmpRunNum.Data(), "hVertexY",event->GetPrimaryVertex()->GetY());
 499     FillHistogram(tmpRunNum.Data(), "hVertexXY",event->GetPrimaryVertex()->GetX(), event->GetPrimaryVertex()->GetY());
 500     FillHistogram(tmpRunNum.Data(), "hVertexZ",event->GetPrimaryVertex()->GetZ());
 501     FillHistogram(tmpRunNum.Data(), "hVertexR",TMath::Sqrt(event->GetPrimaryVertex()->GetX()*event->GetPrimaryVertex()->GetX() + event->GetPrimaryVertex()->GetY()*event->GetPrimaryVertex()->GetY()));
 502     FillHistogram(tmpRunNum.Data(), "hCentralityV0M",centralityPercentileV0M);
 503     FillHistogram(tmpRunNum.Data(), "hCentralityV0A",centralityPercentileV0A);
 504     FillHistogram(tmpRunNum.Data(), "hCentralityV0C",centralityPercentileV0C);
 505
 506     Int_t trackCountAcc = 0;
 507     Int_t nTracks = fTrackArray->GetEntries();
 508     for (Int_t i = 0; i < nTracks; i++)
 509     {
 510       AliVTrack* track = static_cast<AliVTrack*>(fTrackArray->At(i));
 511       if (IsTrackInAcceptance(track))
 512       {
 513         FillHistogram(tmpRunNum.Data(), "hTrackPhiEta", track->Phi(),track->Eta(), 1);
 514         FillHistogram(tmpRunNum.Data(), "hTrackPt", track->Pt());
 515         if(track->Eta() >= 0)
 516           FillHistogram(tmpRunNum.Data(), "hTrackPtPosEta", track->Pt());
 517         else
 518           FillHistogram(tmpRunNum.Data(), "hTrackPtNegEta", track->Pt());
 519
 520         FillHistogram(tmpRunNum.Data(), "hTrackEta", track->Eta());
 521         FillHistogram(tmpRunNum.Data(), "hTrackPhi", track->Phi());
 522         FillHistogram(tmpRunNum.Data(), "hTrackPhiLabel", track->Phi(), track->GetLabel());
 523         for(Int_t j=0;j<20;j++)
 524           if(track->Pt() > j)
 525             FillHistogram(tmpRunNum.Data(), "hTrackPhiPtCut", track->Phi(), track->Pt());
 526
 527         FillHistogram(tmpRunNum.Data(), "hTrackCharge", track->Charge());
 528         trackCountAcc++;
 529       }
 530     }
 531     FillHistogram(tmpRunNum.Data(), "hTrackCountAcc", trackCountAcc, centralityPercentileV0M);
 532
 533   }
 534   #ifdef DEBUGMODE
 535     AliInfo("Calculate()::QA done.");
 536   #endif
 537
 538   ////////////////////// NOTE: Jet analysis and calculations
 539
 540   if (fAnalyzeJets)
 541   {
 542     FillHistogram(tmpRunNum.Data(), "hJetCountAll", fJetArray->GetEntries());
 543     FillHistogram(tmpRunNum.Data(), "hJetCountAccepted", fNumberSignalJets);
 544     // SIGNAL JET ANALYSIS
 545     for (Int_t i = 0; i<fNumberSignalJets; i++)
 546     {
 547       AliEmcalJet* tmpJet = fSignalJets[i];
 548
 549       FillHistogram(tmpRunNum.Data(), "hJetArea", tmpJet->Area());
 550       FillHistogram(tmpRunNum.Data(), "hJetAreaVsPt", tmpJet->Area(), tmpJet->Pt());
 551       FillHistogram(tmpRunNum.Data(), "hJetPt", tmpJet->Pt());
 552       FillHistogram(tmpRunNum.Data(), "hJetPtVsConstituentCount", tmpJet->Pt(),tmpJet->GetNumberOfTracks());
 553       FillHistogram(tmpRunNum.Data(), "hJetPhiEta", tmpJet->Phi(),tmpJet->Eta());
 554
 555
 556     }
 557   } //endif AnalyzeJets
 558
 559   #ifdef DEBUGMODE
 560     AliInfo("Calculate()::Jets done.");
 561   #endif
 562 }
 563
 564 //________________________________________________________________________
 565 Bool_t AliAnalysisTaskQualityAssurancePA::Notify()
 566 {
 567   // Implemented Notify() to read the cross sections
 568   // and number of trials from pyxsec.root
 569   //
 570   #ifdef DEBUGMODE
 571     AliInfo("Notify started.");
 572   #endif
 573
 574   if(fAnalyzePythia)
 575   {
 576     TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 577     TFile *currFile = tree->GetCurrentFile();
 578
 579     TString file(currFile->GetName());
 580
 581     if(file.Contains("root_archive.zip#")){
 582       Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
 583       Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
 584       file.Replace(pos+1,20,"");
 585     }
 586     else {
 587       // not an archive take the basename....
 588       file.ReplaceAll(gSystem->BaseName(file.Data()),"");
 589     }
 590
 591     TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
 592     if(!fxsec){
 593       // next trial fetch the histgram file
 594       fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
 595       if(!fxsec){
 596           // not a severe condition but inciate that we have no information
 597         return kFALSE;
 598       }
 599       else{
 600         // find the tlist we want to be independtent of the name so use the Tkey
 601         TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
 602         if(!key){
 603           fxsec->Close();
 604           return kFALSE;
 605         }
 606         TList *list = dynamic_cast<TList*>(key->ReadObj());
 607         if(!list){
 608           fxsec->Close();
 609           return kFALSE;
 610         }
 611         fCrossSection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 612         fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 613         fxsec->Close();
 614       }
 615     } // no tree pyxsec.root
 616     else {
 617       TTree *xtree = (TTree*)fxsec->Get("Xsection");
 618       if(!xtree){
 619         fxsec->Close();
 620         return kFALSE;
 621       }
 622       UInt_t   ntrials  = 0;
 623       Double_t  xsection  = 0;
 624       xtree->SetBranchAddress("xsection",&xsection);
 625       xtree->SetBranchAddress("ntrials",&ntrials);
 626       xtree->GetEntry(0);
 627       fTrials = ntrials;
 628       fCrossSection = xsection;
 629       fxsec->Close();
 630     }
 631     #ifdef DEBUGMODE
 632       AliInfo("Notify ended.");
 633     #endif
 634   }
 635   return kTRUE;
 636 }
 637
 638 //________________________________________________________________________
 639 inline void AliAnalysisTaskQualityAssurancePA::FillHistogram(const char* runNumber, const char * key, Double_t x)
 640 {
 641   if(strcmp(runNumber,"ALL") != 0)
 642     FillHistogram("ALL", key, x);
 643
 644   TH1* tmpHist = static_cast<TH1*>(fOutputList->FindObject(GetHistoName(runNumber, key)));
 645   if(!tmpHist)
 646   {
 647     AliError(Form("Cannot find histogram <%s> ",GetHistoName(runNumber, key))) ;
 648     return;
 649   }
 650
 651   tmpHist->Fill(x);
 652 }
 653
 654 //________________________________________________________________________
 655 inline void AliAnalysisTaskQualityAssurancePA::FillHistogram(const char* runNumber, const char * key, Double_t x, Double_t y)
 656 {
 657   if(strcmp(runNumber,"ALL") != 0)
 658     FillHistogram("ALL", key, x, y);
 659
 660   TH1* tmpHist = static_cast<TH1*>(fOutputList->FindObject(GetHistoName(runNumber, key)));
 661   if(!tmpHist)
 662   {
 663     AliError(Form("Cannot find histogram <%s> ",GetHistoName(runNumber, key)));
 664     return;
 665   }
 666
 667   if (tmpHist->IsA()->GetBaseClass("TH1"))
 668     static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
 669   else if (tmpHist->IsA()->GetBaseClass("TH2"))
 670     static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
 671 }
 672
 673 //________________________________________________________________________
 674 inline void AliAnalysisTaskQualityAssurancePA::FillHistogram(const char* runNumber, const char * key, Double_t x, Double_t y, Double_t add)
 675 {
 676   if(strcmp(runNumber,"ALL") != 0)
 677     FillHistogram("ALL", key, x, y, add);
 678
 679
 680   TH2* tmpHist = static_cast<TH2*>(fOutputList->FindObject(GetHistoName(runNumber, key)));
 681   if(!tmpHist)
 682   {
 683     AliError(Form("Cannot find histogram <%s> ",GetHistoName(runNumber, key)));
 684     return;
 685   }
 686
 687   tmpHist->Fill(x,y,add);
 688 }
 689 //________________________________________________________________________
 690 template <class T> T* AliAnalysisTaskQualityAssurancePA::AddHistogram1D(const char* runNumber, const char* name, const char* title, const char* options, Int_t xBins, Double_t xMin, Double_t xMax, const char* xTitle, const char* yTitle)
 691 {
 692   T* tmpHist = new T(GetHistoName(runNumber, name), GetHistoName(runNumber, title), xBins, xMin, xMax);
 693
 694   tmpHist->GetXaxis()->SetTitle(xTitle);
 695   tmpHist->GetYaxis()->SetTitle(yTitle);
 696   tmpHist->SetOption(options);
 697   tmpHist->SetMarkerStyle(kFullCircle);
 698   tmpHist->Sumw2();
 699
 700   fHistList->Add(tmpHist);
 701   fHistCount++;
 702
 703   return tmpHist;
 704 }
 705
 706 //________________________________________________________________________
 707 template <class T> T* AliAnalysisTaskQualityAssurancePA::AddHistogram2D(const char* runNumber, const char* name, const char* title, const char* options, Int_t xBins, Double_t xMin, Double_t xMax, Int_t yBins, Double_t yMin, Double_t yMax, const char* xTitle, const char* yTitle, const char* zTitle)
 708 {
 709   T* tmpHist = new T(GetHistoName(runNumber, name), GetHistoName(runNumber, title), xBins, xMin, xMax, yBins, yMin, yMax);
 710   tmpHist->GetXaxis()->SetTitle(xTitle);
 711   tmpHist->GetYaxis()->SetTitle(yTitle);
 712   tmpHist->GetZaxis()->SetTitle(zTitle);
 713   tmpHist->SetOption(options);
 714   tmpHist->SetMarkerStyle(kFullCircle);
 715   tmpHist->Sumw2();
 716
 717   fHistList->Add(tmpHist);
 718   fHistCount++;
 719
 720   return tmpHist;
 721 }
 722
 723 //________________________________________________________________________
 724 void AliAnalysisTaskQualityAssurancePA::Terminate(Option_t *)
 725 {
 726   PostData(1, fOutputList);
 727
 728   // Mandatory
 729   fOutputList = dynamic_cast<TList*> (GetOutputData(1)); // '1' refers to the output slot
 730   if (!fOutputList) {
 731     printf("ERROR: Output list not available\n");
 732     return;
 733   }
 734 }
 735
 736 //________________________________________________________________________
 737 AliAnalysisTaskQualityAssurancePA::~AliAnalysisTaskQualityAssurancePA()
 738 {
 739   // Destructor. Clean-up the output list, but not the histograms that are put inside
 740   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
 741   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 742     delete fOutputList;
 743   }
 744 }
 745
 746 //________________________________________________________________________
 747 void AliAnalysisTaskQualityAssurancePA::UserCreateOutputObjects()
 748 {
 749   // called once to create user defined output objects like histograms, plots etc.
 750   // and to put it on the output list.
 751   // Note: Saving to file with e.g. OpenFile(0) is must be before creating other objects.
 752
 753   fRandom = new TRandom3(0);
 754
 755   fOutputList = new TList();
 756   fOutputList->SetOwner(); // otherwise it produces leaks in merging
 757
 758   PostData(1, fOutputList);
 759 }
 760
 761 //________________________________________________________________________
 762 void AliAnalysisTaskQualityAssurancePA::UserExec(Option_t *)
 763 {
 764   #ifdef DEBUGMODE
 765     AliInfo("UserExec() started.");
 766   #endif
 767
 768   Calculate(InputEvent());
 769
 770   PostData(1, fOutputList);
 771 }