PWGCF/Correlations/DPhi/FourierDecomposition/AliMuonEffMC.cxx

   1 // MUON track QA referring AliMuonEffMC.cxx
   2 // Author : Saehanseul Oh
   3
   4 #include "AliMuonEffMC.h"
   5
   6 #include <TList.h>
   7 #include <TH1D.h>
   8 #include <TH2F.h>
   9 #include <TH3F.h>
  10 #include <THn.h>
  11 #include <TChain.h>
  12 #include <TFile.h>
  13 #include <TParticle.h>
  14
  15 #include "AliStack.h"
  16 #include "AliAnalysisManager.h"
  17 #include "AliAnalysisTask.h"
  18 #include "AliESDVertex.h"
  19 #include "AliESDEvent.h"
  20 #include "AliAODEvent.h"
  21 #include "AliESDMuonTrack.h"
  22 #include "AliAODTrack.h"
  23 #include "AliESDVertex.h"
  24 #include "AliAODVertex.h"
  25 #include "AliCentrality.h"
  26 #include "AliVParticle.h"
  27 #include "AliMCParticle.h"
  28 #include "AliMCEvent.h"
  29 #include "AliAnalysisHelperJetTasks.h"
  30 #include "AliAODMCParticle.h"
  31
  32 using std::cout;
  33 using std::endl;
  34
  35 ClassImp(AliMuonEffMC)
  36
  37 //________________________________________________________________________
  38 AliMuonEffMC::AliMuonEffMC() :
  39   AliAnalysisTaskSE(), fESD(0), fAOD(0), fMC(0), fStack(0), fCentrality(99), fZVertex(99), fOutputList(0x0),
  40   fHEventStat(0), fHXsec(0), fHTrials(0), fHEvt(0x0), fIsMc(kTRUE), fIsPythia(kFALSE), fMDProcess(kFALSE), fFeynmanX(kFALSE), fScatFX(kFALSE), fZvProcess(kTRUE),
  41   fIsCutStudy(kFALSE), fIsFPM(kTRUE), fCentralityEstimator("V0M"), fNEtaBins(15), fNpTBins(50), fNCentBins(1), fNZvtxBins(1), fNPhiBins(12), fNPBins(150), fChiSquareNormCut(5.0),
  42   fHMuonParGen(0x0), fHMuonParGenSec(0x0), fHMuonParGenFPM(0x0), fHMuonParGenSecFPM(0x0), fHMuonParGenP(0x0), fHMuonDetGenP(0x0), fHMuonDetRecP(0x0),
  43   fHFXu(0), fHFXantiu(0), fHFXd(0), fHFXantid(0), fHFXg(0), fHFXetc(0),
  44   fHaFx(0), fHbFx(0), fHabFxRatio(0), fHabDeltaFx(0), fHabRelDeltaFx(0),
  45   fHZvRvPrim(0x0), fHZvRvSec(0x0), fHXvYvPrim(0x0), fHXvYvSec(0x0)
  46 {
  47   // Constructor
  48   //DefineInput(0, TChain::Class());
  49   //DefineOutput(1, TList::Class());
  50   for(Int_t i=0; i<6; i++)
  51   {
  52     fHMuonDetGen[i] = NULL;
  53     fHMuonDetRec[i] = NULL;
  54     fHHadDetRec[i] = NULL;
  55     fHSecDetRec[i] = NULL;
  56   }
  57   for(Int_t i=0; i<4; i++)
  58   {
  59     fHMuonParGenV[i] = NULL;
  60     fHMuonParGenSecV[i] = NULL;
  61     fHMuonDetRecV[i] = NULL;
  62     fHMuMotherGenPt[i] = NULL;
  63     fHMuMotherRecPt[i] = NULL;
  64     fHMuMotherGenPhi[i] = NULL;
  65     fHMuMotherRecPhi[i] = NULL;
  66     fHMuMotherGenEta[i] = NULL;
  67     fHMuMotherRecEta[i] = NULL;
  68     fHMuDCA[i] = NULL;
  69     fHMuMotherGenPtSec[i] = NULL;
  70     fHMuMotherRecPtSec[i] = NULL;
  71     fHMuMotherGenPhiSec[i] = NULL;
  72     fHMuMotherRecPhiSec[i] = NULL;
  73     fHMuMotherGenEtaSec[i] = NULL;
  74     fHMuMotherRecEtaSec[i] = NULL;
  75     fHMuDCASec[i] = NULL;
  76     fHMuZv[i] = NULL;
  77     fHMuRelZv[i] = NULL;
  78     for(Int_t j=0; j<3; j++)
  79     {
  80       fHMuMohterPhiDifGen[i][j] = NULL;
  81       fHMuMohterPhiDifRec[i][j] = NULL;
  82       fHMuMohterEtaDifGen[i][j] = NULL;
  83       fHMuMohterEtaDifRec[i][j] = NULL;
  84       fHMuMohterPhiDifGenSec[i][j] = NULL;
  85       fHMuMohterPhiDifRecSec[i][j] = NULL;
  86       fHMuMohterEtaDifGenSec[i][j] = NULL;
  87       fHMuMohterEtaDifRecSec[i][j] = NULL;
  88     }
  89   }
  90   for(Int_t i=0; i<2; i++)
  91   {
  92     for(Int_t j=0; j<3; j++)
  93     {
  94       fHFXmuonP[i][j] = NULL;
  95       fHFXmuonM[i][j] = NULL;
  96     }
  97   }
  98   for(Int_t i=0; i<3; i++)
  99   {
 100     fHabFxMu[i] = NULL;
 101     fHabFxRatioMu[i] = NULL;
 102     fHabDeltaFxMu[i] = NULL;
 103     fHabRelDeltaFxMu[i] = NULL;
 104   }
 105 }
 106
 107 //________________________________________________________________________
 108 AliMuonEffMC::AliMuonEffMC(const char *name) :
 109   AliAnalysisTaskSE(name), fESD(0), fAOD(0), fMC(0), fStack(0), fCentrality(99), fZVertex(99), fOutputList(0x0),
 110   fHEventStat(0), fHXsec(0), fHTrials(0), fHEvt(0x0),  fIsMc(kTRUE), fIsPythia(kFALSE), fMDProcess(kFALSE), fFeynmanX(kFALSE), fScatFX(kFALSE), fZvProcess(kTRUE),
 111   fIsCutStudy(kFALSE), fIsFPM(kTRUE), fCentralityEstimator("V0M"), fNEtaBins(15), fNpTBins(50), fNCentBins(1), fNZvtxBins(1), fNPhiBins(12), fNPBins(150), fChiSquareNormCut(5.0),
 112   fHMuonParGen(0x0), fHMuonParGenSec(0x0),  fHMuonParGenFPM(0x0), fHMuonParGenSecFPM(0x0), fHMuonParGenP(0x0), fHMuonDetGenP(0x0), fHMuonDetRecP(0x0),
 113   fHFXu(0), fHFXantiu(0), fHFXd(0),  fHFXantid(0), fHFXg(0), fHFXetc(0),
 114   fHaFx(0), fHbFx(0), fHabFxRatio(0), fHabDeltaFx(0), fHabRelDeltaFx(0),
 115   fHZvRvPrim(0x0), fHZvRvSec(0x0), fHXvYvPrim(0x0), fHXvYvSec(0x0)
 116 {
 117   // Constructor
 118   for(Int_t i=0; i<6; i++)
 119   {
 120     fHMuonDetGen[i] = NULL;
 121     fHMuonDetRec[i] = NULL;
 122     fHHadDetRec[i] = NULL;
 123     fHSecDetRec[i] = NULL;
 124   }
 125   for(Int_t i=0; i<4; i++)
 126   {
 127     fHMuonParGenV[i] = NULL;
 128     fHMuonParGenSecV[i] = NULL;
 129     fHMuonDetRecV[i] = NULL;
 130     fHMuMotherGenPt[i] = NULL;
 131     fHMuMotherRecPt[i] = NULL;
 132     fHMuMotherGenPhi[i] = NULL;
 133     fHMuMotherRecPhi[i] = NULL;
 134     fHMuMotherGenEta[i] = NULL;
 135     fHMuMotherRecEta[i] = NULL;
 136     fHMuDCA[i] = NULL;
 137     fHMuMotherGenPtSec[i] = NULL;
 138     fHMuMotherRecPtSec[i] = NULL;
 139     fHMuMotherGenPhiSec[i] = NULL;
 140     fHMuMotherRecPhiSec[i] = NULL;
 141     fHMuMotherGenEtaSec[i] = NULL;
 142     fHMuMotherRecEtaSec[i] = NULL;
 143     fHMuDCASec[i] = NULL;
 144     fHMuZv[i] = NULL;
 145     fHMuRelZv[i] = NULL;
 146     for(Int_t j=0; j<3; j++)
 147     {
 148       fHMuMohterPhiDifGen[i][j] = NULL;
 149       fHMuMohterPhiDifRec[i][j] = NULL;
 150       fHMuMohterEtaDifGen[i][j] = NULL;
 151       fHMuMohterEtaDifRec[i][j] = NULL;
 152       fHMuMohterPhiDifGenSec[i][j] = NULL;
 153       fHMuMohterPhiDifRecSec[i][j] = NULL;
 154       fHMuMohterEtaDifGenSec[i][j] = NULL;
 155       fHMuMohterEtaDifRecSec[i][j] = NULL;
 156     }
 157   }
 158   for(Int_t i=0; i<2; i++)
 159   {
 160     for(Int_t j=0; j<3; j++)
 161     {
 162       fHFXmuonP[i][j] = NULL;
 163       fHFXmuonM[i][j] = NULL;
 164     }
 165   }
 166   for(Int_t i=0; i<3; i++)
 167   {
 168     fHabFxMu[i] = NULL;
 169     fHabFxRatioMu[i] = NULL;
 170     fHabDeltaFxMu[i] = NULL;
 171     fHabRelDeltaFxMu[i] = NULL;
 172   }
 173   DefineInput(0, TChain::Class());
 174   DefineOutput(1, TList::Class());
 175 }
 176
 177 //________________________________________________________________________
 178 AliMuonEffMC::~AliMuonEffMC()
 179 {
 180   //Destructor
 181   if(fOutputList) delete fOutputList;
 182 }
 183
 184 //________________________________________________________________________
 185 void AliMuonEffMC::UserCreateOutputObjects()
 186 {
 187   // Create histograms
 188   // Called once (per slave on PROOF!)
 189   fOutputList = new TList();
 190   fOutputList->SetOwner(1);
 191
 192   fHEventStat = new TH1D("fHEventStat","Event statistics for analysis",18,0,18);
 193   fHEventStat->GetXaxis()->SetBinLabel(1,"Event");
 194   fHEventStat->GetXaxis()->SetBinLabel(2,"SelectedEvent");
 195   fHEventStat->GetXaxis()->SetBinLabel(3,"File");
 196   fHEventStat->GetXaxis()->SetBinLabel(4,"fSPHighpt");  //!Global Trigger Single plus High p_T
 197   fHEventStat->GetXaxis()->SetBinLabel(5,"fSPAllpt");   //!Global Trigger Single plus All p_T
 198   fHEventStat->GetXaxis()->SetBinLabel(6,"fSMLowpt");   //!Global Trigger Single minus Low p_T
 199   fHEventStat->GetXaxis()->SetBinLabel(7,"fSMHighpt");  //!Global Trigger Single minus High p_T
 200   fHEventStat->GetXaxis()->SetBinLabel(8,"fSMAllpt");   //!Global Trigger Single minus All p_T
 201   fHEventStat->GetXaxis()->SetBinLabel(9,"fSULowpt");   //!Global Trigger Single undefined Low p_T
 202   fHEventStat->GetXaxis()->SetBinLabel(10,"fSUHighpt"); //!Global Trigger Single undefined High p_T
 203   fHEventStat->GetXaxis()->SetBinLabel(11,"fSUAllpt");  //!Global Trigger Single undefined All p_T
 204   fHEventStat->GetXaxis()->SetBinLabel(12,"fUSLowpt");  //!Global Trigger UnlikeSign pair Low p_T
 205   fHEventStat->GetXaxis()->SetBinLabel(13,"fUSHighpt"); //!Global Trigger UnlikeSign pair High p_T
 206   fHEventStat->GetXaxis()->SetBinLabel(14,"fUSAllpt");  //!Global Trigger UnlikeSign pair All p_T
 207   fHEventStat->GetXaxis()->SetBinLabel(15,"fLSLowpt");  //!Global Trigger LikeSign pair pair Low p_T
 208   fHEventStat->GetXaxis()->SetBinLabel(16,"fLSHighpt"); //!Global Trigger LikeSign pair pair High p_T
 209   fHEventStat->GetXaxis()->SetBinLabel(17,"fLSAllpt");  //!Global Trigger LikeSign pair pair All p_T
 210   fHEventStat->GetXaxis()->SetBinLabel(18,"fSPLowpt");   //!Global Trigger Single plus Low p_T
 211   fOutputList->Add(fHEventStat);
 212
 213   if(fIsPythia)
 214   {
 215     fHXsec = new TH1F("fHXsec", "Cross section from pyxsec.root", 1, 0, 1);
 216     fHXsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 217     fOutputList->Add(fHXsec);
 218
 219     fHTrials = new TH1F("fHTrials", "Number of Trials", 1, 0, 1);
 220     fHTrials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 221     fOutputList->Add(fHTrials);
 222   }
 223
 224   fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", 30, -15, 15, 103, -2, 101);
 225   fOutputList->Add(fHEvt);
 226
 227   Int_t iTrackBin[6];
 228   Double_t* trackBins[6];
 229   const char* trackAxisTitle[6];
 230
 231   Int_t iTrackBinP[6];
 232   Double_t* trackBinsP[6];
 233   const char* trackAxisTitleP[6];
 234
 235   // eta
 236   Double_t etaBins[fNEtaBins+1];
 237   for(Int_t i=0; i<=fNEtaBins; i++) { etaBins[i] = (Double_t)(-4.0 + 1.5/fNEtaBins*i); }
 238   iTrackBin[0] = fNEtaBins;
 239   iTrackBinP[0] = fNEtaBins;
 240   trackBins[0] = etaBins;
 241   trackBinsP[0] = etaBins;
 242   trackAxisTitle[0] = "#eta";
 243   trackAxisTitleP[0] = "#eta";
 244
 245   // p_T
 246   Double_t pTBins[fNpTBins+1];
 247   for(Int_t i=0; i<=fNpTBins; i++) { pTBins[i] = (Double_t)(5.0/fNpTBins * i); }
 248   iTrackBin[1] = fNpTBins;
 249   trackBins[1] = pTBins;
 250   trackAxisTitle[1] = "p_{T} (GeV/c)";
 251
 252   // P
 253   Double_t PBins[fNPBins+1];
 254   for(Int_t i=0; i<=fNPBins; i++) { PBins[i] = (Double_t)(150.0/fNPBins * i); }
 255   iTrackBinP[1] = fNPBins;
 256   trackBinsP[1] = PBins;
 257   trackAxisTitleP[1] = "P (GeV/c)";
 258
 259   // centrality
 260   Double_t CentBins[fNCentBins+1];
 261   for (Int_t i=0; i<=fNCentBins; i++) { CentBins[i] = (Double_t)(100.0/fNCentBins * i); }
 262   iTrackBin[2] = fNCentBins;
 263   iTrackBinP[2] = fNCentBins;
 264   trackBins[2] = CentBins;
 265   trackBinsP[2] = CentBins;
 266   trackAxisTitle[2] = "Cent";
 267   trackAxisTitleP[2] = "Cent";
 268
 269   // Z-vertex
 270   Double_t ZvtxBins[fNZvtxBins+1];
 271   for(Int_t i=0; i<=fNZvtxBins; i++) { ZvtxBins[i] = (Double_t)(-10.0 + 20.0/fNZvtxBins * i); }
 272   iTrackBin[3] = fNZvtxBins;
 273   iTrackBinP[3] = fNZvtxBins;
 274   trackBins[3] = ZvtxBins;
 275   trackBinsP[3] = ZvtxBins;
 276   trackAxisTitle[3] = "Zvtx";
 277   trackAxisTitleP[3] = "Zvtx";
 278
 279   // phi
 280   Double_t phiBins[fNPhiBins+1];
 281   for(Int_t i=0; i<=fNPhiBins; i++) { phiBins[i] = (Double_t)(TMath::TwoPi()/fNPhiBins * i); }
 282   iTrackBin[4] = fNPhiBins;
 283   iTrackBinP[4] = fNPhiBins;
 284   trackBins[4] = phiBins;
 285   trackBinsP[4] = phiBins;
 286   trackAxisTitle[4] = "#phi";
 287   trackAxisTitleP[4] = "#phi";
 288
 289   // charge
 290   Double_t chargeBins[3] = {-10.0, 0, 10.0};
 291   iTrackBin[5] = 2;
 292   iTrackBinP[5] = 2;
 293   trackBins[5] = chargeBins;
 294   trackBinsP[5] = chargeBins;
 295   trackAxisTitle[5] = "charge";
 296   trackAxisTitleP[5] = "charge";
 297
 298   const char *cutlabel[6] = {"NoCut", "Muon", "Trigger", "eta", "ThetaAbs", "ChiS"};
 299   if(fIsMc)
 300   {
 301     // THn for tracking efficiency
 302     fHMuonParGen = new THnF("fHMuonParGen", "", 6, iTrackBin, 0, 0);
 303     for (Int_t i=0; i<6; i++)
 304     {
 305       fHMuonParGen->SetBinEdges(i, trackBins[i]);
 306       fHMuonParGen->GetAxis(i)->SetTitle(trackAxisTitle[i]);
 307     }
 308     fHMuonParGen->Sumw2();
 309     fOutputList->Add(fHMuonParGen);
 310
 311     fHMuonParGenSec = (THnF*) fHMuonParGen->Clone("fHMuonParGenSec");
 312     fHMuonParGenSec->Sumw2();
 313     fOutputList->Add(fHMuonParGenSec);
 314     if(fIsFPM)
 315     {
 316       fHMuonParGenFPM = (THnF*) fHMuonParGen->Clone("fHMuonParGenFPM");
 317       fHMuonParGenFPM->Sumw2();
 318       fOutputList->Add(fHMuonParGenFPM);
 319
 320       fHMuonParGenSecFPM = (THnF*) fHMuonParGen->Clone("fHMuonParGenSecFPM");
 321       fHMuonParGenSecFPM->Sumw2();
 322       fOutputList->Add(fHMuonParGenSecFPM);
 323     }
 324     if(fIsCutStudy)
 325     {
 326       for(Int_t i=0; i<6; i++)
 327       {
 328         fHMuonDetGen[i] = (THnF*)fHMuonParGen->Clone(Form("fHMuonDetGen_%s",cutlabel[i]));
 329         fHMuonDetGen[i]->Sumw2();
 330         fOutputList->Add(fHMuonDetGen[i]);
 331
 332         fHMuonDetRec[i] = (THnF*) fHMuonParGen->Clone(Form("fHMuonDetRec_%s",cutlabel[i]));
 333         fHMuonDetRec[i]->Sumw2();
 334         fOutputList->Add(fHMuonDetRec[i]);
 335
 336         fHHadDetRec[i] = (THnF*) fHMuonParGen->Clone(Form("fHHadDetRec_%s",cutlabel[i]));
 337         fHHadDetRec[i]->Sumw2();
 338         fOutputList->Add(fHHadDetRec[i]);
 339
 340         fHSecDetRec[i] = (THnF*) fHMuonParGen->Clone(Form("fHSecDetRec_%s",cutlabel[i]));
 341         fHSecDetRec[i]->Sumw2();
 342         fOutputList->Add(fHSecDetRec[i]);
 343       }
 344     }
 345     else
 346     {
 347       fHMuonDetGen[0] = (THnF*)fHMuonParGen->Clone(Form("fHMuonDetGen"));
 348       fHMuonDetGen[0]->Sumw2();
 349       fOutputList->Add(fHMuonDetGen[0]);
 350
 351       fHMuonDetRec[0] = (THnF*) fHMuonParGen->Clone(Form("fHMuonDetRec"));
 352       fHMuonDetRec[0]->Sumw2();
 353       fOutputList->Add(fHMuonDetRec[0]);
 354
 355       fHHadDetRec[0] = (THnF*) fHMuonParGen->Clone(Form("fHHadDetRec"));
 356       fHHadDetRec[0]->Sumw2();
 357       fOutputList->Add(fHHadDetRec[0]);
 358
 359       fHSecDetRec[0] = (THnF*) fHMuonParGen->Clone(Form("fHSecDetRec"));
 360       fHSecDetRec[0]->Sumw2();
 361       fOutputList->Add(fHSecDetRec[0]);
 362     }
 363
 364     const char *vertexlabel[4] = {"0_10", "10_90", "90_503", "503"};
 365     for(Int_t i=0; i<4; i++)
 366     {
 367       fHMuonParGenV[i] = (THnF*)fHMuonParGen->Clone(Form("fHMuonParGenV_%s",vertexlabel[i]));
 368       fHMuonParGenV[i]->Sumw2();
 369       fOutputList->Add(fHMuonParGenV[i]);
 370
 371       fHMuonParGenSecV[i] = (THnF*)fHMuonParGen->Clone(Form("fHMuonParGenSecV_%s",vertexlabel[i]));
 372       fHMuonParGenSecV[i]->Sumw2();
 373       fOutputList->Add(fHMuonParGenSecV[i]);
 374
 375       fHMuonDetRecV[i] = (THnF*)fHMuonParGen->Clone(Form("fHMuonDetRecV_%s",vertexlabel[i]));
 376       fHMuonDetRecV[i]->Sumw2();
 377       fOutputList->Add(fHMuonDetRecV[i]);
 378     }
 379
 380     fHMuonParGenP = new THnF("fHMuonParGenP", "", 6, iTrackBinP, 0, 0);
 381     for (Int_t i=0; i<6; i++)
 382     {
 383       fHMuonParGenP->SetBinEdges(i, trackBinsP[i]);
 384       fHMuonParGenP->GetAxis(i)->SetTitle(trackAxisTitleP[i]);
 385     }
 386     fHMuonParGenP->Sumw2();
 387     fOutputList->Add(fHMuonParGenP);
 388
 389     fHMuonDetGenP = (THnF*) fHMuonParGenP->Clone("fHMuonDetGenP");
 390     fHMuonDetGenP->Sumw2();
 391     fOutputList->Add(fHMuonDetGenP);
 392
 393     fHMuonDetRecP = (THnF*) fHMuonParGenP->Clone("fHMuonDetRecP");
 394     fHMuonDetRecP->Sumw2();
 395     fOutputList->Add(fHMuonDetRecP);
 396
 397     const char* MotherSpecies[4] = {"Pion","Kaon","D", "Etc"};
 398     const char *MuPt[3] = {"0005","0520","2040"};
 399     if(fZvProcess)
 400     {
 401       for(Int_t i=0; i<4; i++)
 402       {
 403         fHMuZv[i] = new TH1F(Form("fHMuZv_%s",MotherSpecies[i]), ";Z_{V}(cm)", 600, -500, 100);
 404         fOutputList->Add(fHMuZv[i]);
 405         fHMuRelZv[i] = new TH1F(Form("fHMuRelZv_%s",MotherSpecies[i]), ";|Z_{V,muon}-Zvtx|(cm)", 500, 0, 500);
 406         fOutputList->Add(fHMuRelZv[i]);
 407       }
 408     }
 409     if(fMDProcess)
 410     {
 411       for(Int_t i=0; i<4; i++)
 412       {
 413         fHMuMotherGenPt[i] = new TH2F(Form("fHMuMotherGenPt_%s",MotherSpecies[i]),";p_{T,muon}^{gen} (GeV/c);p_{T,mother}^{Truth} (GeV/c);",500, 0, 50, 500, 0, 50);
 414         fOutputList->Add(fHMuMotherGenPt[i]);
 415         fHMuMotherRecPt[i] = new TH2F(Form("fHMuMotherRecPt_%s",MotherSpecies[i]),";p_{T,muon}^{rec} (GeV/c);p_{T,mother}^{Truth} (GeV/c);",500, 0, 50, 500, 0, 50);
 416         fOutputList->Add(fHMuMotherRecPt[i]);
 417         fHMuMotherGenPhi[i] = new TH2F(Form("fHMuMotherGenPhi_%s",MotherSpecies[i]),";#phi_{gen};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
 418         fOutputList->Add(fHMuMotherGenPhi[i]);
 419         fHMuMotherRecPhi[i] = new TH2F(Form("fHMuMotherRecPhi_%s",MotherSpecies[i]),";#phi_{rec};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
 420         fOutputList->Add(fHMuMotherRecPhi[i]);
 421         fHMuMotherGenEta[i] = new TH2F(Form("fHMuMotherGenEta_%s",MotherSpecies[i]),";#eta_{gen};mother #eta;",100, -5., -1., 100, -5., -1.);
 422         fOutputList->Add(fHMuMotherGenEta[i]);
 423         fHMuMotherRecEta[i] = new TH2F(Form("fHMuMotherRecEta_%s",MotherSpecies[i]),";#eta_{rec};mother #eta;",100, -5., -1., 100, -5., -1.);
 424         fOutputList->Add(fHMuMotherRecEta[i]);
 425         fHMuDCA[i] =  new TH1F(Form("fHMuDCA_%s",MotherSpecies[i]), ";DCA", 100, 0, 50);
 426         fOutputList->Add(fHMuDCA[i]);
 427
 428         fHMuMotherGenPtSec[i] = new TH2F(Form("fHMuMotherGenPtSec_%s",MotherSpecies[i]),";p_{T,muon}^{gen} (GeV/c);p_{T,mother}^{Truth} (GeV/c);",500, 0, 50, 500, 0, 50);
 429         fOutputList->Add(fHMuMotherGenPtSec[i]);
 430         fHMuMotherRecPtSec[i] = new TH2F(Form("fHMuMotherRecPtSec_%s",MotherSpecies[i]),";p_{T,muon}^{rec} (GeV/c);p_{T,mother}^{Truth} (GeV/c);",500, 0, 50, 500, 0, 50);
 431         fOutputList->Add(fHMuMotherRecPtSec[i]);
 432         fHMuMotherGenPhiSec[i] = new TH2F(Form("fHMuMotherGenPhiSec_%s",MotherSpecies[i]),";#phi_{gen};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
 433         fOutputList->Add(fHMuMotherGenPhiSec[i]);
 434         fHMuMotherRecPhiSec[i] = new TH2F(Form("fHMuMotherRecPhiSec_%s",MotherSpecies[i]),";#phi_{rec};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
 435         fOutputList->Add(fHMuMotherRecPhiSec[i]);
 436         fHMuMotherGenEtaSec[i] = new TH2F(Form("fHMuMotherGenEtaSec_%s",MotherSpecies[i]),";#eta_{gen};mother #eta;",100, -5., -1., 100, -5., -1.);
 437         fOutputList->Add(fHMuMotherGenEtaSec[i]);
 438         fHMuMotherRecEtaSec[i] = new TH2F(Form("fHMuMotherRecEtaSec_%s",MotherSpecies[i]),";#eta_{rec};mother #eta;",100, -5., -1., 100, -5., -1.);
 439         fOutputList->Add(fHMuMotherRecEtaSec[i]);
 440         fHMuDCASec[i] =  new TH1F(Form("fHMuDCASec_%s",MotherSpecies[i]), ";DCA", 100, 0, 50);
 441         fOutputList->Add(fHMuDCASec[i]);
 442
 443         for(Int_t j=0; j<3; j++)
 444         {
 445           fHMuMohterPhiDifGen[i][j] = new TH1F(Form("fHMuMohterPhiDifGen_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#phi",100, -1.0*TMath::Pi(), TMath::Pi());
 446           fOutputList->Add(fHMuMohterPhiDifGen[i][j]);
 447           fHMuMohterPhiDifRec[i][j] = new TH1F(Form("fHMuMohterPhiDifRec_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#phi",100, -1.0*TMath::Pi(), TMath::Pi());
 448           fOutputList->Add(fHMuMohterPhiDifRec[i][j]);
 449           fHMuMohterEtaDifGen[i][j] = new TH1F(Form("fHMuMohterEtaDifGen_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#eta",100, -5.0, 5.0);
 450           fOutputList->Add(fHMuMohterEtaDifGen[i][j]);
 451           fHMuMohterEtaDifRec[i][j] = new TH1F(Form("fHMuMohterEtaDifRec_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#eta",100, -5.0, 5.0);
 452           fOutputList->Add(fHMuMohterEtaDifRec[i][j]);
 453
 454           fHMuMohterPhiDifGenSec[i][j] = new TH1F(Form("fHMuMohterPhiDifGenSec_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#phi",100, -1.0*TMath::Pi(), TMath::Pi());
 455           fOutputList->Add(fHMuMohterPhiDifGenSec[i][j]);
 456           fHMuMohterPhiDifRecSec[i][j] = new TH1F(Form("fHMuMohterPhiDifRecSec_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#phi",100, -1.0*TMath::Pi(), TMath::Pi());
 457           fOutputList->Add(fHMuMohterPhiDifRecSec[i][j]);
 458           fHMuMohterEtaDifGenSec[i][j] = new TH1F(Form("fHMuMohterEtaDifGenSec_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#eta",100, -5.0, 5.0);
 459           fOutputList->Add(fHMuMohterEtaDifGenSec[i][j]);
 460           fHMuMohterEtaDifRecSec[i][j] = new TH1F(Form("fHMuMohterEtaDifRecSec_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#eta",100, -5.0, 5.0);
 461           fOutputList->Add(fHMuMohterEtaDifRecSec[i][j]);
 462         }
 463       }
 464       fHZvRvPrim = new TH2F("fHZvRvPrim", "", 300, -500, 100, 200, 0, 800);
 465       fOutputList->Add(fHZvRvPrim);
 466
 467       fHZvRvSec = new TH2F("fHZvRvSec", "", 300, -500, 100, 200, 0, 800);
 468       fOutputList->Add(fHZvRvSec);
 469
 470       fHXvYvPrim = new TH2F("fHXvYvPrim", "", 200, -500, 500, 200, -500, 500);
 471       fOutputList->Add(fHXvYvPrim);
 472
 473       fHXvYvSec = new TH2F("fHXvYvSec", "", 200, -500, 500, 200, -500, 500);
 474       fOutputList->Add(fHXvYvSec);
 475     }
 476     if(fFeynmanX)
 477     {
 478       fHFXu = new TH1F("fHFXu",";x_{F}",1000, 0.0, 1.0);
 479       fOutputList->Add(fHFXu);
 480
 481       fHFXantiu = new TH1F("fHFXantiu",";x_{F}",1000, 0.0, 1.0);
 482       fOutputList->Add(fHFXantiu);
 483
 484       fHFXd = new TH1F("fHFXd",";x_{F}",1000, 0.0, 1.0);
 485       fOutputList->Add(fHFXd);
 486
 487       fHFXantid = new TH1F("fHFXantid",";x_{F}",1000, 0.0, 1.0);
 488       fOutputList->Add(fHFXantid);
 489
 490       fHFXg = new TH1F("fHFXg",";x_{F}",1000, 0.0, 1.0);
 491       fOutputList->Add(fHFXg);
 492
 493       fHFXetc = new TH1F("fHFXetc",";x_{F}",1000, 0.0, 1.0);
 494       fOutputList->Add(fHFXetc);
 495
 496       const char* muonptbin[2] = {"02","2"};
 497       const char *muonetabin[3] = {"4035","3530","3025"};
 498       for(Int_t i=0; i<2; i++)
 499       {
 500         for(Int_t j=0; j<3; j++)
 501         {
 502           fHFXmuonP[i][j] = new TH1F(Form("fHFXmuonP_%s_%s",muonptbin[i], muonetabin[j]),";x_{F}",1000, 0.0, 1.0);
 503           fOutputList->Add(fHFXmuonP[i][j]);
 504           fHFXmuonM[i][j] = new TH1F(Form("fHFXmuonM_%s_%s",muonptbin[i], muonetabin[j]),";x_{F}",1000, 0.0, 1.0);
 505           fOutputList->Add(fHFXmuonM[i][j]);
 506         }
 507       }
 508     }
 509     if(fScatFX)
 510     {
 511       fHaFx = new TH1F("fHaFx",";x_{F}",1000, 0.0, 1.0);
 512       fOutputList->Add(fHaFx);
 513       fHbFx = new TH1F("fHbFx",";x_{F}",1000, 0.0, 1.0);
 514       fOutputList->Add(fHbFx);
 515       fHabFxRatio = new TH1F("fHabFxRatio",";x_{F,a}/x_{F,b}",1000, 0.0, 5.0);
 516       fOutputList->Add(fHabFxRatio);
 517       fHabDeltaFx = new TH1F("fHabDeltaFx",";#Deltax_{F}",1000, -2.0, 2.0);
 518       fOutputList->Add(fHabDeltaFx);
 519       fHabRelDeltaFx = new TH1F("fHabRelDeltaFx",";(x_{a}-x_{b})/(x_{a}+x{b})",1000, -1.0, 1.0);
 520       fOutputList->Add(fHabRelDeltaFx);
 521
 522       const char* MuonPt[3] = {"Mu01","Mu12", "Mu2"};
 523       for(Int_t i=0; i<3; i++)
 524       {
 525         Int_t NabBin = 22;
 526         Double_t abBins[22+1] = { 0.0, 0.000001, 0.000002, 0.000004, 0.00001, 0.00002, 0.00004, 0.0001, 0.0002, 0.0004, 0.001, 0.002, 0.004, 0.01, 0.02, 0.04, 0.1, 0.2, 0.3, 0.4, 0.5, 0.7, 1};
 527
 528         fHabFxMu[i] = new TH2F(Form("fHabFxMu_%s",MuonPt[i]),";x_{F,a};x_{F,b}",NabBin, abBins, NabBin, abBins);
 529         fOutputList->Add(fHabFxMu[i]);
 530
 531         fHabFxRatioMu[i] = new TH1F(Form("fHabFxRatioMu_%s",MuonPt[i]),";x_{F,a}/x_{F,b}",1000, 0.0, 5.0);
 532         fOutputList->Add(fHabFxRatioMu[i]);
 533
 534         fHabDeltaFxMu[i] = new TH1F(Form("fHabDeltaFxMu_%s",MuonPt[i]),";#Deltax_{F}",1000, -2.0, 2.0);
 535         fOutputList->Add(fHabDeltaFxMu[i]);
 536
 537         fHabRelDeltaFxMu[i] = new TH1F(Form("fHabRelDeltaFxMu_%s",MuonPt[i]),";(x_{a}-x_{b})/(x_{a}+x{b})",1000, -1.0, 1.0);
 538         fOutputList->Add(fHabRelDeltaFxMu[i]);
 539       }
 540     }
 541   }
 542   else
 543   {
 544     fHMuonDetRec[0] = new THnF(Form("fHMuonDetRec_%s", cutlabel[0]), "", 6, iTrackBin, 0, 0);
 545     for (Int_t i=0; i<6; i++)
 546     {
 547       fHMuonDetRec[0]->SetBinEdges(i, trackBins[i]);
 548       fHMuonDetRec[0]->GetAxis(i)->SetTitle(trackAxisTitle[i]);
 549     }
 550     fHMuonDetRec[0]->Sumw2();
 551     fOutputList->Add(fHMuonDetRec[0]);
 552
 553     for(Int_t i=1; i<6; i++)
 554     {
 555
 556       fHMuonDetRec[i] = (THnF*) fHMuonDetRec[0]->Clone(Form("fHMuonDetRec_%s",cutlabel[i]));
 557       fHMuonDetRec[i]->Sumw2();
 558       fOutputList->Add(fHMuonDetRec[i]);
 559     }
 560
 561    fHMuonDetRecP = new THnF("fHMuonDetRecP", "", 6, iTrackBinP, 0, 0);
 562     for (Int_t i=0; i<6; i++)
 563     {
 564       fHMuonDetRecP->SetBinEdges(i, trackBinsP[i]);
 565       fHMuonDetRecP->GetAxis(i)->SetTitle(trackAxisTitleP[i]);
 566     }
 567     fHMuonDetRecP->Sumw2();
 568     fOutputList->Add(fHMuonDetRecP);
 569   }
 570   PostData(1, fOutputList);
 571 }
 572
 573 //________________________________________________________________________
 574 Bool_t AliMuonEffMC::Notify()
 575 {
 576   // Implemented Notify() to read the cross sections
 577   // and number of trials from pyxsec.root
 578   if(fIsPythia)
 579   {
 580     fHEventStat->Fill(2.5);
 581
 582     TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 583     Float_t xsection = 0;
 584     Float_t ftrials  = 1;
 585
 586     if(tree){
 587       TFile *curfile = tree->GetCurrentFile();
 588       if (!curfile) {
 589         Error("Notify","No current file");
 590         return kFALSE;
 591       }
 592
 593       AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 594       fHXsec->Fill("<#sigma>",xsection);
 595       fHTrials->Fill("#sum{ntrials}",ftrials);
 596     }
 597   }
 598   return kTRUE;
 599 }
 600 //________________________________________________________________________
 601 void AliMuonEffMC::UserExec(Option_t *)
 602 {
 603   // Main loop, Called for each event
 604   Int_t ntrks = 0; // number of tracks in an event
 605
 606   if(((TString)InputEvent()->IsA()->GetName())=="AliAODEvent")
 607   {
 608     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 609     if (!fAOD) { AliError("AOD event not found. Nothing done!"); return; }
 610     ntrks = fAOD->GetNTracks();
 611   }
 612   else
 613   {
 614     fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 615     if (!fESD) { AliError("ESD event not found. Nothing done!"); return; }
 616     ntrks = fESD->GetNumberOfMuonTracks();
 617   }
 618
 619   fHEventStat->Fill(0.5);
 620   if(fIsMc)
 621   {
 622     fMC = MCEvent();
 623     if (!fMC) { AliError("MC event not avaliable."); return; }
 624     fStack = fMC->Stack();
 625   }
 626
 627   // Centrality, vertex, other event variables...
 628   if(fAOD)
 629   {
 630     const AliAODVertex* vertex = fAOD->GetPrimaryVertex();
 631     fZVertex = vertex->GetZ();
 632     if(fAOD->GetCentrality())  fCentrality = fAOD->GetCentrality()->GetCentralityPercentile(fCentralityEstimator);
 633   }
 634   else if(fESD)
 635   {
 636     const AliESDVertex* vertex = fESD->GetPrimaryVertex();
 637     fZVertex = vertex->GetZ();
 638     if(fESD->GetCentrality()) fCentrality = fESD->GetCentrality()->GetCentralityPercentile(fCentralityEstimator);
 639   }
 640
 641   if ((fESD && !VertexOk(fESD)) || (fAOD && !VertexOk(fAOD))) { //AliInfo(Form("Event REJECTED. z = %.1f", fZVertex));
 642     return; }
 643   if (fCentrality > 100. || fCentrality < -1.5) { //AliInfo(Form("Event REJECTED. fCentrality = %.1f", fCentrality));
 644     return; }
 645
 646   if(fCentrality < 0) fCentrality = 0.5; //ad hoc centrality for pp
 647   // Fill Event histogram
 648   fHEvt->Fill(fZVertex, fCentrality);
 649   fHEventStat->Fill(1.5);
 650
 651   ULong64_t trigword = 0;
 652   if(fAOD) trigword=fAOD->GetTriggerMask();
 653   else if(fESD) trigword=fESD->GetTriggerMask();
 654
 655   if (trigword & 0x01) fHEventStat->Fill(17.5);
 656   if (trigword & 0x02) fHEventStat->Fill(3.5);
 657   if (trigword & 0x04) fHEventStat->Fill(4.5);
 658   if (trigword & 0x08) fHEventStat->Fill(5.5);
 659   if (trigword & 0x010) fHEventStat->Fill(6.5);
 660   if (trigword & 0x020) fHEventStat->Fill(7.5);
 661   if (trigword & 0x040) fHEventStat->Fill(8.5);
 662   if (trigword & 0x080) fHEventStat->Fill(9.5);
 663   if (trigword & 0x100) fHEventStat->Fill(10.5);
 664   if (trigword & 0x200) fHEventStat->Fill(11.5);
 665   if (trigword & 0x400) fHEventStat->Fill(12.5);
 666   if (trigword & 0x800) fHEventStat->Fill(13.5);
 667   if (trigword & 0x1000) fHEventStat->Fill(14.5);
 668   if (trigword & 0x2000) fHEventStat->Fill(15.5);
 669   if (trigword & 0x4000) fHEventStat->Fill(16.5);
 670
 671   if(fIsMc)
 672   {
 673     // generated level loop
 674     for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)
 675     {
 676       Int_t FirstPrimMother = 0;
 677       Double_t FirstPrimZvtx = 0.0;
 678       Int_t Zvtxbin = 0;
 679       if(fAOD)
 680       {
 681         AliAODMCParticle *AodMcParticle  = (AliAODMCParticle*)fMC->GetTrack(ipart);
 682         if(AodMcParticle->Eta() < -4.0 || AodMcParticle->Eta() > -2.5) continue;
 683
 684         if(TMath::Abs(AodMcParticle->PdgCode())==13)
 685         {
 686           FirstPrimMother = GetFirstPrimaryMother(ipart);
 687           FirstPrimZvtx = GetFirstPrimaryVertex(ipart);
 688           Zvtxbin = GetZVertexBin(FirstPrimZvtx);
 689
 690           Double_t fillArrayParGen[6] = { AodMcParticle->Eta(), AodMcParticle->Pt(), fCentrality, fZVertex, AodMcParticle->Phi(), AodMcParticle->Charge() };
 691           if(AodMcParticle->IsPrimary())
 692           {
 693             fHMuonParGen->Fill(fillArrayParGen);
 694             fHMuonParGenV[Zvtxbin]->Fill(fillArrayParGen);
 695           }
 696           else
 697           {
 698             fHMuonParGenSec->Fill(fillArrayParGen);
 699             fHMuonParGenSecV[Zvtxbin]->Fill(fillArrayParGen);
 700           }
 701           Double_t fillArrayParGenP[6] = { AodMcParticle->Eta(), AodMcParticle->P(), fCentrality, fZVertex, AodMcParticle->Phi(), AodMcParticle->Charge() };
 702           if(AodMcParticle->IsPrimary()) fHMuonParGenP->Fill(fillArrayParGenP);
 703         }
 704       }
 705       else if(fESD)
 706       {
 707         AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(ipart);
 708         if(McParticle->Eta() < -4.0 || McParticle->Eta() > -2.5) continue;
 709
 710         if(TMath::Abs(McParticle->PdgCode())==13)
 711         {
 712           FirstPrimMother = GetFirstPrimaryMother(ipart);
 713           FirstPrimZvtx = GetFirstPrimaryVertex(ipart);
 714           Zvtxbin = GetZVertexBin(FirstPrimZvtx);
 715
 716           Double_t fillArrayParGen[6] = { McParticle->Eta(), McParticle->Pt(), fCentrality, fZVertex, McParticle->Phi(), McParticle->Charge() };
 717           AliMCParticle *FPrimaryMother  = (AliMCParticle*)fMC->GetTrack(FirstPrimMother);
 718           Double_t fillArrayParGenFPM[6] = { FPrimaryMother->Eta(), FPrimaryMother->Pt(), fCentrality, fZVertex, FPrimaryMother->Phi(), FPrimaryMother->Charge() };
 719
 720           if(McParticle->GetMother()< fStack->GetNprimary())
 721           {
 722             fHMuonParGen->Fill(fillArrayParGen);
 723             fHMuonParGenV[Zvtxbin]->Fill(fillArrayParGen);
 724             if(fIsFPM) fHMuonParGenFPM->Fill(fillArrayParGenFPM);
 725           }
 726           else
 727           {
 728             fHMuonParGenSec->Fill(fillArrayParGen);
 729             fHMuonParGenSecV[Zvtxbin]->Fill(fillArrayParGen);
 730             if(fIsFPM) fHMuonParGenSecFPM->Fill(fillArrayParGenFPM);
 731           }
 732           Double_t fillArrayParGenP[6] = { McParticle->Eta(), McParticle->P(), fCentrality, fZVertex, McParticle->Phi(), McParticle->Charge() };
 733           if(McParticle->GetMother()< fStack->GetNprimary()) fHMuonParGenP->Fill(fillArrayParGenP);
 734         }
 735       }
 736     }
 737   }
 738
 739   // reconstructed level loop
 740   for (Int_t iTrack = 0; iTrack<ntrks; iTrack++)
 741   {
 742     Int_t label = 0;
 743     // reconstructed track variables
 744     Double_t trackpt = 0;
 745     Double_t tracketa = 0;
 746     Double_t trackphi = 0;
 747     Double_t trackcharge = 0;
 748     Double_t trackp = 0;
 749     Double_t dcavalue = 0;
 750     Int_t cutNum = 0;
 751     Bool_t isprimary = kFALSE;
 752     // reconstructed track's matched truth particle variables
 753     Double_t mcpt = 0;
 754     Double_t mceta = 0;
 755     Double_t mcphi = 0;
 756     Double_t mcp = 0;
 757     Double_t mccharge = 0;
 758     Double_t mcZv = 0;
 759     // first primary mother variables
 760     Int_t motherlabel =0;
 761     Int_t motherpdg = 0;
 762     Double_t motherpt = 0;
 763     Double_t mothereta = 0;
 764     Double_t motherphi = 0;
 765     Double_t motherXv = 0;
 766     Double_t motherYv = 0;
 767     Double_t motherZv = 0;
 768
 769     if(fAOD)
 770     {
 771       AliAODTrack* muonTrack = (AliAODTrack*)fAOD->GetTrack(iTrack);
 772       if(muonTrack)
 773       {
 774         if(!(muonTrack->IsMuonTrack())) continue;
 775         cutNum = IsGoodMUONtrack(*muonTrack);
 776         trackpt = muonTrack->Pt();
 777         tracketa = muonTrack->Eta();
 778         trackphi = muonTrack->Phi();
 779         trackcharge = muonTrack->Charge();
 780         trackp = muonTrack->P();
 781         label =  TMath::Abs(muonTrack->GetLabel());
 782         if (label>=fMC->GetNumberOfTracks()) {
 783           AliError(Form("Label %d larger than number of particles on stack %d\n",label,fMC->GetNumberOfTracks()));
 784           continue;
 785         }
 786       }
 787     }
 788     else if(fESD)
 789     {
 790       AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iTrack);
 791       if(muonTrack)
 792       {
 793         cutNum = IsGoodMUONtrack(*muonTrack);
 794         trackpt = muonTrack->Pt();
 795         tracketa = muonTrack->Eta();
 796         trackphi = muonTrack->Phi();
 797         trackcharge = muonTrack->Charge();
 798         trackp = muonTrack->P();
 799         label =  TMath::Abs(muonTrack->GetLabel());
 800         dcavalue = muonTrack->GetDCA();
 801       }
 802     }
 803     Double_t fillArrayDetRec[6] = { tracketa, trackpt, fCentrality, fZVertex, trackphi, trackcharge };
 804     Double_t fillArrayDetRecP[6] = { tracketa, trackp, fCentrality, fZVertex, trackphi, trackcharge };
 805     if(fIsCutStudy) { for(Int_t icut = 0; icut <= cutNum; icut++) { fHMuonDetRec[icut]->Fill(fillArrayDetRec); }}
 806     else { if(cutNum==5) fHMuonDetRec[0]->Fill(fillArrayDetRec); }
 807     if(cutNum==5) fHMuonDetRecP->Fill(fillArrayDetRecP);
 808
 809     if(fIsMc)
 810     {
 811       if(fAOD)
 812       {
 813         AliAODMCParticle *aodMcParticle  = (AliAODMCParticle*)fMC->GetTrack(label);
 814         isprimary = (aodMcParticle->IsPrimary()) ? kTRUE : kFALSE;
 815         if(TMath::Abs(aodMcParticle->PdgCode()) != 13)
 816         {
 817           if(fIsCutStudy){ for(Int_t icut=0; icut<=cutNum; icut++) fHHadDetRec[icut]->Fill(fillArrayDetRec); }
 818           else { if(cutNum==5) fHHadDetRec[0]->Fill(fillArrayDetRec); }
 819         }
 820         else if(!isprimary)
 821         {
 822           if(fIsCutStudy){ for(Int_t icut=0; icut<=cutNum; icut++) fHSecDetRec[icut]->Fill(fillArrayDetRec); }
 823           else { if(cutNum==5) fHSecDetRec[0]->Fill(fillArrayDetRec); }
 824         }
 825         mcpt = aodMcParticle->Pt();
 826         mceta = aodMcParticle->Eta();
 827         mcphi = aodMcParticle->Phi();
 828         mccharge = aodMcParticle->Charge();
 829         mcp = aodMcParticle->P();
 830         mcZv = aodMcParticle->Zv();
 831         motherlabel = aodMcParticle->GetMother();
 832       }
 833       else if(fESD)
 834       {
 835         AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(label);
 836         isprimary = (McParticle->GetMother()<fStack->GetNprimary()) ? kTRUE : kFALSE;
 837         if(TMath::Abs(McParticle->PdgCode())!=13)
 838         {
 839           if(fIsCutStudy) { for(Int_t icut=0; icut<=cutNum; icut++) fHHadDetRec[icut]->Fill(fillArrayDetRec); }
 840           else{ if(cutNum==5) fHHadDetRec[0]->Fill(fillArrayDetRec); }
 841         }
 842         if(!isprimary)
 843         {
 844           if(fIsCutStudy) { for(Int_t icut=0; icut<=cutNum; icut++) fHSecDetRec[icut]->Fill(fillArrayDetRec); }
 845           else{ if(cutNum==5)  fHSecDetRec[0]->Fill(fillArrayDetRec); }
 846         }
 847         mcpt = McParticle->Pt();
 848         mceta = McParticle->Eta();
 849         mcphi = McParticle->Phi();
 850         mccharge = McParticle->Charge();
 851         mcp = McParticle->P();
 852         mcZv = McParticle->Zv();
 853
 854         motherlabel = GetFirstPrimaryMother(label);
 855         Double_t primzvtx = GetFirstPrimaryVertex(label);
 856         Int_t reczvtxbin = GetZVertexBin(primzvtx);
 857         if(cutNum==5) fHMuonDetRecV[reczvtxbin]->Fill(fillArrayDetRec);
 858
 859         if(motherlabel > -1)
 860         {
 861           AliMCParticle *FPrimaryMother  = (AliMCParticle*)fMC->GetTrack(motherlabel);
 862           motherpdg = TMath::Abs(FPrimaryMother->PdgCode());
 863           motherpt = FPrimaryMother->Pt();
 864           mothereta = FPrimaryMother->Eta();
 865           motherphi = FPrimaryMother->Phi();
 866           AliMCParticle *DaughtParticle  = (AliMCParticle*)fMC->GetTrack(FPrimaryMother->GetFirstDaughter());
 867           motherXv = DaughtParticle->Xv();
 868           motherYv = DaughtParticle->Yv();
 869           motherZv = DaughtParticle->Zv();
 870         }
 871       }
 872       Double_t fillArrayDetGen[6] = { mceta, mcpt, fCentrality, fZVertex, mcphi, mccharge };
 873       if(fIsCutStudy){ for(Int_t icut=0; icut <= cutNum; icut++) fHMuonDetGen[icut]->Fill(fillArrayDetGen); }
 874       else {if(cutNum==5) fHMuonDetGen[0]->Fill(fillArrayDetGen); }
 875
 876       Double_t fillArrayDetGenP[6] = { mceta, mcp, fCentrality, fZVertex, mcphi, mccharge };
 877       if(cutNum==5)fHMuonDetGenP->Fill(fillArrayDetGenP);
 878
 879       Int_t motherbin = GetMotherBin(motherpdg);
 880       // muon Z-Vertex process
 881       if(fZvProcess && cutNum==5)
 882       {
 883         fHMuZv[motherbin]->Fill(mcZv);
 884         fHMuRelZv[motherbin]->Fill(TMath::Abs(mcZv-fZVertex));
 885       }
 886       // mother-daughter kinematic relation
 887       if(fMDProcess && cutNum==5 && motherlabel>0)
 888       {
 889         if(isprimary)
 890         {
 891           fHZvRvPrim->Fill(motherZv, TMath::Sqrt(motherXv*motherXv + motherYv*motherYv));
 892           fHXvYvPrim->Fill(motherXv, motherYv);
 893         }
 894         else
 895         {
 896           fHZvRvSec->Fill(motherZv, TMath::Sqrt(motherXv*motherXv + motherYv*motherYv));
 897           fHXvYvSec->Fill(motherXv, motherYv);
 898         }
 899         MDProcess(isprimary, motherpdg, mcpt, mcphi, mceta, trackpt, trackphi, tracketa, motherpt, motherphi, mothereta, dcavalue);
 900       }
 901     }
 902   }
 903
 904   if(fIsMc)
 905   {
 906     if(fFeynmanX) FeynmanX();
 907     if(fScatFX) ScatFX();
 908   }
 909
 910   PostData(1, fOutputList);
 911   return;
 912 }
 913
 914 //________________________________________________________________________
 915 void AliMuonEffMC::Terminate(Option_t *)
 916 {
 917   // Draw result to the screen
 918   // Called once at the end of the query
 919 }
 920 //________________________________________________________________________
 921 Bool_t AliMuonEffMC::VertexOk(TObject* obj) const
 922 {
 923   // Modified from AliAnalyseLeadingTrackUE::VertexSelection()
 924
 925   Int_t nContributors  = 0;
 926   Double_t zVertex     = 999;
 927   TString name("");
 928
 929   if (obj->InheritsFrom("AliESDEvent")) {
 930     AliESDEvent* esdevt = (AliESDEvent*) obj;
 931     const AliESDVertex* vtx = esdevt->GetPrimaryVertex();
 932     if (!vtx)
 933       return 0;
 934     nContributors = vtx->GetNContributors();
 935     zVertex       = vtx->GetZ();
 936     name          = vtx->GetName();
 937   }
 938   else if (obj->InheritsFrom("AliAODEvent")) {
 939     AliAODEvent* aodevt = (AliAODEvent*) obj;
 940     if (aodevt->GetNumberOfVertices() < 1)
 941       return 0;
 942     const AliAODVertex* vtx = aodevt->GetPrimaryVertex();
 943     nContributors = vtx->GetNContributors();
 944     zVertex       = vtx->GetZ();
 945     name          = vtx->GetName();
 946   }
 947
 948   // Reject if TPC-only vertex
 949   if (name.CompareTo("TPCVertex")==0)
 950     return kFALSE;
 951
 952   // Check # contributors and range...
 953   if( nContributors < 1 || TMath::Abs(zVertex) > 10 ) {
 954     return kFALSE;
 955   }
 956
 957   return kTRUE;
 958 }
 959
 960 //________________________________________________________________________
 961 Int_t AliMuonEffMC::IsGoodMUONtrack(AliESDMuonTrack &track)
 962 {
 963   // Applying track cuts for MUON tracks
 964   Int_t cutnum = 0;
 965   if(track.ContainTrackerData()) cutnum = 1;
 966   else return cutnum;
 967
 968   if(track.GetMatchTrigger() > 0) cutnum = 2;
 969   else return cutnum;
 970
 971   Double_t thetaTrackAbsEnd = TMath::ATan(track.GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
 972   Double_t eta = track.Eta();
 973
 974   // Eta cut
 975   if(eta > -4. && eta < -2.5) cutnum = 3;
 976   else return cutnum;
 977
 978   // Theta cut at absorber end
 979   if(thetaTrackAbsEnd > 2. && thetaTrackAbsEnd < 10.) cutnum = 4;
 980   else return cutnum;
 981
 982   // Chi-square cut
 983   if(track.GetNormalizedChi2() < fChiSquareNormCut) cutnum = 5;
 984   else return cutnum;
 985
 986   return cutnum;
 987 }
 988
 989 //________________________________________________________________________
 990 Int_t AliMuonEffMC::IsGoodMUONtrack(AliAODTrack &track)
 991 {
 992   Int_t cutnum = 0;
 993   if(track.IsMuonTrack()) cutnum = 1;
 994   else return cutnum;
 995
 996   if (track.GetMatchTrigger() > 0) cutnum = 2;
 997   else return cutnum;
 998
 999   Double_t thetaTrackAbsEnd = TMath::ATan(track.GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
1000   Double_t eta = track.Eta();
1001
1002   // Eta cut
1003   if(eta > -4. && eta < -2.5) cutnum = 3;
1004   else return cutnum;
1005
1006   // Theta cut at absorber end
1007   if(thetaTrackAbsEnd > 2. && thetaTrackAbsEnd < 10.) cutnum = 4;
1008   else return cutnum;
1009
1010   // Chi-square cut
1011   if(track.Chi2perNDF() < fChiSquareNormCut) cutnum = 5;
1012   else return cutnum;
1013
1014   return cutnum;
1015 }
1016
1017 //________________________________________________________________________
1018 void AliMuonEffMC::MDProcess(Bool_t isprimary, Int_t motherpdg, Double_t mcpt, Double_t mcphi, Double_t mceta, Double_t trackpt, Double_t trackphi, Double_t tracketa, Double_t motherpt, Double_t motherphi, Double_t mothereta, Double_t dcavalue)
1019 {
1020   Int_t motherbin = GetMotherBin(motherpdg);
1021   Int_t genptbin = -1;
1022   Int_t recptbin = -1;
1023
1024   if((0. <= mcpt) && (mcpt < 0.5)) genptbin = 0;
1025   else if((0.5 <= mcpt) && (mcpt < 2.0)) genptbin = 1;
1026   else genptbin = 2;
1027
1028   if((0. <= trackpt) && (trackpt < 0.5)) recptbin = 0;
1029   else if((0.5 <= trackpt) && (trackpt < 2.0)) recptbin = 1;
1030   else recptbin = 2;
1031   if(isprimary)
1032   {
1033     fHMuMotherGenPt[motherbin]->Fill(mcpt, motherpt);
1034     fHMuMotherRecPt[motherbin]->Fill(trackpt, motherpt);
1035     fHMuMotherGenPhi[motherbin]->Fill(mcphi, motherphi);
1036     fHMuMotherRecPhi[motherbin]->Fill(trackphi, motherphi);
1037     fHMuMotherGenEta[motherbin]->Fill(mceta, mothereta);
1038     fHMuMotherRecEta[motherbin]->Fill(tracketa, mothereta);
1039
1040     // generated
1041     fHMuMohterPhiDifGen[motherbin][genptbin]->Fill(deltaphi(motherphi-mcphi));
1042     fHMuMohterEtaDifGen[motherbin][genptbin]->Fill(mothereta-mceta);
1043
1044     // reconstructed
1045     fHMuMohterPhiDifRec[motherbin][recptbin]->Fill(deltaphi(motherphi-trackphi));
1046     fHMuMohterEtaDifRec[motherbin][recptbin]->Fill(mothereta-tracketa);
1047
1048     // DCA
1049     if(fESD) fHMuDCA[motherbin]->Fill(dcavalue);
1050   }
1051   else
1052   {
1053     fHMuMotherGenPtSec[motherbin]->Fill(mcpt, motherpt);
1054     fHMuMotherRecPtSec[motherbin]->Fill(trackpt, motherpt);
1055     fHMuMotherGenPhiSec[motherbin]->Fill(mcphi, motherphi);
1056     fHMuMotherRecPhiSec[motherbin]->Fill(trackphi, motherphi);
1057     fHMuMotherGenEtaSec[motherbin]->Fill(mceta, mothereta);
1058     fHMuMotherRecEtaSec[motherbin]->Fill(tracketa, mothereta);
1059
1060     // generated
1061     fHMuMohterPhiDifGenSec[motherbin][genptbin]->Fill(deltaphi(motherphi-mcphi));
1062     fHMuMohterEtaDifGenSec[motherbin][genptbin]->Fill(mothereta-mceta);
1063
1064     // reconstructed
1065     fHMuMohterPhiDifRecSec[motherbin][recptbin]->Fill(deltaphi(motherphi-trackphi));
1066     fHMuMohterEtaDifRecSec[motherbin][recptbin]->Fill(mothereta-tracketa);
1067
1068     // DCA
1069     if(fESD) fHMuDCASec[motherbin]->Fill(dcavalue);
1070   }
1071 }
1072
1073 //________________________________________________________________________
1074 Double_t AliMuonEffMC::deltaphi(Double_t phi)
1075 {
1076   if(phi < -1.0*TMath::Pi()) { return (phi + TMath::TwoPi()); }
1077   else if(phi > TMath::Pi()) { return (phi - TMath::TwoPi()); }
1078   else { return phi; }
1079 }
1080
1081 //________________________________________________________________________
1082 Int_t AliMuonEffMC::GetMotherBin(Int_t motherpdg)
1083 {
1084   Int_t motherbin = -1;
1085
1086   if(motherpdg==211) motherbin = 0;
1087   else if(motherpdg==321) motherbin = 1;
1088   else if(motherpdg==411 || motherpdg==413 || motherpdg==421 || motherpdg==423 || motherpdg==431 || motherpdg==433 || motherpdg==10413 || motherpdg==10411 || motherpdg==10423 || motherpdg==10421 || motherpdg==10433 || motherpdg==10431 || motherpdg==20413 || motherpdg==415 || motherpdg==20423 || motherpdg==425 || motherpdg==20433 || motherpdg==435) motherbin = 2;
1089   else motherbin = 3;
1090
1091   return motherbin;
1092 }
1093
1094 //________________________________________________________________________
1095 void AliMuonEffMC::FeynmanX()
1096 {
1097   for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)
1098   {
1099     TParticle *particle = fStack->Particle(ipart);
1100     if(particle->GetFirstMother()==0 || particle->GetFirstMother()==1)
1101     {
1102       Int_t pdgcode = particle->GetPdgCode();
1103       if(pdgcode==2) fHFXu->Fill(TMath::Abs(particle->Pz())/1380.0);
1104       else if(pdgcode==-2) fHFXantiu->Fill(TMath::Abs(particle->Pz())/1380.0);
1105       else if(pdgcode==1) fHFXd->Fill(TMath::Abs(particle->Pz())/1380.0);
1106       else if(pdgcode==-1) fHFXantid->Fill(TMath::Abs(particle->Pz())/1380.0);
1107       else if(TMath::Abs(pdgcode==21)) fHFXg->Fill(TMath::Abs(particle->Pz())/1380.0);
1108       else fHFXetc->Fill(TMath::Abs(particle->Pz())/1380.0);
1109     }
1110   }
1111
1112   if(fESD)
1113   {
1114     for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfMuonTracks(); iTrack++)
1115     {
1116       AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iTrack);
1117       if(muonTrack)
1118       {
1119         if(!IsGoodMUONtrack(*muonTrack)) continue;
1120         Int_t label = muonTrack->GetLabel();
1121         if(TMath::Abs(label) > fMC->GetNumberOfTracks()) continue;
1122
1123         TParticle *particle = fStack->Particle(TMath::Abs(label));
1124         Int_t motherlabel = particle->GetFirstMother();
1125         if(!(motherlabel==-1) && !(motherlabel==0) && !(motherlabel==1))
1126         {
1127           while(1)
1128           {
1129             particle = fStack->Particle(motherlabel);
1130             motherlabel = particle->GetFirstMother();
1131             if(motherlabel==-1 || motherlabel==0 || motherlabel==1) break;
1132           }
1133         }
1134         if(motherlabel==-1) continue;
1135
1136         Int_t ptbin = -1;
1137         Int_t etabin = -1;
1138
1139         if(muonTrack->Pt() < 2.0) ptbin=0;
1140         else ptbin=1;
1141
1142         if(muonTrack->Eta() <= -3.5) etabin=0;
1143         else if(-3.5<muonTrack->Eta() && muonTrack->Eta()<=-3.0) etabin=1;
1144         else etabin=2;
1145
1146         if(motherlabel==0) fHFXmuonP[ptbin][etabin]->Fill(TMath::Abs(particle->Pz())/1380.0);
1147         else if(motherlabel==1) fHFXmuonM[ptbin][etabin]->Fill(TMath::Abs(particle->Pz())/1380.0);
1148       }
1149     }
1150   }
1151   else if(fAOD)
1152   {
1153     for (Int_t iTrack = 0; iTrack<fAOD->GetNTracks(); iTrack++)
1154     {
1155       AliAODTrack* muonTrack = (AliAODTrack*)fAOD->GetTrack(iTrack);
1156       if(muonTrack)
1157       {
1158         if(!(IsGoodMUONtrack(*muonTrack)) || !(muonTrack->IsMuonTrack())) continue;
1159         Int_t label = muonTrack->GetLabel();
1160         if(TMath::Abs(label) > fMC->GetNumberOfTracks()) continue;
1161
1162         TParticle *particle = fStack->Particle(TMath::Abs(label));
1163         Int_t motherlabel = particle->GetFirstMother();
1164         if(!(motherlabel==-1) && !(motherlabel==0) && !(motherlabel==1))
1165         {
1166           while(1)
1167           {
1168             particle = fStack->Particle(motherlabel);
1169             motherlabel = particle->GetFirstMother();
1170             if(motherlabel==-1 || motherlabel==0 || motherlabel==1) break;
1171           }
1172         }
1173         if(motherlabel==-1) continue;
1174
1175         Int_t ptbin = -1;
1176         Int_t etabin = -1;
1177         if(muonTrack->Pt() < 2.0) ptbin=0;
1178         else ptbin=1;
1179
1180         if(muonTrack->Eta() <= -3.5) etabin=0;
1181         else if(-3.5<muonTrack->Eta() && muonTrack->Eta()<=-3.0) etabin=1;
1182         else etabin=2;
1183
1184         if(motherlabel==0) fHFXmuonP[ptbin][etabin]->Fill(TMath::Abs(particle->Pz())/1380.0);
1185         else if(motherlabel==1) fHFXmuonM[ptbin][etabin]->Fill(TMath::Abs(particle->Pz())/1380.0);
1186       }
1187     }
1188   }
1189 }
1190 //________________________________________________________________________
1191 void AliMuonEffMC::ScatFX()
1192 {
1193   TParticle *parta = fStack->Particle(6);
1194   TParticle *partb = fStack->Particle(7);
1195   if(parta->GetFirstMother()!=-1 || parta->GetFirstMother()!=-1) return;
1196
1197   Double_t xa = TMath::Abs(0.5*(parta->Energy() + partb->Energy() + parta->Pz() + partb->Pz())/1380.0);
1198   Double_t xb = TMath::Abs(0.5*(parta->Energy() + partb->Energy() - parta->Pz() - partb->Pz())/1380.0);
1199
1200   fHaFx->Fill(xa);
1201   fHbFx->Fill(xb);
1202   fHabFxRatio->Fill(xa/xb);
1203   fHabDeltaFx->Fill(xa-xb);
1204   fHabRelDeltaFx->Fill((xa-xb)/(xa+xb));
1205
1206   if(fESD)
1207   {
1208     for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfMuonTracks(); iTrack++)
1209     {
1210       Int_t MuonPtBin = 0;
1211       AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iTrack);
1212       if(muonTrack)
1213       {
1214         if(!IsGoodMUONtrack(*muonTrack)) continue;
1215
1216         if(muonTrack->Pt() < 1) MuonPtBin = 0;
1217         else if(muonTrack->Pt() < 2) MuonPtBin = 1;
1218         else MuonPtBin = 2;
1219
1220         fHabFxMu[MuonPtBin]->Fill(xa, xb);
1221         fHabFxRatioMu[MuonPtBin]->Fill(xa/xb);
1222         fHabDeltaFxMu[MuonPtBin]->Fill(xa-xb);
1223         fHabRelDeltaFxMu[MuonPtBin]->Fill((xa-xb)/(xa+xb));
1224       }
1225     }
1226   }
1227 }
1228 //________________________________________________________________________
1229 Int_t AliMuonEffMC::GetZVertexBin(Double_t zvertex)
1230 {
1231   if(TMath::Abs(zvertex) < 10.0) { return 0; }
1232   else if(TMath::Abs(zvertex) < 90.0) { return 1; }
1233   else if(TMath::Abs(zvertex) < 503.0) { return 2; }
1234   else { return 3; }
1235 }
1236
1237 //________________________________________________________________________
1238 Int_t AliMuonEffMC::GetFirstPrimaryMother(Int_t muonlabel)
1239 {
1240   if(fAOD) return 1;
1241   else if(fESD)
1242   {
1243     AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(muonlabel);
1244     if(McParticle->GetMother()<fStack->GetNprimary()) return McParticle->GetMother();
1245     else
1246     {
1247       Int_t motherlabel = McParticle->GetMother();
1248       while(motherlabel > -1)
1249       {
1250         AliMCParticle *MotherParticle  = (AliMCParticle*)fMC->GetTrack(motherlabel);
1251         if(MotherParticle->GetMother()<fStack->GetNprimary()) break;
1252         else motherlabel = MotherParticle->GetMother();
1253       }
1254       AliMCParticle *FirstSecondaryMotherParticle  = (AliMCParticle*)fMC->GetTrack(motherlabel);
1255       return FirstSecondaryMotherParticle->GetMother();
1256     }
1257   }
1258   else return -1;
1259 }
1260
1261 //________________________________________________________________________
1262 Double_t AliMuonEffMC::GetFirstPrimaryVertex(Int_t muonlabel)
1263 {
1264   if(fAOD) return 1.0;
1265   else if(fESD)
1266   {
1267     AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(muonlabel);
1268     if(McParticle->GetMother()<fStack->GetNprimary()) return McParticle->Zv();
1269     else
1270     {
1271      Int_t motherlabel = McParticle->GetMother();
1272       while(motherlabel > -1)
1273      {
1274        AliMCParticle *MotherParticle  = (AliMCParticle*)fMC->GetTrack(motherlabel);
1275        if(MotherParticle->GetMother()<fStack->GetNprimary()) break;
1276        else motherlabel = MotherParticle->GetMother();
1277      }
1278      AliMCParticle *FirtSecondaryMotherParticle  = (AliMCParticle*)fMC->GetTrack(motherlabel);
1279      return FirtSecondaryMotherParticle->Zv();
1280     }
1281   }
1282   else return -1;
1283 }