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
  14 #include "TParticle.h"
  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 "AliAnalysisHelperJetTasks.h"
  29 #include "AliMCEvent.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), fHEvt(0x0), fIsMc(kFALSE), fMDProcess(kFALSE), fCentralityEstimator("V0M"),
  41   fNEtaBins(15), fNpTBins(100), fNCentBins(5), fNZvtxBins(10), fNPhiBins(36),
  42   fHMuonParGen(0x0), fHMuonDetGen(0x0), fHMuonDetRec(0x0), fHEtcDetRec(0x0), fHFXu(0), fHFXantiu(0), fHFXd(0), fHFXantid(0), fHFXg(0), fHFXetc(0), fHFXmuonP(0), fHFXmuonM(0)
  43 {
  44   // Constructor
  45   //DefineInput(0, TChain::Class());
  46   //DefineOutput(1, TList::Class());
  47
  48
  49 }
  50
  51 //________________________________________________________________________
  52 AliMuonEffMC::AliMuonEffMC(const char *name) :
  53   AliAnalysisTaskSE(name), fESD(0), fAOD(0), fMC(0), fStack(0), fCentrality(99), fZVertex(99), fOutputList(0x0),
  54   fHEventStat(0), fHEvt(0x0),  fIsMc(kFALSE), fMDProcess(kFALSE), fCentralityEstimator("V0M"),
  55   fNEtaBins(15), fNpTBins(100), fNCentBins(5), fNZvtxBins(10), fNPhiBins(36),
  56   fHMuonParGen(0x0), fHMuonDetGen(0x0), fHMuonDetRec(0x0), fHEtcDetRec(0x0), fHFXu(0), fHFXantiu(0), fHFXd(0), fHFXantid(0), fHFXg(0), fHFXetc(0), fHFXmuonP(0), fHFXmuonM(0)
  57 {
  58   // Constructor
  59   DefineInput(0, TChain::Class());
  60   DefineOutput(1, TList::Class());
  61
  62 }
  63
  64 //________________________________________________________________________
  65 AliMuonEffMC::~AliMuonEffMC()
  66 {
  67   //Destructor
  68   if(fOutputList) delete fOutputList;
  69 }
  70
  71 //________________________________________________________________________
  72 void AliMuonEffMC::UserCreateOutputObjects()
  73 {
  74   // Create histograms
  75   // Called once (per slave on PROOF!)
  76
  77   fOutputList = new TList();
  78   fOutputList->SetOwner(1);
  79
  80   for(Int_t i=0; i<4; i++)
  81   {
  82     fHMuMotherGenPt[i] = 0x0;
  83     fHMuMotherRecPt[i] = 0x0;
  84     fHMuMotherGenPhi[i] = 0x0;
  85     fHMuMotherRecPhi[i] = 0x0;
  86     fHMuMotherGenEta[i] = 0x0;
  87     fHMuMotherRecEta[i] = 0x0;
  88     fHMuDCA[i] = 0x0;
  89     for(Int_t j=0; j<3; j++)
  90     {
  91       fHMuMohterPhiDifGen[i][j] = 0x0;
  92       fHMuMohterPhiDifRec[i][j] = 0x0;
  93       fHMuMohterEtaDifGen[i][j] = 0x0;
  94       fHMuMohterEtaDifRec[i][j] = 0x0;
  95     }
  96   }
  97
  98   fHEventStat = new TH1D("fHEventStat","Event statistics for analysis",18,0,18);
  99   fHEventStat->GetXaxis()->SetBinLabel(1,"Event");
 100   fHEventStat->GetXaxis()->SetBinLabel(2,"SelectedEvent");
 101   fHEventStat->GetXaxis()->SetBinLabel(3,"File");
 102   fHEventStat->GetXaxis()->SetBinLabel(4,"fSPHighpt");  //!Global Trigger Single plus High p_T
 103   fHEventStat->GetXaxis()->SetBinLabel(5,"fSPAllpt");   //!Global Trigger Single plus All p_T
 104   fHEventStat->GetXaxis()->SetBinLabel(6,"fSMLowpt");   //!Global Trigger Single minus Low p_T
 105   fHEventStat->GetXaxis()->SetBinLabel(7,"fSMHighpt");  //!Global Trigger Single minus High p_T
 106   fHEventStat->GetXaxis()->SetBinLabel(8,"fSMAllpt");   //!Global Trigger Single minus All p_T
 107   fHEventStat->GetXaxis()->SetBinLabel(9,"fSULowpt");   //!Global Trigger Single undefined Low p_T
 108   fHEventStat->GetXaxis()->SetBinLabel(10,"fSUHighpt"); //!Global Trigger Single undefined High p_T
 109   fHEventStat->GetXaxis()->SetBinLabel(11,"fSUAllpt");  //!Global Trigger Single undefined All p_T
 110   fHEventStat->GetXaxis()->SetBinLabel(12,"fUSLowpt");  //!Global Trigger UnlikeSign pair Low p_T
 111   fHEventStat->GetXaxis()->SetBinLabel(13,"fUSHighpt"); //!Global Trigger UnlikeSign pair High p_T
 112   fHEventStat->GetXaxis()->SetBinLabel(14,"fUSAllpt");  //!Global Trigger UnlikeSign pair All p_T
 113   fHEventStat->GetXaxis()->SetBinLabel(15,"fLSLowpt");  //!Global Trigger LikeSign pair pair Low p_T
 114   fHEventStat->GetXaxis()->SetBinLabel(16,"fLSHighpt"); //!Global Trigger LikeSign pair pair High p_T
 115   fHEventStat->GetXaxis()->SetBinLabel(17,"fLSAllpt");  //!Global Trigger LikeSign pair pair All p_T
 116   fHEventStat->GetXaxis()->SetBinLabel(18,"fSPLowpt");   //!Global Trigger Single plus Low p_T
 117   fOutputList->Add(fHEventStat);
 118
 119   fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", 30, -15, 15, 103, -2, 101);
 120   fOutputList->Add(fHEvt);
 121
 122
 123   Int_t iTrackBin[5];
 124   Double_t* trackBins[5];
 125   const char* trackAxisTitle[5];
 126
 127   // eta
 128   Double_t etaBins[fNEtaBins+1];
 129   for (Int_t i=0; i<=fNEtaBins; i++) { etaBins[i] = (Double_t)(-4.0 + 1.5/fNEtaBins*i); }
 130   iTrackBin[0] = fNEtaBins;
 131   trackBins[0] = etaBins;
 132   trackAxisTitle[0] = "#eta";
 133
 134   // p_T
 135   Double_t pTBins[fNpTBins+1];
 136   for (Int_t i=0; i<=fNpTBins; i++) { pTBins[i] = (Double_t)(10.0/fNpTBins*i); }
 137   iTrackBin[1] = fNpTBins;
 138   trackBins[1] = pTBins;
 139   trackAxisTitle[1] = "p_{T} (GeV/c)";
 140
 141   // centrality
 142   Double_t CentBins[fNCentBins+1];
 143   for (Int_t i=0; i<=fNCentBins; i++) { CentBins[i] = (Double_t)(100.0/fNCentBins*i); }
 144   iTrackBin[2] = fNCentBins;
 145   trackBins[2] = CentBins;
 146   trackAxisTitle[2] = "Cent";
 147
 148   // Z-vertex
 149   Double_t ZvtxBins[fNZvtxBins+1];
 150   for (Int_t i=0; i<=fNZvtxBins; i++) { ZvtxBins[i] = (Double_t)(-10.0 + 20.0/fNZvtxBins*i); }
 151   iTrackBin[3] = fNZvtxBins;
 152   trackBins[3] = ZvtxBins;
 153   trackAxisTitle[3] = "Zvtx";
 154
 155   // phi
 156   Double_t phiBins[fNPhiBins+1];
 157   for (Int_t i=0; i<=fNPhiBins; i++) { phiBins[i] = (Double_t)(TMath::TwoPi()/fNPhiBins * i); }
 158   iTrackBin[4] = fNPhiBins;
 159   trackBins[4] = phiBins;
 160   trackAxisTitle[4] = "#phi";
 161
 162   // THn for tracking efficiency
 163   fHMuonParGen = new THnF("fHMuonParGen", "", 5, iTrackBin, 0, 0);
 164   for (Int_t i=0; i<5; i++)
 165   {
 166     fHMuonParGen->SetBinEdges(i, trackBins[i]);
 167     fHMuonParGen->GetAxis(i)->SetTitle(trackAxisTitle[i]);
 168   }
 169   fHMuonParGen->Sumw2();
 170   fOutputList->Add(fHMuonParGen);
 171
 172   fHMuonDetGen = (THnF*) fHMuonParGen->Clone("fHMuonDetGen");
 173   fHMuonDetGen->Sumw2();
 174   fOutputList->Add(fHMuonDetGen);
 175
 176   fHMuonDetRec = (THnF*) fHMuonParGen->Clone("fHMuonDetRec");
 177   fHMuonDetRec->Sumw2();
 178   fOutputList->Add(fHMuonDetRec);
 179
 180   fHEtcDetRec = (THnF*) fHMuonParGen->Clone("fHEtcDetRec");
 181   fHEtcDetRec->Sumw2();
 182   fOutputList->Add(fHEtcDetRec);
 183
 184   if(fIsMc)
 185   {
 186     const char* MotherSpecies[4] = {"Pion","Kaon","D", "Etc"};
 187     const char *MuPt[3] = {"0510","1020","2040"};
 188     if(fMDProcess)
 189     {
 190       for(Int_t i=0; i<4; i++)
 191       {
 192         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);
 193         fOutputList->Add(fHMuMotherGenPt[i]);
 194         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);
 195         fOutputList->Add(fHMuMotherRecPt[i]);
 196         fHMuMotherGenPhi[i] = new TH2F(Form("fHMuMotherGenPhi_%s",MotherSpecies[i]),";#phi_{gen};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
 197         fOutputList->Add(fHMuMotherGenPhi[i]);
 198         fHMuMotherRecPhi[i] = new TH2F(Form("fHMuMotherRecPhi_%s",MotherSpecies[i]),";#phi_{rec};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
 199         fOutputList->Add(fHMuMotherRecPhi[i]);
 200         fHMuMotherGenEta[i] = new TH2F(Form("fHMuMotherGenEta_%s",MotherSpecies[i]),";#eta_{gen};mother #eta;",100, -5., -1., 100, -5., -1.);
 201         fOutputList->Add(fHMuMotherGenEta[i]);
 202         fHMuMotherRecEta[i] = new TH2F(Form("fHMuMotherRecEta_%s",MotherSpecies[i]),";#eta_{rec};mother #eta;",100, -5., -1., 100, -5., -1.);
 203         fOutputList->Add(fHMuMotherRecEta[i]);
 204         fHMuDCA[i] =  new TH1F(Form("fHMuDCA_%s",MotherSpecies[i]), ";DCA", 100, 0, 50);
 205         fOutputList->Add(fHMuDCA[i]);
 206
 207         for(Int_t j=0; j<3; j++)
 208         {
 209           fHMuMohterPhiDifGen[i][j] = new TH1F(Form("fHMuMohterPhiDifGen_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#phi",100, -1.0*TMath::Pi(), TMath::Pi());
 210           fOutputList->Add(fHMuMohterPhiDifGen[i][j]);
 211           fHMuMohterPhiDifRec[i][j] = new TH1F(Form("fHMuMohterPhiDifRec_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#phi",100, -1.0*TMath::Pi(), TMath::Pi());
 212           fOutputList->Add(fHMuMohterPhiDifRec[i][j]);
 213           fHMuMohterEtaDifGen[i][j] = new TH1F(Form("fHMuMohterEtaDifGen_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#eta",100, -5.0, 5.0);
 214           fOutputList->Add(fHMuMohterEtaDifGen[i][j]);
 215           fHMuMohterEtaDifRec[i][j] = new TH1F(Form("fHMuMohterEtaDifRec_%s_%s",MotherSpecies[i], MuPt[j]),";#Delta#eta",100, -5.0, 5.0);
 216           fOutputList->Add(fHMuMohterEtaDifRec[i][j]);
 217         }
 218       }
 219     }
 220
 221     fHFXu = new TH1F("fHFXu",";x_{F}",1000, 0.0, 1.0);
 222     fOutputList->Add(fHFXu);
 223
 224     fHFXantiu = new TH1F("fHFXantiu",";x_{F}",1000, 0.0, 1.0);
 225     fOutputList->Add(fHFXantiu);
 226
 227     fHFXd = new TH1F("fHFXd",";x_{F}",1000, 0.0, 1.0);
 228     fOutputList->Add(fHFXd);
 229
 230     fHFXantid = new TH1F("fHFXantid",";x_{F}",1000, 0.0, 1.0);
 231     fOutputList->Add(fHFXantid);
 232
 233     fHFXg = new TH1F("fHFXg",";x_{F}",1000, 0.0, 1.0);
 234     fOutputList->Add(fHFXg);
 235
 236     fHFXetc = new TH1F("fHFXetc",";x_{F}",1000, 0.0, 1.0);
 237     fOutputList->Add(fHFXetc);
 238
 239     fHFXmuonP = new TH1F("fHFXmuonP",";x_{F}",1000, 0.0, 1.0);
 240     fOutputList->Add(fHFXmuonP);
 241
 242     fHFXmuonM = new TH1F("fHFXmuonM",";x_{F}",1000, 0.0, 1.0);
 243     fOutputList->Add(fHFXmuonM);
 244   }
 245   PostData(1, fOutputList);
 246 }
 247
 248 //________________________________________________________________________
 249 void AliMuonEffMC::UserExec(Option_t *)
 250 {
 251   // Main loop, Called for each event
 252   Int_t ntrks = 0; // number of tracks in an event
 253
 254   if(((TString)InputEvent()->IsA()->GetName())=="AliAODEvent")
 255   {
 256     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 257     if (!fAOD) { AliError("AOD event not found. Nothing done!"); return; }
 258     ntrks = fAOD->GetNTracks();
 259   }
 260   else
 261   {
 262     fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 263     if (!fESD) { AliError("ESD event not found. Nothing done!"); return; }
 264     ntrks = fESD->GetNumberOfMuonTracks();
 265   }
 266
 267   fHEventStat->Fill(0.5);
 268   if(fIsMc)
 269   {
 270     fMC = MCEvent();
 271   }
 272
 273   // Centrality, vertex, other event variables...
 274   if(fAOD)
 275   {
 276     const AliAODVertex* vertex = fAOD->GetPrimaryVertex();
 277     fZVertex = vertex->GetZ();
 278     if(fAOD->GetCentrality())  fCentrality = fAOD->GetCentrality()->GetCentralityPercentile(fCentralityEstimator);
 279   }
 280   else if(fESD)
 281   {
 282     const AliESDVertex* vertex = fESD->GetPrimaryVertex();
 283     fZVertex = vertex->GetZ();
 284     if(fESD->GetCentrality()) fCentrality = fESD->GetCentrality()->GetCentralityPercentile(fCentralityEstimator);
 285   }
 286
 287   if ((fESD && !VertexOk(fESD)) || (fAOD && !VertexOk(fAOD))) { //AliInfo(Form("Event REJECTED. z = %.1f", fZVertex));
 288     return; }
 289   if (fCentrality > 100. || fCentrality < -1.5) { //AliInfo(Form("Event REJECTED. fCentrality = %.1f", fCentrality));
 290     return; }
 291
 292   if(fCentrality < 0) fCentrality = 0.0; //ad hoc centrality for pp
 293   // Fill Event histogram
 294   fHEvt->Fill(fZVertex, fCentrality);
 295   fHEventStat->Fill(1.5);
 296
 297   ULong64_t trigword = 0;
 298   if(fAOD) trigword=fAOD->GetTriggerMask();
 299   else if(fESD) trigword=fESD->GetTriggerMask();
 300
 301   if (trigword & 0x01) fHEventStat->Fill(17.5);
 302   if (trigword & 0x02) fHEventStat->Fill(3.5);
 303   if (trigword & 0x04) fHEventStat->Fill(4.5);
 304   if (trigword & 0x08) fHEventStat->Fill(5.5);
 305   if (trigword & 0x010) fHEventStat->Fill(6.5);
 306   if (trigword & 0x020) fHEventStat->Fill(7.5);
 307   if (trigword & 0x040) fHEventStat->Fill(8.5);
 308   if (trigword & 0x080) fHEventStat->Fill(9.5);
 309   if (trigword & 0x100) fHEventStat->Fill(10.5);
 310   if (trigword & 0x200) fHEventStat->Fill(11.5);
 311   if (trigword & 0x400) fHEventStat->Fill(12.5);
 312   if (trigword & 0x800) fHEventStat->Fill(13.5);
 313   if (trigword & 0x1000) fHEventStat->Fill(14.5);
 314   if (trigword & 0x2000) fHEventStat->Fill(15.5);
 315   if (trigword & 0x4000) fHEventStat->Fill(16.5);
 316
 317   // generated level loop
 318   if(fIsMc)
 319   {
 320     for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)
 321     {
 322       if(fAOD)
 323       {
 324         AliAODMCParticle *AodMcParticle  = (AliAODMCParticle*)fMC->GetTrack(ipart);
 325         if(TMath::Abs(AodMcParticle->PdgCode())==13)
 326         {
 327           Double_t fillArrayParGen[5] = { AodMcParticle->Eta(), AodMcParticle->Pt(), fCentrality, fZVertex, AodMcParticle->Phi() };
 328           fHMuonParGen->Fill(fillArrayParGen);
 329         }
 330       }
 331       else if(fESD)
 332       {
 333         AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(ipart);
 334         if(TMath::Abs(McParticle->PdgCode())==13)
 335         {
 336           Double_t fillArrayParGen[5] = { McParticle->Eta(), McParticle->Pt(), fCentrality, fZVertex, McParticle->Phi() };
 337           fHMuonParGen->Fill(fillArrayParGen);
 338         }
 339       }
 340     }
 341   }
 342
 343   // reconstructed level loop
 344   for (Int_t iTrack = 0; iTrack<ntrks; iTrack++)
 345   {
 346     Int_t label = 0;
 347     Double_t trackpt = 0;
 348     Double_t tracketa = 0;
 349     Double_t trackphi = 0;
 350     Double_t dcavalue = 0;
 351
 352     Double_t mcpt = 0;
 353     Double_t mceta = 0;
 354     Double_t mcphi = 0;
 355
 356     Int_t motherlabel =0;
 357     Int_t motherpdg = 0;
 358     Double_t motherpt = 0;
 359     Double_t mothereta = 0;
 360     Double_t motherphi = 0;
 361
 362     if(fAOD)
 363     {
 364       AliAODTrack* muonTrack = (AliAODTrack*)fAOD->GetTrack(iTrack);
 365       if(muonTrack)
 366       {
 367         if(!(IsGoodMUONtrack(*muonTrack)) || !(muonTrack->IsMuonTrack())) continue;
 368         trackpt = muonTrack->Pt();
 369         tracketa = muonTrack->Eta();
 370         trackphi = muonTrack->Phi();
 371         label =  TMath::Abs(muonTrack->GetLabel());
 372       }
 373     }
 374     else if(fESD)
 375     {
 376       AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iTrack);
 377       if(muonTrack)
 378       {
 379         if(!IsGoodMUONtrack(*muonTrack)) continue;
 380         trackpt = muonTrack->Pt();
 381         tracketa = muonTrack->Eta();
 382         trackphi = muonTrack->Phi();
 383         label =  TMath::Abs(muonTrack->GetLabel());
 384         dcavalue = muonTrack->GetDCA();
 385       }
 386     }
 387     Double_t fillArrayDetRec[5] = { tracketa, trackpt, fCentrality, fZVertex, trackphi };
 388     fHMuonDetRec->Fill(fillArrayDetRec);
 389
 390     if(fIsMc)
 391     {
 392       if(fAOD)
 393       {
 394         AliAODMCParticle *aodMcParticle  = (AliAODMCParticle*)fMC->GetTrack(label);
 395         if(TMath::Abs(aodMcParticle->PdgCode()) != 13)
 396         {
 397           fHEtcDetRec->Fill(fillArrayDetRec);
 398           continue;
 399         }
 400         mcpt = aodMcParticle->Pt();
 401         mceta = aodMcParticle->Eta();
 402         mcphi = aodMcParticle->Phi();
 403         motherlabel = aodMcParticle->GetMother();
 404         if(motherlabel > 0)
 405         {
 406           AliAODMCParticle *aodMotherParticle  = (AliAODMCParticle*)fMC->GetTrack(motherlabel);
 407           motherpdg = TMath::Abs(aodMotherParticle->PdgCode());
 408           motherpt = aodMotherParticle->Pt();
 409           mothereta = aodMotherParticle->Eta();
 410           motherphi = aodMotherParticle->Phi();
 411         }
 412       }
 413       else if(fESD)
 414       {
 415         AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(label);
 416         if(TMath::Abs(McParticle->PdgCode()) != 13)
 417         {
 418           fHEtcDetRec->Fill(fillArrayDetRec);
 419           continue;
 420         }
 421         mcpt = McParticle->Pt();
 422         mceta = McParticle->Eta();
 423         mcphi = McParticle->Phi();
 424         motherlabel = McParticle->GetMother();
 425         if(motherlabel > 0)
 426         {
 427           AliMCParticle *MotherParticle  = (AliMCParticle*)fMC->GetTrack(motherlabel);
 428           motherpdg = TMath::Abs(MotherParticle->PdgCode());
 429           motherpt = MotherParticle->Pt();
 430           mothereta = MotherParticle->Eta();
 431           motherphi = MotherParticle->Phi();
 432         }
 433       }
 434       Double_t fillArrayDetGen[5] = { mceta, mcpt, fCentrality, fZVertex, mcphi };
 435       fHMuonDetGen->Fill(fillArrayDetGen);
 436
 437       // mother-daughter kinematic relation
 438       if(fMDProcess && motherlabel>0) MDProcess(motherpdg, mcpt, mcphi, mceta, trackpt, trackphi, tracketa, motherpt, motherphi, mothereta, dcavalue);
 439     }
 440   }
 441   PostData(1, fOutputList);
 442   return;
 443 }
 444
 445 //________________________________________________________________________
 446 void AliMuonEffMC::Terminate(Option_t *)
 447 {
 448   // Draw result to the screen
 449   // Called once at the end of the query
 450 }
 451 //________________________________________________________________________
 452 Bool_t AliMuonEffMC::VertexOk(TObject* obj) const
 453 {
 454   // Modified from AliAnalyseLeadingTrackUE::VertexSelection()
 455
 456   Int_t nContributors  = 0;
 457   Double_t zVertex     = 999;
 458   TString name("");
 459
 460   if (obj->InheritsFrom("AliESDEvent")) {
 461     AliESDEvent* esdevt = (AliESDEvent*) obj;
 462     const AliESDVertex* vtx = esdevt->GetPrimaryVertex();
 463     if (!vtx)
 464       return 0;
 465     nContributors = vtx->GetNContributors();
 466     zVertex       = vtx->GetZ();
 467     name          = vtx->GetName();
 468   }
 469   else if (obj->InheritsFrom("AliAODEvent")) {
 470     AliAODEvent* aodevt = (AliAODEvent*) obj;
 471     if (aodevt->GetNumberOfVertices() < 1)
 472       return 0;
 473     const AliAODVertex* vtx = aodevt->GetPrimaryVertex();
 474     nContributors = vtx->GetNContributors();
 475     zVertex       = vtx->GetZ();
 476     name          = vtx->GetName();
 477   }
 478
 479   // Reject if TPC-only vertex
 480   if (name.CompareTo("TPCVertex")==0)
 481     return kFALSE;
 482
 483   // Check # contributors and range...
 484   if( nContributors < 1 || TMath::Abs(zVertex) > 10 ) {
 485     return kFALSE;
 486   }
 487
 488   return kTRUE;
 489 }
 490
 491
 492 //________________________________________________________________________
 493 Bool_t AliMuonEffMC::IsGoodMUONtrack(AliESDMuonTrack &track)
 494 {
 495   // Applying track cuts for MUON tracks
 496   if(!track.ContainTrackerData()) return kFALSE;
 497
 498   Double_t thetaTrackAbsEnd = TMath::ATan(track.GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
 499   Double_t eta = track.Eta();
 500
 501   // Theta cut at absorber end
 502   if(thetaTrackAbsEnd <= 2. || thetaTrackAbsEnd >= 10.) return kFALSE;
 503   // Eta cut
 504   if(eta <= -4. || eta >= -2.5) return kFALSE;
 505   if(track.GetMatchTrigger() <= 0) return kFALSE;
 506
 507   return kTRUE;
 508 }
 509
 510 //________________________________________________________________________
 511 Bool_t AliMuonEffMC::IsGoodMUONtrack(AliAODTrack &track)
 512 {
 513   if (!track.IsMuonTrack()) return kFALSE;
 514
 515   Double_t dThetaAbs = TMath::ATan(track.GetRAtAbsorberEnd()/505.)
 516                      * TMath::RadToDeg();
 517   if ((dThetaAbs<2.) || (dThetaAbs>10.)) return kFALSE;
 518
 519   Double_t dEta = track.Eta();
 520   if ((dEta<-4.) || (dEta>2.5)) return kFALSE;
 521
 522   if (track.GetMatchTrigger()<0.5) return kFALSE;
 523
 524   return kTRUE;
 525 }
 526
 527 //________________________________________________________________________
 528 void AliMuonEffMC::MDProcess(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)
 529 {
 530   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)
 531   {
 532     fHMuMotherGenPt[2]->Fill(mcpt, motherpt);
 533     fHMuMotherRecPt[2]->Fill(trackpt, motherpt);
 534     fHMuMotherGenPhi[2]->Fill(mcphi, motherphi);
 535     fHMuMotherRecPhi[2]->Fill(trackphi, motherphi);
 536     fHMuMotherGenEta[2]->Fill(mceta, mothereta);
 537     fHMuMotherRecEta[2]->Fill(tracketa, mothereta);
 538
 539     // generated
 540     if((0.5 <= mcpt) && (mcpt < 1.0))
 541     {
 542       fHMuMohterPhiDifGen[2][0]->Fill(deltaphi(motherphi-mcphi));
 543       fHMuMohterEtaDifGen[2][0]->Fill(mothereta-mceta);
 544     }
 545     else if((1.0 <= mcpt) && (mcpt < 2.0))
 546     {
 547       fHMuMohterPhiDifGen[2][1]->Fill(deltaphi(motherphi-mcphi));
 548       fHMuMohterEtaDifGen[2][1]->Fill(mothereta-mceta);
 549     }
 550     else if((2.0 <= mcpt) && (mcpt < 4.0))
 551     {
 552       fHMuMohterPhiDifGen[2][2]->Fill(deltaphi(motherphi-mcphi));
 553       fHMuMohterEtaDifGen[2][2]->Fill(mothereta-mceta);
 554     }
 555     // reconstructed
 556     if((0.5 <= trackpt) && (trackpt < 1.0))
 557     {
 558       fHMuMohterPhiDifRec[2][0]->Fill(deltaphi(motherphi-trackphi));
 559       fHMuMohterEtaDifRec[2][0]->Fill(mothereta-tracketa);
 560     }
 561     else if((1.0 <= trackpt) && (trackpt < 2.0))
 562     {
 563       fHMuMohterPhiDifRec[2][1]->Fill(deltaphi(motherphi-trackphi));
 564       fHMuMohterEtaDifRec[2][1]->Fill(mothereta-tracketa);
 565     }
 566     else if((2.0 <= trackpt) && (trackpt < 4.0))
 567     {
 568       fHMuMohterPhiDifRec[2][2]->Fill(deltaphi(motherphi-trackphi));
 569       fHMuMohterEtaDifRec[2][2]->Fill(mothereta-tracketa);
 570     }
 571
 572     // DCA
 573     if(fESD) fHMuDCA[2]->Fill(dcavalue);
 574   }
 575
 576   else if(motherpdg==211)
 577   {
 578     fHMuMotherGenPt[0]->Fill(mcpt, motherpt);
 579     fHMuMotherRecPt[0]->Fill(trackpt, motherpt);
 580     fHMuMotherGenPhi[0]->Fill(mcphi, motherphi);
 581     fHMuMotherRecPhi[0]->Fill(trackphi, motherphi);
 582     fHMuMotherGenEta[0]->Fill(mceta, mothereta);
 583     fHMuMotherRecEta[0]->Fill(tracketa, mothereta);
 584     // generated
 585     if((0.5 <= mcpt) && (mcpt < 1.0))
 586     {
 587       fHMuMohterPhiDifGen[0][0]->Fill(deltaphi(motherphi-mcphi));
 588       fHMuMohterEtaDifGen[0][0]->Fill(mothereta-mceta);
 589     }
 590     else if((1.0 <= mcpt) && (mcpt < 2.0))
 591     {
 592       fHMuMohterPhiDifGen[0][1]->Fill(deltaphi(motherphi-mcphi));
 593       fHMuMohterEtaDifGen[0][1]->Fill(mothereta-mceta);
 594     }
 595     else if((2.0 <= mcpt) && (mcpt < 4.0))
 596     {
 597       fHMuMohterPhiDifGen[0][2]->Fill(deltaphi(motherphi-mcphi));
 598       fHMuMohterEtaDifGen[0][2]->Fill(mothereta-mceta);
 599     }
 600     // reconstructed
 601     if((0.5 <= trackpt) && (trackpt < 1.0))
 602     {
 603       fHMuMohterPhiDifRec[0][0]->Fill(deltaphi(motherphi-trackphi));
 604       fHMuMohterEtaDifRec[0][0]->Fill(mothereta-tracketa);
 605     }
 606     else if((1.0 <= trackpt) && (trackpt < 2.0))
 607     {
 608       fHMuMohterPhiDifRec[0][1]->Fill(deltaphi(motherphi-trackphi));
 609       fHMuMohterEtaDifRec[0][1]->Fill(mothereta-tracketa);
 610     }
 611     else if((2.0 <= trackpt) && (trackpt < 4.0))
 612     {
 613       fHMuMohterPhiDifRec[0][2]->Fill(deltaphi(motherphi-trackphi));
 614       fHMuMohterEtaDifRec[0][2]->Fill(mothereta-tracketa);
 615     }
 616     if(fESD) fHMuDCA[0]->Fill(dcavalue);
 617   }
 618   else if(motherpdg==321)
 619   {
 620     fHMuMotherGenPt[1]->Fill(mcpt, motherpt);
 621     fHMuMotherRecPt[1]->Fill(trackpt, motherpt);
 622     fHMuMotherGenPhi[1]->Fill(mcphi, motherphi);
 623     fHMuMotherRecPhi[1]->Fill(trackphi, motherphi);
 624     fHMuMotherGenEta[1]->Fill(mceta, mothereta);
 625     fHMuMotherRecEta[1]->Fill(tracketa, mothereta);
 626     // generated
 627     if((0.5 <= mcpt) && (mcpt < 1.0))
 628     {
 629       fHMuMohterPhiDifGen[1][0]->Fill(deltaphi(motherphi-mcphi));
 630       fHMuMohterEtaDifGen[1][0]->Fill(mothereta-mceta);
 631     }
 632     else if((1.0 <= mcpt) && (mcpt < 2.0))
 633     {
 634       fHMuMohterPhiDifGen[1][1]->Fill(deltaphi(motherphi-mcphi));
 635       fHMuMohterEtaDifGen[1][1]->Fill(mothereta-mceta);
 636     }
 637     else if((2.0 <= mcpt) && (mcpt < 4.0))
 638     {
 639       fHMuMohterPhiDifGen[1][2]->Fill(deltaphi(motherphi-mcphi));
 640       fHMuMohterEtaDifGen[1][2]->Fill(mothereta-mceta);
 641     }
 642     // reconstructed
 643     if((0.5 <= trackpt) && (trackpt < 1.0))
 644     {
 645       fHMuMohterPhiDifRec[1][0]->Fill(deltaphi(motherphi-trackphi));
 646       fHMuMohterEtaDifRec[1][0]->Fill(mothereta-tracketa);
 647     }
 648     else if((1.0 <= trackpt) && (trackpt < 2.0))
 649     {
 650       fHMuMohterPhiDifRec[1][1]->Fill(deltaphi(motherphi-trackphi));
 651       fHMuMohterEtaDifRec[1][1]->Fill(mothereta-tracketa);
 652     }
 653     else if((2.0 <= trackpt) && (trackpt < 4.0))
 654     {
 655       fHMuMohterPhiDifRec[1][2]->Fill(deltaphi(motherphi-trackphi));
 656       fHMuMohterEtaDifRec[1][2]->Fill(mothereta-tracketa);
 657     }
 658     if(fESD) fHMuDCA[1]->Fill(dcavalue);
 659   }
 660   else
 661   {
 662     fHMuMotherGenPt[3]->Fill(mcpt, motherpt);
 663     fHMuMotherRecPt[3]->Fill(trackpt, motherpt);
 664     fHMuMotherGenPhi[3]->Fill(mcphi, motherphi);
 665     fHMuMotherRecPhi[3]->Fill(trackphi, motherphi);
 666     fHMuMotherGenEta[3]->Fill(mceta, mothereta);
 667     fHMuMotherRecEta[3]->Fill(tracketa, mothereta);
 668     // generated
 669     if((0.5 <= mcpt) && (mcpt < 1.0))
 670     {
 671       fHMuMohterPhiDifGen[3][0]->Fill(deltaphi(motherphi-mcphi));
 672       fHMuMohterEtaDifGen[3][0]->Fill(mothereta-mceta);
 673     }
 674     else if((1.0 <= mcpt) && (mcpt < 2.0))
 675     {
 676       fHMuMohterPhiDifGen[3][1]->Fill(deltaphi(motherphi-mcphi));
 677       fHMuMohterEtaDifGen[3][1]->Fill(mothereta-mceta);
 678     }
 679     else if((2.0 <= mcpt) && (mcpt < 4.0))
 680     {
 681       fHMuMohterPhiDifGen[3][2]->Fill(deltaphi(motherphi-mcphi));
 682       fHMuMohterEtaDifGen[3][2]->Fill(mothereta-mceta);
 683     }
 684     // reconstructed
 685     if((0.5 <= trackpt) && (trackpt < 1.0))
 686     {
 687       fHMuMohterPhiDifRec[3][0]->Fill(deltaphi(motherphi-trackphi));
 688       fHMuMohterEtaDifRec[3][0]->Fill(mothereta-tracketa);
 689     }
 690     else if((1.0 <= trackpt) && (trackpt < 2.0))
 691     {
 692       fHMuMohterPhiDifRec[3][1]->Fill(deltaphi(motherphi-trackphi));
 693       fHMuMohterEtaDifRec[3][1]->Fill(mothereta-tracketa);
 694     }
 695     else if((2.0 <= trackpt) && (trackpt < 4.0))
 696     {
 697       fHMuMohterPhiDifRec[3][2]->Fill(deltaphi(motherphi-trackphi));
 698       fHMuMohterEtaDifRec[3][2]->Fill(mothereta-tracketa);
 699     }
 700     if(fESD) fHMuDCA[3]->Fill(dcavalue);
 701   }
 702 }
 703 //________________________________________________________________________
 704 Double_t AliMuonEffMC::deltaphi(Double_t phi)
 705 {
 706   if(phi < -1.0*TMath::Pi()) { return (phi + TMath::TwoPi()); }
 707   else if(phi > TMath::Pi()) { return (phi - TMath::TwoPi()); }
 708   else { return phi; }
 709 }