PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskQualityAssurancePA.cxx

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