PWGLF/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   if(!hMult) {
 403     AliWarning("no hist - returning");
 404     return;
 405   }
 406   hMult->Reset();
 407
 408   // Create generic names for retrieving histograms
 409
 410   TList *list = 0; // List for getting either ESD or MC "hit map"
 411
 412   TString sEtaPhiMap;
 413
 414   if (sType.Contains("ESD")) {
 415
 416     list = fInputList;
 417
 418     sEtaPhiMap = "dNdetadphiHistogramSPDTrVtx";
 419   }
 420
 421   if (sType.Contains("MC")) {
 422
 423     list = fInternalList;
 424
 425     sEtaPhiMap = "hEtaPhiParticleMap";
 426   }
 427
 428   TString sMult("hSERebinMult");
 429   TString sMultMirror("hSERebinMultMirror");
 430   TString sMultW("hSERebinMultW");
 431   TString sMultMirrorW("hSERebinMultMirrorW");
 432
 433   sType += "_binning";
 434
 435   sMult        += sType;
 436   sMultMirror  += sType;
 437   sMultW       += sType;
 438   sMultMirrorW += sType;
 439
 440   // Get the 2D histogram containing the particles of the event being analyzed
 441
 442   TH2F *hEtaPhiMap = dynamic_cast<TH2F*>(list->FindObject(sEtaPhiMap));
 443
 444   TH1D *hProjection = hEtaPhiMap->ProjectionX("hTemporary");
 445   hMult->Add(hProjection);
 446
 447   // Loop over the 3 binnings
 448
 449   for (Int_t i = 1; i<=4; i++) {
 450     if (i == 3) continue;
 451
 452     // Project the 2D "hit map" onto the eta-axis and rebin
 453
 454     TH1D *hProjRebin = hEtaPhiMap->ProjectionX("hProjRebin");
 455     hProjRebin->Rebin(i*5);
 456
 457     // Retrieve the histograms to store the Singe Event information and reset
 458
 459     TH1D *hSEMult        = static_cast<TH1D*>(fInternalList->FindObject(Form("%s%d", sMult.Data(),        i)));
 460     hSEMult->Reset();
 461
 462     TH1D *hSEMultMirror  = static_cast<TH1D*>(fInternalList->FindObject(Form("%s%d", sMultMirror.Data(),  i)));
 463     hSEMultMirror->Reset();
 464
 465     TH1D *hSEMultW       = static_cast<TH1D*>(fInternalList->FindObject(Form("%s%d", sMultW.Data(),       i)));
 466     hSEMultW->Reset();
 467
 468     TH1D *hSEMultMirrorW = static_cast<TH1D*>(fInternalList->FindObject(Form("%s%d", sMultMirrorW.Data(), i)));
 469     hSEMultMirrorW->Reset();
 470
 471     // Fill the histograms with the Single Event information
 472
 473     hSEMult->Add(hProjRebin);
 474
 475     for (Int_t bin = 1; bin <= hSEMult->GetNbinsX(); bin++) {
 476
 477       hSEMultMirror->SetBinContent(hSEMultMirror->FindBin(-hProjRebin->GetBinCenter(bin)), hProjRebin->GetBinContent(bin));
 478       hSEMultMirror->SetBinError(hSEMultMirror->FindBin(-hProjRebin->GetBinCenter(bin)), hProjRebin->GetBinError(bin));
 479
 480       // hMultDist->Fill(bin, hMultTemp->GetBinContent(bin));
 481     }
 482     hProjRebin->Delete();
 483   }
 484 }
 485
 486 //_____________________________________________________________________
 487 void AliFMDAnalysisTaskBFCorrelation::CalculateValues(TString sType) {
 488
 489   sType.ToUpper();
 490
 491   if (!sType.Contains("ESD") && !sType.Contains("MC")) {
 492     std::cout << "Wrong type specification for 'CalculateValues'" << std::endl;
 493     return;
 494   }
 495
 496   TString sMult("hSERebinMult");
 497   TString sMultMirror("hSERebinMultMirror");
 498   //  TString_t *sMultW("hSEMultW");
 499   //  TString_t *sMultMirrorW("hSEMultMirrorW");
 500
 501   TString sn("hRebinn");
 502   TString snMirror("hRebinnMirror");
 503   TString sn2("hRebinn2");
 504   TString sn2Mirror("hRebinn2Mirror");
 505   TString snfb("hRebinnfb");
 506
 507   sType += "_binning";
 508
 509   sMult       += sType;
 510   sMultMirror += sType;
 511
 512   sn          += sType;
 513   snMirror    += sType;
 514   sn2         += sType;
 515   sn2Mirror   += sType;
 516   snfb        += sType;
 517
 518   for (Int_t i = 1; i <= 4; i++) {
 519     if (i == 3) continue;
 520
 521     TH1D *hSEMult          = (TH1D*)fInternalList->FindObject(Form("%s%d", sMult.Data(),       i));
 522     TH1D *hSEMultMirror    = (TH1D*)fInternalList->FindObject(Form("%s%d", sMultMirror.Data(), i));
 523     /*
 524     TH1D *hSEMultW         = (TH1D*)fInternalList->FindObject(Form("%s%d", cMultW,       cType, i));
 525     TH1D *hSEMultMirrorW   = (TH1D*)fInternalList->FindObject(Form("%s%d", cMultMirrorW, cType, i));
 526     */
 527
 528     TH1F *hn               = (TH1F*)fOutputList->FindObject(Form("%s%d", sn.Data(),        i));
 529     TH1F *hnMirror         = (TH1F*)fOutputList->FindObject(Form("%s%d", snMirror.Data(),  i));
 530     TH1F *hn2              = (TH1F*)fOutputList->FindObject(Form("%s%d", sn2.Data(),       i));
 531     TH1F *hn2Mirror        = (TH1F*)fOutputList->FindObject(Form("%s%d", sn2Mirror.Data(), i));
 532     TH1F *hnfb             = (TH1F*)fOutputList->FindObject(Form("%s%d", snfb.Data(),      i));
 533     /*
 534     TH1F *hnW              = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cn, cType, i));
 535     TH1F *hnMirrorW        = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cnMirror, cType, i));
 536     TH1F *hn2W             = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cn2, cType, i));
 537     TH1F *hn2MirrorW       = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cn2Mirror, cType, i));
 538     TH1F *hnfbW            = (TH1F*)fOutputList->FindObject(Form("%sW%s_binning%d", cnfb, cType, i));
 539     */
 540     TH1D *hTemp            = (TH1D*)fInternalList->FindObject(Form("hRebinTemp_binning%d", i));
 541
 542     hn->Add(hSEMult);
 543
 544     hnMirror->Add(hSEMultMirror);
 545
 546     hTemp->Reset();
 547     hTemp->Add(hSEMult);
 548     hTemp->Multiply(hSEMult);
 549     hn2->Add(hTemp);
 550
 551     hTemp->Reset();
 552     hTemp->Add(hSEMultMirror);
 553     hTemp->Multiply(hSEMultMirror);
 554     hn2Mirror->Add(hTemp);
 555
 556     hSEMultMirror->Multiply(hSEMult);
 557     hnfb->Add(hSEMultMirror);
 558   }
 559 }
 560
 561 //_____________________________________________________________________
 562 void AliFMDAnalysisTaskBFCorrelation::MultiplicityVsEta(TString sType) {
 563
 564   sType.ToUpper();
 565
 566   if (!sType.Contains("ESD") && !sType.Contains("MC")) {
 567     std::cout << "Wrong type specification for 'MultiplicityVsEta'" << std::endl;
 568     return;
 569   }
 570 }
 571
 572
 573 //_____________________________________________________________________
 574 void AliFMDAnalysisTaskBFCorrelation::CreateResponseMatrix() {
 575
 576   TH2F *hResponseMatrix = (TH2F*)fOutputList->FindObject("hResponseMatrix");
 577
 578   TH1D *hSEMultESD = (TH1D*)fInternalList->FindObject("hSEMultESD");
 579   TH1D *hSEMultMC  = (TH1D*)fInternalList->FindObject("hSEMultMC");
 580
 581   Float_t HitsESD = 0;
 582   Float_t HitsMC  = 0;
 583
 584   for (Int_t bin = 1; bin<= fnBinsX; bin++) {
 585
 586     if ((hSEMultMC->GetBinLowEdge(bin)  > -3.4) &&
 587         (hSEMultMC->GetBinLowEdge(bin+1) <  5.)) {
 588
 589       HitsESD += hSEMultESD->GetBinContent(bin);
 590       HitsMC  += hSEMultMC->GetBinContent(bin);
 591     }
 592   }
 593
 594   hResponseMatrix->Fill(HitsMC, HitsESD);
 595 }
 596
 597 //_____________________________________________________________________
 598 void AliFMDAnalysisTaskBFCorrelation::Terminate(Option_t */*option*/) {
 599   /*
 600   std::cout << "Terminating !" << std::endl;
 601
 602   TH1I *hnEvents   = (TH1I*)fOutputList->FindObject("nEvents");
 603   TH1I *hnMCEvents = (TH1I*)fOutputList->FindObject("nEvents");
 604
 605   Int_t nEvents   = hnEvents->GetEntries();
 606   Int_t nMCEvents = hnMCEvents->GetEntries();
 607
 608   for (Int_t i = 1; i <= 4; i++) {
 609     if (i == 3) continue;
 610
 611     TH1F *hnFnB = (TH1F*)fOutputList->FindObject(Form("hnFnB_binning%d", i));
 612     TH1F *hnF2  = (TH1F*)fOutputList->FindObject(Form("hnF2_binning%d", i));
 613     TH1F *hnB2  = (TH1F*)fOutputList->FindObject(Form("hnB2_binning%d", i));
 614     TH1F *hnF   = (TH1F*)fOutputList->FindObject(Form("hnF_binning%d", i));
 615     TH1F *hnB   = (TH1F*)fOutputList->FindObject(Form("hnB_binning%d", i));
 616
 617     TH1F *hnFnB_MC = (TH1F*)fOutputList->FindObject(Form("hnFnB_MC_binning%d", i));
 618     TH1F *hnF2_MC  = (TH1F*)fOutputList->FindObject(Form("hnF2_MC_binning%d", i));
 619     TH1F *hnB2_MC  = (TH1F*)fOutputList->FindObject(Form("hnB2_MC_binning%d", i));
 620     TH1F *hnF_MC   = (TH1F*)fOutputList->FindObject(Form("hnF_MC_binning%d", i));
 621     TH1F *hnB_MC   = (TH1F*)fOutputList->FindObject(Form("hnB_MC_binning%d", i));
 622
 623     hnFnB->Scale(1./Float_t(nEvents));
 624     hnF2->Scale(1./Float_t(nEvents));
 625     hnB2->Scale(1./Float_t(nEvents));
 626     hnF->Scale(1./Float_t(nEvents));
 627     hnB->Scale(1./Float_t(nEvents));
 628
 629     hnFnB_MC->Scale(1./Float_t(nMCEvents));
 630     hnF2_MC->Scale(1./Float_t(nMCEvents));
 631     hnB2_MC->Scale(1./Float_t(nMCEvents));
 632     hnF_MC->Scale(1./Float_t(nMCEvents));
 633     hnB_MC->Scale(1./Float_t(nMCEvents));
 634
 635     for (Int_t bin = 1; bin <= hnFnB->GetNbinsX(); bin++) {
 636
 637       Double_t nFnBav = hnFnB->GetBinContent(bin);
 638       Double_t nF2av  = hnF2->GetBinContent(bin);
 639       Double_t nB2av  = hnB2->GetBinContent(bin);
 640       Double_t nFav   = hnF->GetBinContent(bin);
 641       Double_t nBav   = hnB->GetBinContent(bin);
 642
 643       Double_t nFnB_MCav = hnFnB_MC->GetBinContent(bin);
 644       Double_t nF2_MCav  = hnF2_MC->GetBinContent(bin);
 645       Double_t nB2_MCav  = hnB2_MC->GetBinContent(bin);
 646       Double_t nF_MCav   = hnF_MC->GetBinContent(bin);
 647       Double_t nB_MCav   = hnB_MC->GetBinContent(bin);
 648
 649       Double_t bF = ((nF2av-nFav*nFav) == 0 ? 0. : (nFnBav-nFav*nBav)/(nF2av-nFav*nFav));
 650       Double_t bB = ((nF2av-nFav*nFav) == 0 ? 0. : (nFnBav-nFav*nBav)/(nB2av-nBav*nBav));
 651
 652       Double_t bF_MC = ((nF2_MCav-nF_MCav*nF_MCav) == 0 ? 0. :
 653                         (nFnB_MCav-nF_MCav*nB_MCav)/(nF2_MCav-nF_MCav*nF_MCav));
 654       Double_t bB_MC = ((nF2_MCav-nF_MCav*nF_MCav) == 0 ? 0. :
 655                         (nFnB_MCav-nF_MCav*nB_MCav)/(nB2_MCav-nB_MCav*nB_MCav));
 656
 657       TH1D *hBFFcor = (TH1D*)fOutputList->FindObject(Form("hBFFcor_binning%d", i));
 658       TH1D *hBFBcor = (TH1D*)fOutputList->FindObject(Form("hBFBcor_binning%d", i));
 659       TH1D *hBFFcor_MC = (TH1D*)fOutputList->FindObject(Form("hBFFcor_MC_binning%d", i));
 660       TH1D *hBFBcor_MC = (TH1D*)fOutputList->FindObject(Form("hBFBcor_MC_binning%d", i));
 661
 662       hBFFcor->SetBinContent(bin, bF);
 663       hBFBcor->SetBinContent(bin, bB);
 664       hBFFcor_MC->SetBinContent(bin, bF_MC);
 665       hBFBcor_MC->SetBinContent(bin, bB_MC);
 666     }
 667     }*/
 668 }
 669
 670 //_____________________________________________________________________
 671 void AliFMDAnalysisTaskBFCorrelation::ProcessPrimary() {
 672
 673   AliMCEventHandler* eventHandler =
 674     dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()
 675                                       ->GetMCtruthEventHandler());
 676   if (!eventHandler) return;
 677
 678   AliMCEvent* mcEvent = eventHandler->MCEvent();
 679   if(!mcEvent) return;
 680
 681   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 682
 683   AliMCParticle* particle = 0;
 684   AliStack* stack = mcEvent->Stack();
 685
 686   //TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
 687   AliHeader* header            = mcEvent->Header();
 688   AliGenEventHeader* genHeader = header->GenEventHeader();
 689   /*
 690   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 691
 692   if (!pythiaGenHeader) {
 693     std::cout<<" no pythia header!"<<std::endl;
 694     return;
 695   }
 696
 697   Int_t pythiaType = pythiaGenHeader->ProcessType();
 698
 699   if(pythiaType==92||pythiaType==93){
 700     std::cout<<"single diffractive"<<std::endl;
 701     return;
 702   }
 703   if(pythiaType==94){
 704     std::cout<<"double diffractive"<<std::endl;
 705     return;
 706   }
 707   */
 708   TArrayF vertex;
 709   genHeader->PrimaryVertex(vertex);
 710   if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
 711     return;
 712   //  Double_t delta           = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
 713   //  Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
 714   //  Int_t    vertexBin       = (Int_t)vertexBinDouble;
 715
 716   //if (vertexBin != 5) return;
 717
 718   //Bool_t firstTrack = kTRUE;
 719
 720   // we loop over the primaries only unless we need the hits (diagnostics running slowly)
 721   Int_t nTracks = stack->GetNprimary();
 722   //  if(pars->GetProcessHits())
 723   //  nTracks = stack->GetNtrack();
 724
 725   TH2F *hEtaPhiParticleMap = (TH2F*)fInternalList->FindObject("hEtaPhiParticleMap");
 726   hEtaPhiParticleMap->Reset();
 727
 728   for(Int_t i= 0; i < nTracks; i++) {
 729     particle = (AliMCParticle*) mcEvent->GetTrack(i);
 730     if(!particle) continue;
 731
 732     if((stack->IsPhysicalPrimary(i)) && (particle->Charge() != 0)) {
 733       hEtaPhiParticleMap->Fill(particle->Eta(), particle->Phi());
 734       //std::cout << "Hans er nederen" << std::endl;
 735     }
 736   }
 737   ProjectAndMirror("MC");
 738   CalculateValues("MC");
 739 }
 740
 741 //_____________________________________________________________________
 742 //
 743 //
 744 // EOF