PWG2/FORWARD/analysis/AliFMDAnalysisTaskBFCorrelation.cxx

   1 #include <TROOT.h>
   2 #include <TSystem.h>
   3 #include <TInterpreter.h>
   4 #include <TChain.h>
   5 #include <TFile.h>
   6 #include <TList.h>
   7 #include <iostream>
   8 #include <string>
   9 #include <TCanvas.h>
  10 #include "TH1F.h"
  11 #include "TH2F.h"
  12 #include "TProfile.h"
  13 #include "AliFMDAnalysisTaskBFCorrelation.h"
  14 #include "AliAnalysisManager.h"
  15 #include "AliESDFMD.h"
  16 #include "AliESDEvent.h"
  17 #include "AliAODEvent.h"
  18 #include "AliAODHandler.h"
  19 #include "AliMCEventHandler.h"
  20 #include "AliStack.h"
  21 #include "AliLog.h"
  22 #include "AliESDVertex.h"
  23 #include "TMath.h"
  24 #include "AliFMDAnaParameters.h"
  25 //#include "AliFMDGeometry.h"
  26 #include "AliGenEventHeader.h"
  27 #include "AliGenPythiaEventHeader.h"
  28 #include "AliHeader.h"
  29 //#include "TDatabasePDG.h"
  30 //#include "TParticlePDG.h"
  31 #include "AliESDInputHandler.h"
  32
  33 ClassImp(AliFMDAnalysisTaskBFCorrelation)
  34
  35 AliFMDAnalysisTaskBFCorrelation::AliFMDAnalysisTaskBFCorrelation()
  36 : fDebug(0),
  37   fOutputList(0),
  38   fInputList(0),
  39   fInternalList(0),
  40   fVertexString(0x0),
  41   fStandalone(kTRUE),
  42   fEvent(0),
  43   fnBinsX(0),
  44   fXmin(0),
  45   fXmax(0),
  46   fnBinsY(0),
  47   fYmin(0),
  48   fYmax(0)
  49
  50 {
  51   // Default constructor
  52   DefineInput (0, TList::Class());
  53   DefineOutput(0, TList::Class());
  54 }
  55
  56 //_____________________________________________________________________
  57 AliFMDAnalysisTaskBFCorrelation::AliFMDAnalysisTaskBFCorrelation(const char* name, Bool_t SE):
  58   AliAnalysisTask(name,name),
  59   fDebug(0),
  60   fOutputList(0),
  61   fInputList(0),
  62   fInternalList(0),
  63   fVertexString(0x0),
  64   fStandalone(kTRUE),
  65   fEvent(0),
  66   fnBinsX(0),
  67   fXmin(0),
  68   fXmax(0),
  69   fnBinsY(0),
  70   fYmin(0),
  71   fYmax(0)
  72 {
  73   fStandalone = SE;
  74   if(fStandalone) {
  75     DefineInput (0, TList::Class());
  76     DefineInput(1, TObjString::Class());
  77     DefineOutput(0, TList::Class());
  78   }
  79 }
  80
  81 //_____________________________________________________________________
  82 void AliFMDAnalysisTaskBFCorrelation::CreateOutputObjects()
  83 {
  84
  85   // Nomenclature for naming histograms
  86   // SE    - Used in histogram names if used for containing Single Event information
  87   // Rebin - Used in histogram names if they have been rebinned for analysis purposes
  88
  89   // Setup the list for storing results, if it does not exist
  90
  91   if(!fOutputList) {
  92     fOutputList = new TList();
  93     fOutputList->SetName("BackgroundCorrected");
  94   }
  95
  96   // Setup the list for temporary storage during the analysis
  97
  98   if(!fInternalList) {
  99     fInternalList = new TList();
 100     fInternalList->SetName("InternalBFList");
 101   }
 102
 103   // Get the bounds for the input histogram.
 104   // This version is optimized for 200 bins in the eta range (-4, 6)
 105
 106   AliFMDAnaParameters *pars = AliFMDAnaParameters::Instance();
 107   TH2F *hDefault = (TH2F*)pars->GetBackgroundCorrection(1, 'I', 1);
 108
 109   fnBinsX = hDefault->GetNbinsX();
 110   fnBinsY = hDefault->GetNbinsY();
 111   fXmin   = hDefault->GetXaxis()->GetXmin();
 112   fXmax   = hDefault->GetXaxis()->GetXmax();
 113   fYmin   = hDefault->GetYaxis()->GetXmin();
 114   fYmax   = hDefault->GetYaxis()->GetXmax();
 115
 116   // Histogram to contain an event of MC data. Same dimension as ESD event
 117
 118   TH2F *hEtaPhiParticleMap = new TH2F("hEtaPhiParticleMap",
 119                                       "Distributions of MC particles (Eta,Phi)",
 120                                       fnBinsX, fXmin, fXmax,
 121                                       fnBinsY, fYmin, fYmax);
 122   hEtaPhiParticleMap->Sumw2();
 123
 124   fInternalList->Add(hEtaPhiParticleMap);
 125
 126   // Create histograms with same binning as input for Response analysis
 127   // and control histograms. One temporary histogram is also created to
 128   // avoid the new and delete call. Must be reset after use !
 129
 130   TH1D *hSEMultESD      = new TH1D("hSEMultESD", "Multiplicity", fnBinsX, fXmin, fXmax);
 131   TH1D *hSEMultMC       = new TH1D("hSEMultMC",  "Multiplicity", fnBinsX, fXmin, fXmax);
 132
 133   TH1D *hTemp           = new TH1D("hTemp", "Temporary histogram", fnBinsX, fXmin, fXmax);
 134
 135   fInternalList->Add(hSEMultESD);
 136   fInternalList->Add(hSEMultMC);
 137   fInternalList->Add(hTemp);
 138
 139   TH1F *hMultvsEtaESD   = new TH1F("hMultvsEtaESD", "Multiplicity vs Eta (ESD)", fnBinsX, fXmin, fXmax);
 140   TH1F *hMultvsEtaMC    = new TH1F("hMultvsEtaMC",  "Multiplicity vs Eta (MC)",  fnBinsX, fXmin, fXmax);
 141
 142   TH2F *hResponseMatrix = new TH2F("hResponseMatrix", "Response matrix", 151, -0.5, 150.5, 151, -0.5, 150.5);
 143
 144   fOutputList->Add(hMultvsEtaESD);
 145   fOutputList->Add(hMultvsEtaMC);
 146   fOutputList->Add(hResponseMatrix);
 147
 148   // Set up histograms for analysis and storage. 3 different binnings
 149
 150   for (Int_t i = 1; i <= 4; i++) {
 151     if (i == 3) continue;
 152
 153     Int_t nBinsX  = fnBinsX/(5*i);
 154
 155     // Histograms for the event-by-event analysis
 156
 157     TH1D *hSERebinMultESD        = new TH1D(Form("hSERebinMultESD_binning%d", i),
 158                                             Form("Multiplicity per event vs eta ESD (%d bins)", nBinsX),
 159                                             nBinsX, fXmin, fXmax);
 160
 161     TH1D *hSERebinMultMC         = new TH1D(Form("hSERebinMultMC_binning%d", i),
 162                                             Form("Multiplicity per event vs eta MC-truth (%d bins)", nBinsX),
 163                                             nBinsX, fXmin, fXmax);
 164
 165     TH1D *hSERebinMultMirrorESD  = new TH1D(Form("hSERebinMultMirrorESD_binning%d", i),
 166                                             Form("Multiplicity per event vs eta ESD Mirrored (%d bins)", nBinsX),
 167                                             nBinsX, fXmin, fXmax);
 168
 169     TH1D *hSERebinMultMirrorMC   = new TH1D(Form("hSERebinMultMirrorMC_binning%d", i),
 170                                             Form("Multiplicity per event vs eta MC-truth Mirrored (%d bins)", nBinsX),
 171                                             nBinsX, fXmin, fXmax);
 172
 173
 174     TH1D *hSERebinMultWESD       = new TH1D(Form("hSERebinMultWESD_binning%d", i),
 175                                             Form("Multiplicity per event vs eta ESD (%d bins)", nBinsX),
 176                                             nBinsX, fXmin, fXmax);
 177
 178     TH1D *hSERebinMultWMC        = new TH1D(Form("hSERebinMultWMC_binning%d", i),
 179                                             Form("Multiplicity per event vs eta MC-truth (%d bins)", nBinsX),
 180                                             nBinsX, fXmin, fXmax);
 181
 182     TH1D *hSERebinMultMirrorWESD = new TH1D(Form("hSERebinMultMirrorWESD_binning%d", i),
 183                                             Form("Multiplicity per event vs eta ESD Mirrored (%d bins)", nBinsX),
 184                                             nBinsX, fXmin, fXmax);
 185
 186     TH1D *hSERebinMultMirrorWMC  = new TH1D(Form("hSERebinMultMirrorWMC_binning%d", i),
 187                                             Form("Multiplicity per event vs eta MC-truth Mirrored (%d bins)", nBinsX),
 188                                             nBinsX, fXmin, fXmax);
 189
 190     fInternalList->Add(hSERebinMultESD);
 191     fInternalList->Add(hSERebinMultMC);
 192     fInternalList->Add(hSERebinMultMirrorESD);
 193     fInternalList->Add(hSERebinMultMirrorMC);
 194
 195     fInternalList->Add(hSERebinMultWESD);
 196     fInternalList->Add(hSERebinMultWMC);
 197     fInternalList->Add(hSERebinMultMirrorWESD);
 198     fInternalList->Add(hSERebinMultMirrorWMC);
 199
 200     // Histograms for storing the acummulated parameters.
 201
 202     TH1F *hRebinnESD        = new TH1F(Form("hRebinnESD_binning%d", i),
 203                                   Form("Counts ESD (%d bins)", nBinsX),
 204                                   nBinsX, fXmin, fXmax);
 205
 206     TH1F *hRebinnMirrorESD  = new TH1F(Form("hRebinnMirrorESD_binning%d", i),
 207                                   Form("Counts ESD Mirrored (%d bins)", nBinsX),
 208                                   nBinsX, fXmin, fXmax);
 209
 210     TH1F *hRebinn2ESD       = new TH1F(Form("hRebinn2ESD_binning%d", i),
 211                                   Form("Counts^2 ESD (%d bins)", nBinsX),
 212                                   nBinsX, fXmin, fXmax);
 213
 214     TH1F *hRebinn2MirrorESD = new TH1F(Form("hRebinn2MirrorESD_binning%d", i),
 215                                   Form("Counts^2 ESD Mirrored (%d bins)", nBinsX),
 216                                   nBinsX, fXmin, fXmax);
 217
 218     TH1F *hRebinnfbESD      = new TH1F(Form("hRebinnfbESD_binning%d", i),
 219                                   Form("Fwd*bwd ESD (%d bins)", nBinsX),
 220                                   nBinsX, fXmin, fXmax);
 221
 222
 223     TH1F *hRebinnMC         = new TH1F(Form("hRebinnMC_binning%d", i),
 224                                   Form("Counts MC (%d bins)", nBinsX),
 225                                   nBinsX, fXmin, fXmax);
 226
 227     TH1F *hRebinnMirrorMC   = new TH1F(Form("hRebinnMirrorMC_binning%d", i),
 228                                   Form("Counts MC Mirrored (%d bins)", nBinsX),
 229                                   nBinsX, fXmin, fXmax);
 230
 231     TH1F *hRebinn2MC        = new TH1F(Form("hRebinn2MC_binning%d", i),
 232                                   Form("Counts^2 MC (%d bins)", nBinsX),
 233                                   nBinsX, fXmin, fXmax);
 234
 235     TH1F *hRebinn2MirrorMC  = new TH1F(Form("hRebinn2MirrorMC_binning%d", i),
 236                                   Form("Counts^2 MC Mirrored (%d bins)", nBinsX),
 237                                   nBinsX, fXmin, fXmax);
 238
 239     TH1F *hRebinnfbMC       = new TH1F(Form("hRebinnfbMC_binning%d", i),
 240                                   Form("Fwd*bwd MC (%d bins)", nBinsX),
 241                                   nBinsX, fXmin, fXmax);
 242
 243     fOutputList->Add(hRebinnESD);
 244     fOutputList->Add(hRebinnMirrorESD);
 245     fOutputList->Add(hRebinn2ESD);
 246     fOutputList->Add(hRebinn2MirrorESD);
 247     fOutputList->Add(hRebinnfbESD);
 248
 249     fOutputList->Add(hRebinnMC);
 250     fOutputList->Add(hRebinnMirrorMC);
 251     fOutputList->Add(hRebinn2MC);
 252     fOutputList->Add(hRebinn2MirrorMC);
 253     fOutputList->Add(hRebinnfbMC);
 254
 255     // Histograms for storing the weights for the acummulated parameters.
 256
 257     TH1F *hRebinnWESD        = new TH1F(Form("hRebinnWESD_binning%d", i),
 258                                    Form("Counts ESD Weights (%d bins)", nBinsX),
 259                                    nBinsX, fXmin, fXmax);
 260
 261     TH1F *hRebinnMirrorWESD  = new TH1F(Form("hRebinnMirrorWESD_binning%d", i),
 262                                    Form("Counts ESD Weights (%d bins)", nBinsX),
 263                                    nBinsX, fXmin, fXmax);
 264
 265     TH1F *hRebinn2WESD       = new TH1F(Form("hRebinn2WESD_binning%d", i),
 266                                    Form("Counts^2 ESD Weights (%d bins)", nBinsX),
 267                                    nBinsX, fXmin, fXmax);
 268
 269     TH1F *hRebinn2MirrorWESD = new TH1F(Form("hRebinn2MirrorWESD_binning%d", i),
 270                                    Form("Counts^2 ESD Weights (%d bins)", nBinsX),
 271                                    nBinsX, fXmin, fXmax);
 272
 273     TH1F *hRebinnfbWESD      = new TH1F(Form("hRebinnfbWESD_binning%d", i),
 274                                    Form("Fwd*bwd ESD Weights (%d bins)", nBinsX),
 275                                    nBinsX, fXmin, fXmax);
 276
 277
 278     TH1F *hRebinnWMC         = new TH1F(Form("hRebinnWMC_binning%d", i),
 279                                    Form("Counts MC weights (%d bins)", nBinsX),
 280                                    nBinsX, fXmin, fXmax);
 281
 282     TH1F *hRebinnMirrorWMC   = new TH1F(Form("hRebinnMirrorWMC_binning%d", i),
 283                                    Form("Counts MC weights (%d bins)", nBinsX),
 284                                    nBinsX, fXmin, fXmax);
 285
 286     TH1F *hRebinn2WMC        = new TH1F(Form("hRebinn2WMC_binning%d", i),
 287                                    Form("Counts^2 MC weights (%d bins)", nBinsX),
 288                                    nBinsX, fXmin, fXmax);
 289
 290     TH1F *hRebinn2MirrorWMC  = new TH1F(Form("hRebinn2MirrorWMC_binning%d", i),
 291                                    Form("Counts^2 MC weights (%d bins)", nBinsX),
 292                                    nBinsX, fXmin, fXmax);
 293
 294     TH1F *hRebinnfbWMC       = new TH1F(Form("hRebinnfbWMC_binning%d", i),
 295                                    Form("Fwd*bwd MC weights (%d bins)", nBinsX),
 296                                    nBinsX, fXmin, fXmax);
 297
 298     fOutputList->Add(hRebinnWESD);
 299     fOutputList->Add(hRebinnMirrorWESD);
 300     fOutputList->Add(hRebinn2WESD);
 301     fOutputList->Add(hRebinn2MirrorWESD);
 302     fOutputList->Add(hRebinnfbWESD);
 303
 304     fOutputList->Add(hRebinnWMC);
 305     fOutputList->Add(hRebinnMirrorWMC);
 306     fOutputList->Add(hRebinn2WMC);
 307     fOutputList->Add(hRebinn2MirrorWMC);
 308     fOutputList->Add(hRebinnfbWMC);
 309
 310     // Histograms for the final result
 311     /*
 312     TH1D *hBFFcor = new TH1D(Form("hBFFcor_binning%d", i),
 313                              Form("Forward - backward correlations F (%d bins)", nBinsXBF),
 314                              nBinsXBF, 0, fXmax);
 315     hBFFcor->GetXaxis()->SetTitle("#eta");
 316     hBFFcor->GetXaxis()->SetTitle("b");
 317     TH1D *hBFBcor = new TH1D(Form("hBFBcor_binning%d", i),
 318                              Form("Forward - backward correlations B (%d bins)", nBinsXBF),
 319                              nBinsXBF, 0, fXmax);
 320     hBFBcor->GetXaxis()->SetTitle("#eta");
 321     hBFBcor->GetXaxis()->SetTitle("b");
 322
 323     TH1D *hBFFcor_MC = new TH1D(Form("hBFFcor_MC_binning%d", i),
 324                                 Form("Forward - backward correlations F (%d bins) MC", nBinsXBF),
 325                                 nBinsXBF, 0, fXmax);
 326     hBFFcor_MC->GetXaxis()->SetTitle("#eta");
 327     hBFFcor_MC->GetXaxis()->SetTitle("b");
 328     TH1D *hBFBcor_MC = new TH1D(Form("hBFBcor_MC_binning%d", i),
 329                                 Form("Forward - backward correlations B (%d bins) MC", nBinsXBF),
 330                                 nBinsXBF, 0, fXmax);
 331     hBFBcor_MC->GetXaxis()->SetTitle("#eta");
 332     hBFBcor_MC->GetXaxis()->SetTitle("b");
 333
 334     fOutputList->Add(hBFFcor);
 335     fOutputList->Add(hBFBcor);
 336     fOutputList->Add(hBFFcor_MC);
 337     fOutputList->Add(hBFBcor_MC);
 338     */
 339
 340     // Temporary histogram to avoid new-delete
 341
 342     TH1D* hRebinTemp = new TH1D(Form("hRebinTemp_binning%d", i),
 343                                   Form("Temporary histogram (%d bins)", nBinsX),
 344                                   nBinsX, fXmin, fXmax);
 345
 346     fInternalList->Add(hRebinTemp);
 347   }
 348 }
 349
 350 //_____________________________________________________________________
 351 void AliFMDAnalysisTaskBFCorrelation::ConnectInputData(Option_t */*option*/)
 352 {
 353   if(fStandalone) {
 354     fInputList   = (TList*)GetInputData(0);
 355     fVertexString = (TObjString*)GetInputData(1);
 356   }
 357 }
 358
 359 //_____________________________________________________________________
 360 void AliFMDAnalysisTaskBFCorrelation::Exec(Option_t */*option*/) {
 361
 362   fEvent++;
 363   //if (fEvent % 1000 == 0)
 364   //  std::cout << "Event # " << fEvent << std::endl;
 365
 366   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 367
 368   fVertexString = (TObjString*)fInputList->At(0);
 369
 370   //Int_t vtxbin = fVertexString->GetString().Atoi();
 371   //if (vtxbin != 5) return;
 372
 373   ProjectAndMirror("ESD");
 374   CalculateValues("ESD");
 375
 376   if(pars->GetProcessPrimary()) {
 377
 378     ProcessPrimary();
 379
 380     CreateResponseMatrix();
 381   }
 382
 383
 384   if(fStandalone) {
 385     PostData(0, fOutputList);
 386   }
 387 }
 388
 389 //_____________________________________________________________________
 390 void AliFMDAnalysisTaskBFCorrelation::ProjectAndMirror(TString sType) {
 391
 392   sType.ToUpper();
 393
 394   if (!sType.Contains("ESD") && !sType.Contains("MC")) {
 395     std::cout << "Wrong type specification for 'ProjectAndMirror'" << std::endl;
 396     return;
 397   }
 398
 399   // Get Single Event histograms for storing hits without rebinning
 400
 401   TH1D *hMult = dynamic_cast<TH1D*>(fInternalList->FindObject(Form("hSEMult%s", sType.Data())));
 402   hMult->Reset();
 403
 404   // Create generic names for retrieving histograms
 405
 406   TList *list = 0; // List for getting either ESD or MC "hit map"
 407
 408   TString sEtaPhiMap;
 409
 410   if (sType.Contains("ESD")) {
 411
 412     list = fInputList;
 413
 414     sEtaPhiMap = "dNdetadphiHistogramSPDTrVtx";
 415   }
 416
 417   if (sType.Contains("MC")) {
 418
 419     list = fInternalList;
 420
 421     sEtaPhiMap = "hEtaPhiParticleMap";
 422   }
 423
 424   TString sMult("hSERebinMult");
 425   TString sMultMirror("hSERebinMultMirror");
 426   TString sMultW("hSERebinMultW");
 427   TString sMultMirrorW("hSERebinMultMirrorW");
 428
 429   sType += "_binning";
 430
 431   sMult        += sType;
 432   sMultMirror  += sType;
 433   sMultW       += sType;
 434   sMultMirrorW += sType;
 435
 436   // Get the 2D histogram containing the particles of the event being analyzed
 437
 438   TH2F *hEtaPhiMap = dynamic_cast<TH2F*>(list->FindObject(sEtaPhiMap));
 439
 440   TH1D *hProjection = hEtaPhiMap->ProjectionX("hTemporary");
 441   hMult->Add(hProjection);
 442
 443   // Loop over the 3 binnings
 444
 445   for (Int_t i = 1; i<=4; i++) {
 446     if (i == 3) continue;
 447
 448     // Project the 2D "hit map" onto the eta-axis and rebin
 449
 450     TH1D *hProjRebin = hEtaPhiMap->ProjectionX("hProjRebin");
 451     hProjRebin->Rebin(i*5);
 452
 453     // Retrieve the histograms to store the Singe Event information and reset
 454
 455     TH1D *hSEMult        = static_cast<TH1D*>(fInternalList->FindObject(Form("%s%d", sMult.Data(),        i)));
 456     hSEMult->Reset();
 457
 458     TH1D *hSEMultMirror  = static_cast<TH1D*>(fInternalList->FindObject(Form("%s%d", sMultMirror.Data(),  i)));
 459     hSEMultMirror->Reset();
 460
 461     TH1D *hSEMultW       = static_cast<TH1D*>(fInternalList->FindObject(Form("%s%d", sMultW.Data(),       i)));
 462     hSEMultW->Reset();
 463
 464     TH1D *hSEMultMirrorW = static_cast<TH1D*>(fInternalList->FindObject(Form("%s%d", sMultMirrorW.Data(), i)));
 465     hSEMultMirrorW->Reset();
 466
 467     // Fill the histograms with the Single Event information
 468
 469     hSEMult->Add(hProjRebin);
 470
 471     for (Int_t bin = 1; bin <= hSEMult->GetNbinsX(); bin++) {
 472
 473       hSEMultMirror->SetBinContent(hSEMultMirror->FindBin(-hProjRebin->GetBinCenter(bin)), hProjRebin->GetBinContent(bin));
 474       hSEMultMirror->SetBinError(hSEMultMirror->FindBin(-hProjRebin->GetBinCenter(bin)), hProjRebin->GetBinError(bin));
 475
 476       // hMultDist->Fill(bin, hMultTemp->GetBinContent(bin));
 477     }
 478     hProjRebin->Delete();
 479   }
 480 }
 481
 482 //_____________________________________________________________________
 483 void AliFMDAnalysisTaskBFCorrelation::CalculateValues(TString sType) {
 484
 485   sType.ToUpper();
 486
 487   if (!sType.Contains("ESD") && !sType.Contains("MC")) {
 488     std::cout << "Wrong type specification for 'CalculateValues'" << std::endl;
 489     return;
 490   }
 491
 492   TString sMult("hSERebinMult");
 493   TString sMultMirror("hSERebinMultMirror");
 494   //  TString_t *sMultW("hSEMultW");
 495   //  TString_t *sMultMirrorW("hSEMultMirrorW");
 496
 497   TString sn("hRebinn");
 498   TString snMirror("hRebinnMirror");
 499   TString sn2("hRebinn2");
 500   TString sn2Mirror("hRebinn2Mirror");
 501   TString snfb("hRebinnfb");
 502
 503   sType += "_binning";
 504
 505   sMult       += sType;
 506   sMultMirror += sType;
 507
 508   sn          += sType;
 509   snMirror    += sType;
 510   sn2         += sType;
 511   sn2Mirror   += sType;
 512   snfb        += sType;
 513
 514   for (Int_t i = 1; i <= 4; i++) {
 515     if (i == 3) continue;
 516
 517     TH1D *hSEMult          = (TH1D*)fInternalList->FindObject(Form("%s%d", sMult.Data(),       i));
 518     TH1D *hSEMultMirror    = (TH1D*)fInternalList->FindObject(Form("%s%d", sMultMirror.Data(), i));
 519     /*
 520     TH1D *hSEMultW         = (TH1D*)fInternalList->FindObject(Form("%s%d", cMultW,       cType, i));
 521     TH1D *hSEMultMirrorW   = (TH1D*)fInternalList->FindObject(Form("%s%d", cMultMirrorW, cType, i));
 522     */
 523
 524     TH1F *hn               = (TH1F*)fOutputList->FindObject(Form("%s%d", sn.Data(),        i));
 525     TH1F *hnMirror         = (TH1F*)fOutputList->FindObject(Form("%s%d", snMirror.Data(),  i));
 526     TH1F *hn2              = (TH1F*)fOutputList->FindObject(Form("%s%d", sn2.Data(),       i));
 527     TH1F *hn2Mirror        = (TH1F*)fOutputList->FindObject(Form("%s%d", sn2Mirror.Data(), i));
 528     TH1F *hnfb             = (TH1F*)fOutputList->FindObject(Form("%s%d", snfb.Data(),      i));
 529     /*
 530     TH1F *hnW              = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cn, cType, i));
 531     TH1F *hnMirrorW        = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cnMirror, cType, i));
 532     TH1F *hn2W             = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cn2, cType, i));
 533     TH1F *hn2MirrorW       = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cn2Mirror, cType, i));
 534     TH1F *hnfbW            = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cnfb, cType, i));
 535     */
 536     TH1D *hTemp            = (TH1D*)fInternalList->FindObject(Form("hRebinTemp_binning%d", i));
 537
 538     hn->Add(hSEMult);
 539
 540     hnMirror->Add(hSEMultMirror);
 541
 542     hTemp->Reset();
 543     hTemp->Add(hSEMult);
 544     hTemp->Multiply(hSEMult);
 545     hn2->Add(hTemp);
 546
 547     hTemp->Reset();
 548     hTemp->Add(hSEMultMirror);
 549     hTemp->Multiply(hSEMultMirror);
 550     hn2Mirror->Add(hTemp);
 551
 552     hSEMultMirror->Multiply(hSEMult);
 553     hnfb->Add(hSEMultMirror);
 554   }
 555 }
 556
 557 //_____________________________________________________________________
 558 void AliFMDAnalysisTaskBFCorrelation::MultiplicityVsEta(TString sType) {
 559
 560   sType.ToUpper();
 561
 562   if (!sType.Contains("ESD") && !sType.Contains("MC")) {
 563     std::cout << "Wrong type specification for 'MultiplicityVsEta'" << std::endl;
 564     return;
 565   }
 566 }
 567
 568
 569 //_____________________________________________________________________
 570 void AliFMDAnalysisTaskBFCorrelation::CreateResponseMatrix() {
 571
 572   TH2F *hResponseMatrix = (TH2F*)fOutputList->FindObject("hResponseMatrix");
 573
 574   TH1D *hSEMultESD = (TH1D*)fInternalList->FindObject("hSEMultESD");
 575   TH1D *hSEMultMC  = (TH1D*)fInternalList->FindObject("hSEMultMC");
 576
 577   Float_t HitsESD = 0;
 578   Float_t HitsMC  = 0;
 579
 580   for (Int_t bin = 1; bin<= fnBinsX; bin++) {
 581
 582     if ((hSEMultMC->GetBinLowEdge(bin)  > -3.4) &&
 583         (hSEMultMC->GetBinLowEdge(bin+1) <  5.)) {
 584
 585       HitsESD += hSEMultESD->GetBinContent(bin);
 586       HitsMC  += hSEMultMC->GetBinContent(bin);
 587     }
 588   }
 589
 590   hResponseMatrix->Fill(HitsMC, HitsESD);
 591 }
 592
 593 //_____________________________________________________________________
 594 void AliFMDAnalysisTaskBFCorrelation::Terminate(Option_t */*option*/) {
 595   /*
 596   std::cout << "Terminating !" << std::endl;
 597
 598   TH1I *hnEvents   = (TH1I*)fOutputList->FindObject("nEvents");
 599   TH1I *hnMCEvents = (TH1I*)fOutputList->FindObject("nEvents");
 600
 601   Int_t nEvents   = hnEvents->GetEntries();
 602   Int_t nMCEvents = hnMCEvents->GetEntries();
 603
 604   for (Int_t i = 1; i <= 4; i++) {
 605     if (i == 3) continue;
 606
 607     TH1F *hnFnB = (TH1F*)fOutputList->FindObject(Form("hnFnB_binning%d", i));
 608     TH1F *hnF2  = (TH1F*)fOutputList->FindObject(Form("hnF2_binning%d", i));
 609     TH1F *hnB2  = (TH1F*)fOutputList->FindObject(Form("hnB2_binning%d", i));
 610     TH1F *hnF   = (TH1F*)fOutputList->FindObject(Form("hnF_binning%d", i));
 611     TH1F *hnB   = (TH1F*)fOutputList->FindObject(Form("hnB_binning%d", i));
 612
 613     TH1F *hnFnB_MC = (TH1F*)fOutputList->FindObject(Form("hnFnB_MC_binning%d", i));
 614     TH1F *hnF2_MC  = (TH1F*)fOutputList->FindObject(Form("hnF2_MC_binning%d", i));
 615     TH1F *hnB2_MC  = (TH1F*)fOutputList->FindObject(Form("hnB2_MC_binning%d", i));
 616     TH1F *hnF_MC   = (TH1F*)fOutputList->FindObject(Form("hnF_MC_binning%d", i));
 617     TH1F *hnB_MC   = (TH1F*)fOutputList->FindObject(Form("hnB_MC_binning%d", i));
 618
 619     hnFnB->Scale(1./Float_t(nEvents));
 620     hnF2->Scale(1./Float_t(nEvents));
 621     hnB2->Scale(1./Float_t(nEvents));
 622     hnF->Scale(1./Float_t(nEvents));
 623     hnB->Scale(1./Float_t(nEvents));
 624
 625     hnFnB_MC->Scale(1./Float_t(nMCEvents));
 626     hnF2_MC->Scale(1./Float_t(nMCEvents));
 627     hnB2_MC->Scale(1./Float_t(nMCEvents));
 628     hnF_MC->Scale(1./Float_t(nMCEvents));
 629     hnB_MC->Scale(1./Float_t(nMCEvents));
 630
 631     for (Int_t bin = 1; bin <= hnFnB->GetNbinsX(); bin++) {
 632
 633       Double_t nFnBav = hnFnB->GetBinContent(bin);
 634       Double_t nF2av  = hnF2->GetBinContent(bin);
 635       Double_t nB2av  = hnB2->GetBinContent(bin);
 636       Double_t nFav   = hnF->GetBinContent(bin);
 637       Double_t nBav   = hnB->GetBinContent(bin);
 638
 639       Double_t nFnB_MCav = hnFnB_MC->GetBinContent(bin);
 640       Double_t nF2_MCav  = hnF2_MC->GetBinContent(bin);
 641       Double_t nB2_MCav  = hnB2_MC->GetBinContent(bin);
 642       Double_t nF_MCav   = hnF_MC->GetBinContent(bin);
 643       Double_t nB_MCav   = hnB_MC->GetBinContent(bin);
 644
 645       Double_t bF = ((nF2av-nFav*nFav) == 0 ? 0. : (nFnBav-nFav*nBav)/(nF2av-nFav*nFav));
 646       Double_t bB = ((nF2av-nFav*nFav) == 0 ? 0. : (nFnBav-nFav*nBav)/(nB2av-nBav*nBav));
 647
 648       Double_t bF_MC = ((nF2_MCav-nF_MCav*nF_MCav) == 0 ? 0. :
 649                         (nFnB_MCav-nF_MCav*nB_MCav)/(nF2_MCav-nF_MCav*nF_MCav));
 650       Double_t bB_MC = ((nF2_MCav-nF_MCav*nF_MCav) == 0 ? 0. :
 651                         (nFnB_MCav-nF_MCav*nB_MCav)/(nB2_MCav-nB_MCav*nB_MCav));
 652
 653       TH1D *hBFFcor = (TH1D*)fOutputList->FindObject(Form("hBFFcor_binning%d", i));
 654       TH1D *hBFBcor = (TH1D*)fOutputList->FindObject(Form("hBFBcor_binning%d", i));
 655       TH1D *hBFFcor_MC = (TH1D*)fOutputList->FindObject(Form("hBFFcor_MC_binning%d", i));
 656       TH1D *hBFBcor_MC = (TH1D*)fOutputList->FindObject(Form("hBFBcor_MC_binning%d", i));
 657
 658       hBFFcor->SetBinContent(bin, bF);
 659       hBFBcor->SetBinContent(bin, bB);
 660       hBFFcor_MC->SetBinContent(bin, bF_MC);
 661       hBFBcor_MC->SetBinContent(bin, bB_MC);
 662     }
 663     }*/
 664 }
 665
 666 //_____________________________________________________________________
 667 void AliFMDAnalysisTaskBFCorrelation::ProcessPrimary() {
 668
 669   AliMCEventHandler* eventHandler =
 670     dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()
 671                                       ->GetMCtruthEventHandler());
 672   if (!eventHandler) return;
 673
 674   AliMCEvent* mcEvent = eventHandler->MCEvent();
 675   if(!mcEvent) return;
 676
 677   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 678
 679   AliMCParticle* particle = 0;
 680   AliStack* stack = mcEvent->Stack();
 681
 682   //TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
 683   AliHeader* header            = mcEvent->Header();
 684   AliGenEventHeader* genHeader = header->GenEventHeader();
 685   /*
 686   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 687
 688   if (!pythiaGenHeader) {
 689     std::cout<<" no pythia header!"<<std::endl;
 690     return;
 691   }
 692
 693   Int_t pythiaType = pythiaGenHeader->ProcessType();
 694
 695   if(pythiaType==92||pythiaType==93){
 696     std::cout<<"single diffractive"<<std::endl;
 697     return;
 698   }
 699   if(pythiaType==94){
 700     std::cout<<"double diffractive"<<std::endl;
 701     return;
 702   }
 703   */
 704   TArrayF vertex;
 705   genHeader->PrimaryVertex(vertex);
 706   if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
 707     return;
 708   //  Double_t delta           = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 709   //  Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
 710   //  Int_t    vertexBin       = (Int_t)vertexBinDouble;
 711
 712   //if (vertexBin != 5) return;
 713
 714   //Bool_t firstTrack = kTRUE;
 715
 716   // we loop over the primaries only unless we need the hits (diagnostics running slowly)
 717   Int_t nTracks = stack->GetNprimary();
 718   //  if(pars->GetProcessHits())
 719   //  nTracks = stack->GetNtrack();
 720
 721   TH2F *hEtaPhiParticleMap = (TH2F*)fInternalList->FindObject("hEtaPhiParticleMap");
 722   hEtaPhiParticleMap->Reset();
 723
 724   for(Int_t i= 0; i < nTracks; i++) {
 725     particle = (AliMCParticle*) mcEvent->GetTrack(i);
 726     if(!particle) continue;
 727
 728     if((stack->IsPhysicalPrimary(i)) && (particle->Charge() != 0)) {
 729       hEtaPhiParticleMap->Fill(particle->Eta(), particle->Phi());
 730       //std::cout << "Hans er nederen" << std::endl;
 731     }
 732   }
 733   ProjectAndMirror("MC");
 734   CalculateValues("MC");
 735 }
 736
 737 //_____________________________________________________________________
 738 //
 739 //
 740 // EOF