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),
  40   fOutputList(0x0),fHEventStat(0), fHXsec(0), fHTrials(0), fHEvt(0x0), fIsMc(kTRUE), fIsPythia(kFALSE), fMDProcess(kFALSE), fPlotMode(0),
  41   fCentralityEstimator("V0M"), fNEtaBins(100), fNpTBins(50), fNCentBins(1), fNZvtxBins(1), fNPhiBins(12),
  42   fHFPM(0x0), fHPP(0)
  43 {
  44   // Constructor
  45   //DefineInput(0, TChain::Class());
  46   //DefineOutput(1, TList::Class());
  47   for(Int_t i=0; i<2; i++)
  48   {
  49     fHDetRecMu[i] = NULL;
  50     fHDetRecMuFPM[i] = NULL;
  51     fHDetRecMuPP[i] = NULL;
  52     fHMuFPM[i] = NULL;
  53     fHMuPP[i] = NULL;
  54   }
  55   for(Int_t i=0; i<2; i++)
  56   {
  57     for(Int_t j=0; j<3; j++)
  58     {
  59       fHMuMotherRecPt[i][j] = NULL;
  60       fHMuMotherRecPhi[i][j] = NULL;
  61       fHMuMotherRecEta[i][j] = NULL;
  62        for(Int_t k=0; k<3; k++)
  63       {
  64         fHMuMohterPtDifRec[i][j][k] = NULL;
  65         fHMuMohterPhiDifRec[i][j][k] = NULL;
  66         fHMuMohterEtaDifRec[i][j][k] = NULL;
  67       }
  68     }
  69   }
  70   for(Int_t i=0; i<3; i++)
  71   {
  72     fHZvRv[i] = NULL;
  73     fHXvYv[i] = NULL;
  74   }
  75 }
  76
  77 //________________________________________________________________________
  78 AliMuonEffMC::AliMuonEffMC(const char *name) :
  79   AliAnalysisTaskSE(name), fESD(0), fAOD(0), fMC(0), fStack(0), fCentrality(99), fZVertex(99),
  80   fOutputList(0x0),fHEventStat(0), fHXsec(0), fHTrials(0), fHEvt(0x0), fIsMc(kTRUE), fIsPythia(kFALSE), fMDProcess(kFALSE), fPlotMode(0),
  81   fCentralityEstimator("V0M"), fNEtaBins(100), fNpTBins(50), fNCentBins(1), fNZvtxBins(1), fNPhiBins(12),
  82   fHFPM(0x0), fHPP(0)
  83 {
  84   // Constructor
  85   for(Int_t i=0; i<2; i++)
  86   {
  87     fHDetRecMu[i] = NULL;
  88     fHDetRecMuFPM[i] = NULL;
  89     fHDetRecMuPP[i] = NULL;
  90     fHMuFPM[i] = NULL;
  91     fHMuPP[i] = NULL;
  92   }
  93   for(Int_t i=0; i<2; i++)
  94   {
  95     for(Int_t j=0; j<3; j++)
  96     {
  97       fHMuMotherRecPt[i][j] = NULL;
  98       fHMuMotherRecPhi[i][j] = NULL;
  99       fHMuMotherRecEta[i][j] = NULL;
 100       for(Int_t k=0; k<3; k++)
 101       {
 102         fHMuMohterPtDifRec[i][j][k] = NULL;
 103         fHMuMohterPhiDifRec[i][j][k] = NULL;
 104         fHMuMohterEtaDifRec[i][j][k] = NULL;
 105       }
 106     }
 107   }
 108   for(Int_t i=0; i<3; i++)
 109   {
 110     fHZvRv[i] = NULL;
 111     fHXvYv[i] = NULL;
 112   }
 113   DefineInput(0, TChain::Class());
 114   DefineOutput(1, TList::Class());
 115 }
 116
 117 //________________________________________________________________________
 118 AliMuonEffMC::~AliMuonEffMC()
 119 {
 120   //Destructor
 121   if(fOutputList) delete fOutputList;
 122 }
 123
 124 //________________________________________________________________________
 125 void AliMuonEffMC::UserCreateOutputObjects()
 126 {
 127   // Create histograms
 128   // Called once (per slave on PROOF!)
 129   fOutputList = new TList();
 130   fOutputList->SetOwner(1);
 131
 132   fHEventStat = new TH1D("fHEventStat","Event statistics for analysis",18,0,18);
 133   fHEventStat->GetXaxis()->SetBinLabel(1,"Event");
 134   fHEventStat->GetXaxis()->SetBinLabel(2,"SelectedEvent");
 135   fHEventStat->GetXaxis()->SetBinLabel(3,"File");
 136   fHEventStat->GetXaxis()->SetBinLabel(4,"fSPHighpt");  //!Global Trigger Single plus High p_T
 137   fHEventStat->GetXaxis()->SetBinLabel(5,"fSPAllpt");   //!Global Trigger Single plus All p_T
 138   fHEventStat->GetXaxis()->SetBinLabel(6,"fSMLowpt");   //!Global Trigger Single minus Low p_T
 139   fHEventStat->GetXaxis()->SetBinLabel(7,"fSMHighpt");  //!Global Trigger Single minus High p_T
 140   fHEventStat->GetXaxis()->SetBinLabel(8,"fSMAllpt");   //!Global Trigger Single minus All p_T
 141   fHEventStat->GetXaxis()->SetBinLabel(9,"fSULowpt");   //!Global Trigger Single undefined Low p_T
 142   fHEventStat->GetXaxis()->SetBinLabel(10,"fSUHighpt"); //!Global Trigger Single undefined High p_T
 143   fHEventStat->GetXaxis()->SetBinLabel(11,"fSUAllpt");  //!Global Trigger Single undefined All p_T
 144   fHEventStat->GetXaxis()->SetBinLabel(12,"fUSLowpt");  //!Global Trigger UnlikeSign pair Low p_T
 145   fHEventStat->GetXaxis()->SetBinLabel(13,"fUSHighpt"); //!Global Trigger UnlikeSign pair High p_T
 146   fHEventStat->GetXaxis()->SetBinLabel(14,"fUSAllpt");  //!Global Trigger UnlikeSign pair All p_T
 147   fHEventStat->GetXaxis()->SetBinLabel(15,"fLSLowpt");  //!Global Trigger LikeSign pair pair Low p_T
 148   fHEventStat->GetXaxis()->SetBinLabel(16,"fLSHighpt"); //!Global Trigger LikeSign pair pair High p_T
 149   fHEventStat->GetXaxis()->SetBinLabel(17,"fLSAllpt");  //!Global Trigger LikeSign pair pair All p_T
 150   fHEventStat->GetXaxis()->SetBinLabel(18,"fSPLowpt");   //!Global Trigger Single plus Low p_T
 151   fOutputList->Add(fHEventStat);
 152
 153   if(fIsPythia)
 154   {
 155     fHXsec = new TH1F("fHXsec", "Cross section from pyxsec.root", 1, 0, 1);
 156     fHXsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 157     fOutputList->Add(fHXsec);
 158
 159     fHTrials = new TH1F("fHTrials", "Number of Trials", 1, 0, 1);
 160     fHTrials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 161     fOutputList->Add(fHTrials);
 162   }
 163
 164   fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", 30, -15, 15, 103, -2, 101);
 165   fOutputList->Add(fHEvt);
 166
 167   // Define THn's
 168   Int_t iTrackBinFPM[7];
 169   Double_t* trackBinsFPM[7];
 170   const char* trackAxisTitleFPM[7];
 171
 172   Int_t iTrackBinPP[7];
 173   Double_t* trackBinsPP[7];
 174   const char* trackAxisTitlePP[7];
 175
 176   Int_t iTrackBinMu[7];
 177   Double_t* trackBinsMu[7];
 178   const char* trackAxisTitleMu[7];
 179
 180   Int_t iTrackBinMuFPM[6];
 181   Double_t* trackBinsMuFPM[6];
 182   const char* trackAxisTitleMuFPM[6];
 183
 184   Int_t iTrackBinMuPP[6];
 185   Double_t* trackBinsMuPP[6];
 186   const char* trackAxisTitleMuPP[6];
 187
 188   // eta
 189   Double_t etaBins[fNEtaBins+1];
 190   for(Int_t i=0; i<=fNEtaBins; i++) { etaBins[i] = (Double_t)(-5.0 + 10.0/fNEtaBins*i); }
 191   iTrackBinFPM[0] = fNEtaBins;
 192   trackBinsFPM[0] = etaBins;
 193   trackAxisTitleFPM[0] = "#eta";
 194
 195   iTrackBinPP[0] = fNEtaBins;
 196   trackBinsPP[0] = etaBins;
 197   trackAxisTitlePP[0] = "#eta";
 198
 199   iTrackBinMu[0] = fNEtaBins;
 200   trackBinsMu[0] = etaBins;
 201   trackAxisTitleMu[0] = "#eta";
 202
 203   iTrackBinMuFPM[0] = fNEtaBins;
 204   trackBinsMuFPM[0] = etaBins;
 205   trackAxisTitleMuFPM[0] = "#eta";
 206
 207   iTrackBinMuPP[0] = fNEtaBins;
 208   trackBinsMuPP[0] = etaBins;
 209   trackAxisTitleMuPP[0] = "#eta";
 210
 211   // p_T
 212   Double_t pTBins[fNpTBins+1];
 213   for(Int_t i=0; i<=fNpTBins; i++) { pTBins[i] = (Double_t)(5.0/fNpTBins * i); }
 214   iTrackBinFPM[1] = fNpTBins;
 215   trackBinsFPM[1] = pTBins;
 216   trackAxisTitleFPM[1] = "p_{T} (GeV/c)";
 217
 218   iTrackBinPP[1] = fNpTBins;
 219   trackBinsPP[1] = pTBins;
 220   trackAxisTitlePP[1] = "p_{T} (GeV/c)";
 221
 222   iTrackBinMu[1] = fNpTBins;
 223   trackBinsMu[1] = pTBins;
 224   trackAxisTitleMu[1] = "p_{T} (GeV/c)";
 225
 226   iTrackBinMuFPM[1] = fNpTBins;
 227   trackBinsMuFPM[1] = pTBins;
 228   trackAxisTitleMuFPM[1] = "p_{T} (GeV/c)";
 229
 230   iTrackBinMuPP[1] = fNpTBins;
 231   trackBinsMuPP[1] = pTBins;
 232   trackAxisTitleMuPP[1] = "p_{T} (GeV/c)";
 233
 234   // centrality
 235   Double_t CentBins[fNCentBins+1];
 236   for (Int_t i=0; i<=fNCentBins; i++) { CentBins[i] = (Double_t)(100.0/fNCentBins * i); }
 237   iTrackBinFPM[2] = fNCentBins;
 238   trackBinsFPM[2] = CentBins;
 239   trackAxisTitleFPM[2] = "Cent";
 240
 241   iTrackBinPP[2] = fNCentBins;
 242   trackBinsPP[2] = CentBins;
 243   trackAxisTitlePP[2] = "Cent";
 244
 245   iTrackBinMu[2] = fNCentBins;
 246   trackBinsMu[2] = CentBins;
 247   trackAxisTitleMu[2] = "Cent";
 248
 249   // Z-vertex
 250   Double_t ZvtxBins[fNZvtxBins+1];
 251   for(Int_t i=0; i<=fNZvtxBins; i++) { ZvtxBins[i] = (Double_t)(-10.0 + 20.0/fNZvtxBins * i); }
 252   iTrackBinFPM[3] = fNZvtxBins;
 253   trackBinsFPM[3] = ZvtxBins;
 254   trackAxisTitleFPM[3] = "Zvtx";
 255
 256   iTrackBinPP[3] = fNZvtxBins;
 257   trackBinsPP[3] = ZvtxBins;
 258   trackAxisTitlePP[3] = "Zvtx";
 259
 260   iTrackBinMu[3] = fNZvtxBins;
 261   trackBinsMu[3] = ZvtxBins;
 262   trackAxisTitleMu[3] = "Zvtx";
 263
 264   // phi
 265   Double_t phiBins[fNPhiBins+1];
 266   for(Int_t i=0; i<=fNPhiBins; i++) { phiBins[i] = (Double_t)(TMath::TwoPi()/fNPhiBins * i); }
 267   iTrackBinFPM[4] = fNPhiBins;
 268   trackBinsFPM[4] = phiBins;
 269   trackAxisTitleFPM[4] = "#phi";
 270
 271   iTrackBinPP[4] = fNPhiBins;
 272   trackBinsPP[4] = phiBins;
 273   trackAxisTitlePP[4] = "#phi";
 274
 275   iTrackBinMu[4] = fNPhiBins;
 276   trackBinsMu[4] = phiBins;
 277   trackAxisTitleMu[4] = "#phi";
 278
 279   iTrackBinMuFPM[2] = fNPhiBins;
 280   trackBinsMuFPM[2] = phiBins;
 281   trackAxisTitleMuFPM[2] = "#phi";
 282
 283   iTrackBinMuPP[2] = fNPhiBins;
 284   trackBinsMuPP[2] = phiBins;
 285   trackAxisTitleMuPP[2] = "#phi";
 286
 287   // charge
 288   Double_t chargeBins[4] = {-10.0, -0.5, 0.5, 10.0};
 289   iTrackBinFPM[5] = 3;
 290   trackBinsFPM[5] = chargeBins;
 291   trackAxisTitleFPM[5] = "charge";
 292
 293   iTrackBinPP[5] = 3;
 294   trackBinsPP[5] = chargeBins;
 295   trackAxisTitlePP[5] = "charge";
 296
 297   iTrackBinMu[5] = 3;
 298   trackBinsMu[5] = chargeBins;
 299   trackAxisTitleMu[5] = "charge";
 300
 301   iTrackBinMuFPM[3] = 3;
 302   trackBinsMuFPM[3] = chargeBins;
 303   trackAxisTitleMuFPM[3] = "charge";
 304
 305   iTrackBinMuPP[3] = 3;
 306   trackBinsMuPP[3] = chargeBins;
 307   trackAxisTitleMuPP[3] = "charge";
 308
 309   // Muon type
 310   Double_t MuSpeciesBins[4] = {0.0, 1.0, 2.0, 3.0};
 311   iTrackBinMu[6] = 3;
 312   trackBinsMu[6] = MuSpeciesBins;
 313   trackAxisTitleMu[6] = "MUON type";
 314
 315   iTrackBinMuFPM[4] = 3;
 316   trackBinsMuFPM[4] = MuSpeciesBins;
 317   trackAxisTitleMuFPM[4] = "MUON type";
 318
 319   iTrackBinMuPP[4] = 3;
 320   trackBinsMuPP[4] = MuSpeciesBins;
 321   trackAxisTitleMuPP[4] = "MUON type";
 322
 323   // FPM species
 324   Double_t FPMSpecies[8] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0};
 325   iTrackBinFPM[6] = 7;
 326   trackBinsFPM[6] = FPMSpecies;
 327   trackAxisTitleFPM[6] = "FPM species";
 328
 329   iTrackBinMuFPM[5] = 7;
 330   trackBinsMuFPM[5] = FPMSpecies;
 331   trackAxisTitleMuFPM[5] = "FPM species";
 332
 333   // First Physical Primary  species
 334   Double_t PPMSpecies[8] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0};
 335   iTrackBinPP[6] = 7;
 336   trackBinsPP[6] = PPMSpecies;
 337   trackAxisTitlePP[6] = "PPM species";
 338
 339   iTrackBinMuPP[5] = 7;
 340   trackBinsMuPP[5] = PPMSpecies;
 341   trackAxisTitleMuPP[5] = "PPM species";
 342
 343   const char* MuonType[3] = {"Prim","Sec","Had"};
 344   const char *MuPt[3] = {"0005","0520","2040"};
 345   const char *cutlabel[2] = {"cut1", "cut2"};
 346
 347   if(fIsMc)
 348   {
 349     // THn for tracking efficiency
 350     if(fPlotMode==0)
 351     {
 352       fHFPM = new THnF("fHFPM", "", 7, iTrackBinFPM, 0, 0);
 353       for (Int_t i=0; i<7; i++)
 354       {
 355         fHFPM->SetBinEdges(i, trackBinsFPM[i]);
 356         fHFPM->GetAxis(i)->SetTitle(trackAxisTitleFPM[i]);
 357       }
 358       fOutputList->Add(fHFPM);
 359
 360       fHPP = new THnF("fHPP", "", 7, iTrackBinPP, 0, 0);
 361       for (Int_t i=0; i<7; i++)
 362       {
 363         fHPP->SetBinEdges(i, trackBinsPP[i]);
 364         fHPP->GetAxis(i)->SetTitle(trackAxisTitlePP[i]);
 365       }
 366       fOutputList->Add(fHPP);
 367     }
 368
 369     for(Int_t j=0; j<2; j++)
 370     {
 371       if(fPlotMode==1)
 372       {
 373         fHDetRecMu[j] = new THnF(Form("fHDetRecMu_%s",cutlabel[j]),"", 7, iTrackBinMu, 0, 0);
 374         for (Int_t i=0; i<7; i++)
 375         {
 376           fHDetRecMu[j]->SetBinEdges(i, trackBinsMu[i]);
 377           fHDetRecMu[j]->GetAxis(i)->SetTitle(trackAxisTitleMu[i]);
 378         }
 379         fOutputList->Add(fHDetRecMu[j]);
 380       }
 381       if(fPlotMode==2)
 382       {
 383         fHDetRecMuFPM[j] = new THnF(Form("fHDetRecMuFPM_%s",cutlabel[j]),"", 6, iTrackBinMuFPM, 0, 0);
 384         for (Int_t i=0; i<6; i++)
 385         {
 386           fHDetRecMuFPM[j]->SetBinEdges(i, trackBinsMuFPM[i]);
 387           fHDetRecMuFPM[j]->GetAxis(i)->SetTitle(trackAxisTitleMuFPM[i]);
 388         }
 389         fOutputList->Add(fHDetRecMuFPM[j]);
 390
 391         fHDetRecMuPP[j] = new THnF(Form("fHDetRecMuPP_%s",cutlabel[j]),"", 6, iTrackBinMuPP, 0, 0);
 392         for (Int_t i=0; i<6; i++)
 393         {
 394           fHDetRecMuPP[j]->SetBinEdges(i, trackBinsMuPP[i]);
 395           fHDetRecMuPP[j]->GetAxis(i)->SetTitle(trackAxisTitleMuPP[i]);
 396         }
 397         fOutputList->Add(fHDetRecMuPP[j]);
 398       }
 399       if(fPlotMode==3)
 400       {
 401         fHMuFPM[j] = new THnF(Form("fHMuFPM_%s",cutlabel[j]),"", 6, iTrackBinMuFPM, 0, 0);
 402         for (Int_t i=0; i<6; i++)
 403         {
 404           fHMuFPM[j]->SetBinEdges(i, trackBinsMuFPM[i]);
 405           fHMuFPM[j]->GetAxis(i)->SetTitle(trackAxisTitleMuFPM[i]);
 406         }
 407         fOutputList->Add(fHMuFPM[j]);
 408
 409         fHMuPP[j] = new THnF(Form("fHMuPP_%s",cutlabel[j]),"", 6, iTrackBinMuPP, 0, 0);
 410         for (Int_t i=0; i<6; i++)
 411         {
 412           fHMuPP[j]->SetBinEdges(i, trackBinsMuPP[i]);
 413           fHMuPP[j]->GetAxis(i)->SetTitle(trackAxisTitleMuPP[i]);
 414         }
 415         fOutputList->Add(fHMuPP[j]);
 416       }
 417     }
 418
 419     if(fMDProcess)
 420     {
 421       for(Int_t i=0; i<2; i++)
 422       {
 423         for(Int_t j=0; j<3; j++)
 424         {
 425           fHMuMotherRecPt[i][j] = new TH2F(Form("fHMuMotherRecPt_%s_%s",cutlabel[i], MuonType[j]),";p_{T,muon}^{rec} (GeV/c);p_{T,mother}^{Truth} (GeV/c);",500, 0, 50, 500, 0, 50);
 426           fOutputList->Add(fHMuMotherRecPt[i][j]);
 427           fHMuMotherRecPhi[i][j] = new TH2F(Form("fHMuMotherRecPhi_%s_%s",cutlabel[i], MuonType[j]),";#phi_{rec};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
 428           fOutputList->Add(fHMuMotherRecPhi[i][j]);
 429           fHMuMotherRecEta[i][j] = new TH2F(Form("fHMuMotherRecEta_%s_%s",cutlabel[i], MuonType[j]),";#eta_{rec};mother #eta;",100, -5., -1., 100, -5., -1.);
 430           fOutputList->Add(fHMuMotherRecEta[i][j]);
 431
 432           for(Int_t k=0; k<3; k++)
 433           {
 434             fHMuMohterPtDifRec[i][j][k] = new TH1F(Form("fHMuMohterPtDifRec_%s_%s_%s",cutlabel[i], MuonType[j], MuPt[k]),";#Delta#phi",200, -10.0, 10.0);
 435             fOutputList->Add(fHMuMohterPtDifRec[i][j][k]);
 436
 437             fHMuMohterPhiDifRec[i][j][k] = new TH1F(Form("fHMuMohterPhiDifRec_%s_%s_%s",cutlabel[i], MuonType[j], MuPt[k]),";#Delta#phi",100, -1.0*TMath::Pi(), TMath::Pi());
 438             fOutputList->Add(fHMuMohterPhiDifRec[i][j][k]);
 439
 440             fHMuMohterEtaDifRec[i][j][k] = new TH1F(Form("fHMuMohterEtaDifRec_%s_%s_%s",cutlabel[i], MuonType[j], MuPt[k]),";#Delta#eta",100, -5.0, 5.0);
 441             fOutputList->Add(fHMuMohterEtaDifRec[i][j][k]);
 442           }
 443         }
 444       }
 445       for(Int_t i = 0; i<3; i++)
 446       {
 447         fHZvRv[i] = new TH2F(Form("fHZvRv_%s",MuonType[i]), "", 300, -500, 100, 200, 0, 800);
 448         fOutputList->Add(fHZvRv[i]);
 449         fHXvYv[i] = new TH2F(Form("fHXvYv_%s",MuonType[i]), "", 200, -500, 500, 200, -500, 500);
 450         fOutputList->Add(fHXvYv[i]);
 451       }
 452     }
 453   }
 454   else
 455   {
 456     for(Int_t j=0; j<2; j++)
 457     {
 458       fHDetRecMu[j] = new THnF(Form("fHDetRecMu_%s",cutlabel[j]),"", 7, iTrackBinMu, 0, 0);
 459       for (Int_t i=0; i<7; i++)
 460       {
 461         fHDetRecMu[j]->SetBinEdges(i, trackBinsMu[i]);
 462         fHDetRecMu[j]->GetAxis(i)->SetTitle(trackAxisTitleMu[i]);
 463       }
 464       fOutputList->Add(fHDetRecMu[j]);
 465     }
 466   }
 467   PostData(1, fOutputList);
 468 }
 469
 470 //________________________________________________________________________
 471 Bool_t AliMuonEffMC::Notify()
 472 {
 473   // Implemented Notify() to read the cross sections
 474   // and number of trials from pyxsec.root
 475   if(fIsPythia)
 476   {
 477     fHEventStat->Fill(2.5);
 478
 479     TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 480     Float_t xsection = 0;
 481     Float_t ftrials  = 1;
 482
 483     if(tree){
 484       TFile *curfile = tree->GetCurrentFile();
 485       if (!curfile) {
 486         Error("Notify","No current file");
 487         return kFALSE;
 488       }
 489
 490       AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 491       fHXsec->Fill("<#sigma>",xsection);
 492       fHTrials->Fill("#sum{ntrials}",ftrials);
 493     }
 494   }
 495   return kTRUE;
 496 }
 497 //________________________________________________________________________
 498 void AliMuonEffMC::UserExec(Option_t *)
 499 {
 500   // Main loop, Called for each event
 501   Int_t ntrks = 0; // number of tracks in an event
 502
 503   if(((TString)InputEvent()->IsA()->GetName())=="AliAODEvent")
 504   {
 505     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 506     if (!fAOD) { AliError("AOD event not found. Nothing done!"); return; }
 507     ntrks = fAOD->GetNTracks();
 508   }
 509   else
 510   {
 511     fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 512     if (!fESD) { AliError("ESD event not found. Nothing done!"); return; }
 513     ntrks = fESD->GetNumberOfMuonTracks();
 514   }
 515
 516   fHEventStat->Fill(0.5);
 517   if(fIsMc)
 518   {
 519     fMC = MCEvent();
 520     if (!fMC) { AliError("MC event not avaliable."); return; }
 521     fStack = fMC->Stack();
 522   }
 523
 524   // Centrality, vertex, other event variables...
 525   if(fAOD)
 526   {
 527     const AliAODVertex* vertex = fAOD->GetPrimaryVertex();
 528     fZVertex = vertex->GetZ();
 529     if(fAOD->GetCentrality())  fCentrality = fAOD->GetCentrality()->GetCentralityPercentile(fCentralityEstimator);
 530   }
 531   else if(fESD)
 532   {
 533     const AliESDVertex* vertex = fESD->GetPrimaryVertex();
 534     fZVertex = vertex->GetZ();
 535     if(fESD->GetCentrality()) fCentrality = fESD->GetCentrality()->GetCentralityPercentile(fCentralityEstimator);
 536   }
 537
 538   if ((fESD && !VertexOk(fESD)) || (fAOD && !VertexOk(fAOD))) {
 539     //AliInfo(Form("Event REJECTED. z = %.1f", fZVertex));
 540     return;
 541   }
 542   if (fCentrality > 100. || fCentrality < -1.5) {
 543     //AliInfo(Form("Event REJECTED. fCentrality = %.1f", fCentrality));
 544     return;
 545   }
 546
 547   if(fCentrality < 0) fCentrality = 0.5; //ad hoc centrality for pp
 548   // Fill Event histogram
 549   fHEvt->Fill(fZVertex, fCentrality);
 550   fHEventStat->Fill(1.5);
 551
 552   ULong64_t trigword = 0;
 553   if(fAOD) trigword=fAOD->GetTriggerMask();
 554   else if(fESD) trigword=fESD->GetTriggerMask();
 555
 556   if (trigword & 0x01) fHEventStat->Fill(17.5);
 557   if (trigword & 0x02) fHEventStat->Fill(3.5);
 558   if (trigword & 0x04) fHEventStat->Fill(4.5);
 559   if (trigword & 0x08) fHEventStat->Fill(5.5);
 560   if (trigword & 0x010) fHEventStat->Fill(6.5);
 561   if (trigword & 0x020) fHEventStat->Fill(7.5);
 562   if (trigword & 0x040) fHEventStat->Fill(8.5);
 563   if (trigword & 0x080) fHEventStat->Fill(9.5);
 564   if (trigword & 0x100) fHEventStat->Fill(10.5);
 565   if (trigword & 0x200) fHEventStat->Fill(11.5);
 566   if (trigword & 0x400) fHEventStat->Fill(12.5);
 567   if (trigword & 0x800) fHEventStat->Fill(13.5);
 568   if (trigword & 0x1000) fHEventStat->Fill(14.5);
 569   if (trigword & 0x2000) fHEventStat->Fill(15.5);
 570   if (trigword & 0x4000) fHEventStat->Fill(16.5);
 571
 572   if(fIsMc && fPlotMode==0)
 573   {
 574     Double_t PPEta = 0.0;
 575     Double_t PPPt = 0.0;
 576     Double_t PPPhi = 0.0;
 577     Double_t PPCharge = 0.0;
 578     Double_t PPSpecies = 0.0;
 579
 580     Double_t FPMEta = 0.0;
 581     Double_t FPMPt = 0.0;
 582     Double_t FPMPhi = 0.0;
 583     Double_t FPMCharge = 0.0;
 584     Double_t FPMSpecies = 0.0;
 585
 586     // generated level loop
 587     for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)
 588     {
 589       if(fESD)
 590       {
 591         AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(ipart);
 592         if(!fMC->IsPhysicalPrimary(ipart) || McParticle->Charge()==0) continue;
 593
 594         PPEta = McParticle->Eta();
 595         PPPt = McParticle->Pt();
 596         PPPhi = McParticle->Phi();
 597         PPCharge = McParticle->Charge();
 598         PPSpecies = GetSpecies(McParticle->PdgCode());
 599
 600         AliMCParticle *FPMParticle  = (AliMCParticle*)fMC->GetTrack(GetFirstPrimaryMother(ipart));
 601         if(PPSpecies==1.5 || PPSpecies==2.5 || PPSpecies==4.5)
 602         {
 603           if(ipart < fStack->GetNprimary())
 604           {
 605             FPMEta = PPEta;
 606             FPMPt = PPPt;
 607             FPMPhi = PPPhi;
 608             FPMCharge = PPCharge;
 609             FPMSpecies = PPSpecies;
 610           }
 611           else
 612           {
 613             FPMEta = FPMParticle->Eta();
 614             FPMPt = FPMParticle->Pt();
 615             FPMPhi = FPMParticle->Phi();
 616             FPMCharge = FPMParticle->Charge();
 617             FPMSpecies = GetSpecies(FPMParticle->PdgCode());
 618           }
 619         }
 620         else
 621         {
 622           FPMEta = FPMParticle->Eta();
 623           FPMPt = FPMParticle->Pt();
 624           FPMPhi = FPMParticle->Phi();
 625           FPMCharge = FPMParticle->Charge();
 626           FPMSpecies = GetSpecies(FPMParticle->PdgCode());
 627         }
 628       }
 629       if(fAOD)
 630       {
 631         AliAODMCParticle *AodMcParticle  = (AliAODMCParticle*)fMC->GetTrack(ipart);
 632         if(!fMC->IsPhysicalPrimary(ipart) || AodMcParticle->Charge()==0) continue;
 633
 634         PPEta = AodMcParticle->Eta();
 635         PPPt = AodMcParticle->Pt();
 636         PPPhi = AodMcParticle->Phi();
 637         PPCharge = AodMcParticle->Charge();
 638         PPSpecies = GetSpecies(AodMcParticle->PdgCode());
 639
 640         AliAODMCParticle *AODFPMParticle  = (AliAODMCParticle*)fMC->GetTrack(GetFirstPrimaryMother(ipart));
 641         if(PPSpecies==1.5 || PPSpecies==2.5)
 642         {
 643           if(ipart < fStack->GetNprimary())
 644           {
 645             FPMEta = PPEta;
 646             FPMPt = PPPt;
 647             FPMPhi = PPPhi;
 648             FPMCharge = PPCharge;
 649             FPMSpecies = PPSpecies;
 650           }
 651           else
 652           {
 653             FPMEta = AODFPMParticle->Eta();
 654             FPMPt = AODFPMParticle->Pt();
 655             FPMPhi = AODFPMParticle->Phi();
 656             FPMCharge = AODFPMParticle->Charge();
 657             FPMSpecies = GetSpecies(AODFPMParticle->PdgCode());
 658           }
 659         }
 660         else
 661         {
 662           FPMEta = AODFPMParticle->Eta();
 663           FPMPt = AODFPMParticle->Pt();
 664           FPMPhi = AODFPMParticle->Phi();
 665           FPMCharge = AODFPMParticle->Charge();
 666           FPMSpecies = GetSpecies(AODFPMParticle->PdgCode());
 667         }
 668       }
 669       Double_t fillArrayPP[7] = { PPEta, PPPt, fCentrality, fZVertex, PPPhi, PPCharge, PPSpecies };
 670       fHPP->Fill(fillArrayPP);
 671       Double_t fillArrayFPM[7] = { FPMEta, FPMPt, fCentrality, fZVertex, FPMPhi, FPMCharge, FPMSpecies };
 672       fHFPM->Fill(fillArrayFPM);
 673     } // end of generated level loop
 674   }
 675
 676   // reconstructed level loop
 677   for (Int_t iTrack = 0; iTrack<ntrks; iTrack++)
 678   {
 679     Int_t label = 0;
 680     Int_t CutType = 0;
 681     Double_t MuonType = 0.0;
 682     Double_t FPMSpecies = 0.0;
 683     Double_t PPSpecies = 0.0;
 684
 685     Double_t RecEta = 0.0;
 686     Double_t RecPt = 0.0;
 687     Double_t RecPhi = 0.0;
 688     Double_t RecCharge = 0.0;
 689
 690     Double_t MuFPMEta = 0.0;
 691     Double_t MuFPMPt = 0.0;
 692     Double_t MuFPMPhi = 0.0;
 693     Double_t MuFPMCharge = 0.0;
 694
 695     Double_t MuPPEta = 0.0;
 696     Double_t MuPPPt = 0.0;
 697     Double_t MuPPPhi = 0.0;
 698     Double_t MuPPCharge = 0.0;
 699
 700     Double_t motherXv = 0.0;
 701     Double_t motherYv = 0.0;
 702     Double_t motherZv = 0.0;
 703
 704     if(fESD)
 705     {
 706       AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iTrack);
 707       if(muonTrack)
 708       {
 709         CutType = GetMUONCutType(*muonTrack);
 710         if(CutType > 1) continue;
 711         RecEta = muonTrack->Eta();
 712         RecPt = muonTrack->Pt();
 713         RecPhi = muonTrack->Phi();
 714         RecCharge = muonTrack->Charge();
 715         label =  TMath::Abs(muonTrack->GetLabel());
 716         if (label>=fMC->GetNumberOfTracks()) {
 717           AliError(Form("Label %d larger than number of particles on stack %d\n",label,fMC->GetNumberOfTracks()));
 718           continue;
 719         }
 720       }
 721     }
 722     else if(fAOD)
 723     {
 724       AliAODTrack* muonTrack = (AliAODTrack*)fAOD->GetTrack(iTrack);
 725       if(muonTrack)
 726       {
 727         if(!(muonTrack->IsMuonTrack())) continue;
 728         CutType =  GetMUONCutType(*muonTrack);
 729         if(CutType > 1) continue;
 730         RecEta = muonTrack->Eta();
 731         RecPt = muonTrack->Pt();
 732         RecPhi = muonTrack->Phi();
 733         RecCharge = muonTrack->Charge();
 734         label =  TMath::Abs(muonTrack->GetLabel());
 735         if (label>=fMC->GetNumberOfTracks()) {
 736           AliError(Form("Label %d larger than number of particles on stack %d\n",label,fMC->GetNumberOfTracks()));
 737           continue;
 738         }
 739       }
 740     }
 741
 742     if(fIsMc)
 743     {
 744       AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(label);
 745       MuonType = GetMuonTrackType(*McParticle);
 746
 747       if(GetFirstPrimaryMother(label) < 0) continue;
 748       AliMCParticle *MuFPMParticle  = (AliMCParticle*)fMC->GetTrack(GetFirstPrimaryMother(label));
 749       MuFPMEta = MuFPMParticle->Eta();
 750       MuFPMPt = MuFPMParticle->Pt();
 751       MuFPMPhi = MuFPMParticle->Phi();
 752       MuFPMCharge = MuFPMParticle->Charge();
 753       FPMSpecies = GetSpecies(MuFPMParticle->PdgCode());
 754
 755       if(GetFirstPPMother(label) < 0) continue;
 756       AliMCParticle *MuPPParticle  = (AliMCParticle*)fMC->GetTrack(GetFirstPPMother(label));
 757       MuPPEta = MuPPParticle->Eta();
 758       MuPPPt = MuPPParticle->Pt();
 759       MuPPPhi = MuPPParticle->Phi();
 760       MuPPCharge = MuPPParticle->Charge();
 761       PPSpecies = GetSpecies(MuPPParticle->PdgCode());
 762
 763       if(MuFPMParticle->GetFirstDaughter() > 0)
 764       {
 765         AliMCParticle *DaughtParticle  = (AliMCParticle*)fMC->GetTrack(MuFPMParticle->GetFirstDaughter());
 766         motherXv = DaughtParticle->Xv();
 767         motherYv = DaughtParticle->Yv();
 768         motherZv = DaughtParticle->Zv();
 769       }
 770
 771       if(fPlotMode==1)
 772       {
 773         Double_t fillArrayDetRecMu[7] = { RecEta, RecPt, fCentrality, fZVertex, RecPhi, RecCharge, MuonType };
 774         fHDetRecMu[CutType]->Fill(fillArrayDetRecMu);
 775       }
 776       if(fPlotMode==2)
 777       {
 778         Double_t fillArrayDetRecMuFPM[6] = { RecEta, RecPt, RecPhi, RecCharge, MuonType, FPMSpecies };
 779         fHDetRecMuFPM[CutType]->Fill(fillArrayDetRecMuFPM);
 780
 781         Double_t fillArrayDetRecMuPP[6] = { RecEta, RecPt, RecPhi, RecCharge, MuonType, PPSpecies };
 782         fHDetRecMuPP[CutType]->Fill(fillArrayDetRecMuPP);
 783       }
 784       if(fPlotMode==3)
 785       {
 786         Double_t fillArrayMuFPM[6] = { MuFPMEta, MuFPMPt, MuFPMPhi, MuFPMCharge, MuonType, FPMSpecies };
 787         fHMuFPM[CutType]->Fill(fillArrayMuFPM);
 788
 789         Double_t fillArrayMuPP[6] = { MuPPEta, MuPPPt, MuPPPhi, MuPPCharge, MuonType, PPSpecies };
 790         fHMuPP[CutType]->Fill(fillArrayMuPP);
 791       }
 792
 793       // mother-daughter kinematic relation
 794       if(fMDProcess)
 795       {
 796         if(MuFPMParticle->GetFirstDaughter() > 0)
 797         {
 798           if(CutType==0 || CutType==1)
 799           {
 800             fHZvRv[(Int_t)(MuonType - 0.5)]->Fill(motherZv, TMath::Sqrt(motherXv*motherXv + motherYv*motherYv));
 801             fHXvYv[(Int_t)(MuonType - 0.5)]->Fill(motherXv, motherYv);
 802           }
 803         }
 804         MDProcess((Int_t)(MuonType - 0.5), (Int_t)(CutType - 0.5), RecPt, RecPhi, RecEta, MuFPMPt, MuFPMPhi, MuFPMEta);
 805       }
 806     }// end of MC process
 807   }// end of reconstructed loop
 808   PostData(1, fOutputList);
 809   return;
 810 }
 811
 812 //________________________________________________________________________
 813 void AliMuonEffMC::Terminate(Option_t *)
 814 {
 815   // Draw result to the screen
 816   // Called once at the end of the query
 817 }
 818 //________________________________________________________________________
 819 Bool_t AliMuonEffMC::VertexOk(TObject* obj) const
 820 {
 821   // Modified from AliAnalyseLeadingTrackUE::VertexSelection()
 822
 823   Int_t nContributors  = 0;
 824   Double_t zVertex     = 999;
 825   TString name("");
 826
 827   if (obj->InheritsFrom("AliESDEvent")) {
 828     AliESDEvent* esdevt = (AliESDEvent*) obj;
 829     const AliESDVertex* vtx = esdevt->GetPrimaryVertex();
 830     if (!vtx)
 831       return 0;
 832     nContributors = vtx->GetNContributors();
 833     zVertex       = vtx->GetZ();
 834     name          = vtx->GetName();
 835   }
 836   else if (obj->InheritsFrom("AliAODEvent")) {
 837     AliAODEvent* aodevt = (AliAODEvent*) obj;
 838     if (aodevt->GetNumberOfVertices() < 1)
 839       return 0;
 840     const AliAODVertex* vtx = aodevt->GetPrimaryVertex();
 841     nContributors = vtx->GetNContributors();
 842     zVertex       = vtx->GetZ();
 843     name          = vtx->GetName();
 844   }
 845
 846   // Reject if TPC-only vertex
 847   if (name.CompareTo("TPCVertex")==0)
 848     return kFALSE;
 849
 850   // Check # contributors and range...
 851   if( nContributors < 1 || TMath::Abs(zVertex) > 10 ) {
 852     return kFALSE;
 853   }
 854
 855   return kTRUE;
 856 }
 857
 858 //________________________________________________________________________
 859 Int_t AliMuonEffMC::GetMUONCutType(AliESDMuonTrack &track)
 860 {
 861   Int_t cutNum = 4;
 862   Double_t thetaTrackAbsEnd = TMath::ATan(track.GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
 863   Double_t eta = track.Eta();
 864
 865   if(track.ContainTrackerData()) cutNum = 3;
 866   if(track.ContainTrackerData() && eta > -4. && -2.5 > eta) cutNum = 2;
 867   if(track.ContainTrackerData() && eta > -4. && -2.5 > eta && thetaTrackAbsEnd > 2. &&  10. > thetaTrackAbsEnd) cutNum =1;
 868   if(track.ContainTrackerData() && eta > -4. && -2.5 > eta && thetaTrackAbsEnd > 2. &&  10. > thetaTrackAbsEnd && track.GetMatchTrigger() > 0) cutNum = 0;
 869   return cutNum;
 870 }
 871
 872 //________________________________________________________________________
 873 Int_t AliMuonEffMC::GetMUONCutType(AliAODTrack &track)
 874 {
 875   Int_t cutNum = 4;
 876   Double_t thetaTrackAbsEnd = TMath::ATan(track.GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
 877   Double_t eta = track.Eta();
 878
 879   if(track.IsMuonTrack()) cutNum = 3;
 880   if(track.IsMuonTrack() && eta > -4. && -2.5 > eta) cutNum = 2;
 881   if(track.IsMuonTrack() && eta > -4. && -2.5 > eta && thetaTrackAbsEnd > 2. &&  10. > thetaTrackAbsEnd) cutNum = 1;
 882   if(track.IsMuonTrack() && eta > -4. && -2.5 > eta && thetaTrackAbsEnd > 2. &&  10. > thetaTrackAbsEnd && track.GetMatchTrigger() > 0) cutNum = 0;
 883
 884   return cutNum;
 885 }
 886
 887 //________________________________________________________________________
 888 Double_t AliMuonEffMC::GetMuonTrackType(AliMCParticle &track)
 889 {
 890   if(track.GetMother() < fStack->GetNprimary() && track.PdgCode() == 13) return 0.5;
 891   else if(track.GetMother() >= fStack->GetNprimary() && track.PdgCode() == 13) return 1.5;
 892   else return 2.5;
 893 }
 894
 895 //________________________________________________________________________
 896 Double_t AliMuonEffMC::GetSpecies(Int_t PdgCode)
 897 {
 898   Int_t code = TMath::Abs(PdgCode);
 899   if(code==13) return 0.5;
 900   else if(code==211) return 1.5;
 901   else if(code==321) return 2.5;
 902   else if(code==411 || code==413 || code==421 || code==423 || code==431 || code==433 || code==10413 || code==10411 || code==10423 || code==10421 || code==10433 || code==10431 || code==20413 || code==415 || code==20423 || code==425 || code==20433 || code==435) return 3.5;
 903   else if(code==2212) return 4.5;
 904   else if(code==11) return 5.5;
 905   else return 6.5;
 906 }
 907
 908 //________________________________________________________________________
 909 void AliMuonEffMC::MDProcess(Int_t isprimary, Int_t cutNum, Double_t trackpt, Double_t trackphi, Double_t tracketa, Double_t motherpt, Double_t motherphi, Double_t mothereta)
 910 {
 911   Int_t recptbin = -1;
 912
 913   if((0. <= trackpt) && (trackpt < 0.5)) recptbin = 0;
 914   else if((0.5 <= trackpt) && (trackpt < 2.0)) recptbin = 1;
 915   else recptbin = 2;
 916
 917   fHMuMotherRecPt[cutNum][isprimary]->Fill(trackpt, motherpt);
 918   fHMuMotherRecPhi[cutNum][isprimary]->Fill(trackphi, motherphi);
 919   fHMuMotherRecEta[cutNum][isprimary]->Fill(tracketa, mothereta);
 920
 921   fHMuMohterPtDifRec[cutNum][isprimary][recptbin]->Fill(motherpt-trackpt);
 922   fHMuMohterPhiDifRec[cutNum][isprimary][recptbin]->Fill(deltaphi(motherphi-trackphi));
 923   fHMuMohterEtaDifRec[cutNum][isprimary][recptbin]->Fill(mothereta-tracketa);
 924 }
 925
 926 //________________________________________________________________________
 927 Double_t AliMuonEffMC::deltaphi(Double_t phi)
 928 {
 929   if(phi < -1.0*TMath::Pi()) { return (phi + TMath::TwoPi()); }
 930   else if(phi > TMath::Pi()) { return (phi - TMath::TwoPi()); }
 931   else { return phi; }
 932 }
 933
 934 //________________________________________________________________________
 935 Int_t AliMuonEffMC::GetFirstPrimaryMother(Int_t muonlabel)
 936 {
 937   if(fAOD) return 1;
 938   else if(fESD)
 939   {
 940     AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(muonlabel);
 941     if(McParticle->GetMother()<fStack->GetNprimary()) return McParticle->GetMother();
 942     else
 943     {
 944       Int_t motherlabel = McParticle->GetMother();
 945       while(motherlabel > -1)
 946       {
 947         AliMCParticle *MotherParticle  = (AliMCParticle*)fMC->GetTrack(motherlabel);
 948         if(MotherParticle->GetMother()<fStack->GetNprimary()) break;
 949         else motherlabel = MotherParticle->GetMother();
 950       }
 951       AliMCParticle *FirstSecondaryMotherParticle  = (AliMCParticle*)fMC->GetTrack(motherlabel);
 952       return FirstSecondaryMotherParticle->GetMother();
 953     }
 954   }
 955   else return -1;
 956 }
 957
 958 //________________________________________________________________________
 959 Int_t AliMuonEffMC::GetFirstPPMother(Int_t muonlabel)
 960 {
 961   if(fAOD) return 1;
 962   else if(fESD)
 963   {
 964     AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(muonlabel);
 965     if(fMC->IsPhysicalPrimary(muonlabel)) return muonlabel;
 966     else
 967     {
 968       Int_t motherlabel = McParticle->GetMother();
 969       while(motherlabel > -1)
 970       {
 971         AliMCParticle *MotherParticle  = (AliMCParticle*)fMC->GetTrack(motherlabel);
 972         if(fMC->IsPhysicalPrimary(motherlabel)) break;
 973         else motherlabel = MotherParticle->GetMother();
 974       }
 975       return motherlabel;
 976     }
 977   }
 978   else return -1;
 979 }