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 "AliAnalysisManager.h"
  15 #include "AliAnalysisTask.h"
  16 #include "AliESDVertex.h"
  17 #include "AliESDEvent.h"
  18 #include "AliAODEvent.h"
  19 #include "AliESDMuonTrack.h"
  20 #include "AliESDVertex.h"
  21 #include "AliAODVertex.h"
  22 #include "AliCentrality.h"
  23 #include "AliVParticle.h"
  24 #include "AliMCParticle.h"
  25 #include "AliAnalysisHelperJetTasks.h"
  26 #include "AliMCEvent.h"
  27
  28 using std::cout;
  29 using std::endl;
  30
  31 ClassImp(AliMuonEffMC)
  32
  33 //________________________________________________________________________
  34 AliMuonEffMC::AliMuonEffMC() :
  35   AliAnalysisTaskSE(), fESD(0), fAOD(0), fMC(0), fCentrality(99), fZVertex(99), fOutputList(0x0),
  36   fHEventStat(0), fHEvt(0x0), fIsMc(kTRUE), fMDProcess(kTRUE), fCentralityEstimator("V0M"),
  37   fNEtaBins(15), fNpTBins(100), fNCentBins(5), fNZvtxBins(10), fNPhiBins(24),
  38   fHMuonParGen(0x0), fHMuonDetGen(0x0), fHMuonDetRec(0x0), fHEtcDetRec(0x0)
  39 {
  40   // Constructor
  41   DefineInput(0, TChain::Class());
  42   DefineOutput(1, TList::Class());
  43
  44   for(Int_t i=0; i<4; i++)
  45   {
  46     fHMuMotherGenPt[i] = 0x0;
  47     fHMuMotherRecPt[i] = 0x0;
  48     fHMuMotherGenPhi[i] = 0x0;
  49     fHMuMotherRecPhi[i] = 0x0;
  50     fHMuMotherGenEta[i] = 0x0;
  51     fHMuMotherRecEta[i] = 0x0;
  52     fHMuDCA[i] = 0x0;
  53   }
  54 }
  55
  56 //________________________________________________________________________
  57 AliMuonEffMC::AliMuonEffMC(const char *name) :
  58   AliAnalysisTaskSE(name), fESD(0), fAOD(0), fMC(0), fCentrality(99), fZVertex(99), fOutputList(0x0),
  59   fHEventStat(0), fHEvt(0x0),  fIsMc(kTRUE), fMDProcess(kTRUE), fCentralityEstimator("V0M"),
  60   fNEtaBins(15), fNpTBins(100), fNCentBins(5), fNZvtxBins(10), fNPhiBins(24),
  61   fHMuonParGen(0x0), fHMuonDetGen(0x0), fHMuonDetRec(0x0), fHEtcDetRec(0x0)
  62 {
  63   // Constructor
  64   DefineInput(0, TChain::Class());
  65   DefineOutput(1, TList::Class());
  66
  67   for(Int_t i=0; i<4; i++)
  68   {
  69     fHMuMotherGenPt[i] = 0x0;
  70     fHMuMotherRecPt[i] = 0x0;
  71     fHMuMotherGenPhi[i] = 0x0;
  72     fHMuMotherRecPhi[i] = 0x0;
  73     fHMuMotherGenEta[i] = 0x0;
  74     fHMuMotherRecEta[i] = 0x0;
  75     fHMuDCA[i] = 0x0;
  76   }
  77 }
  78
  79 //________________________________________________________________________
  80 AliMuonEffMC::~AliMuonEffMC()
  81 {
  82   //Destructor
  83   if(fOutputList) delete fOutputList;
  84 }
  85
  86 //________________________________________________________________________
  87 void AliMuonEffMC::UserCreateOutputObjects()
  88 {
  89   // Create histograms
  90   // Called once (per slave on PROOF!)
  91
  92   fOutputList = new TList();
  93   fOutputList->SetOwner(1);
  94
  95   fHEventStat = new TH1D("fHEventStat","Event statistics for analysis",18,0,18);
  96   fHEventStat->GetXaxis()->SetBinLabel(1,"Event");
  97   fHEventStat->GetXaxis()->SetBinLabel(2,"SelectedEvent");
  98   fHEventStat->GetXaxis()->SetBinLabel(3,"File");
  99   fHEventStat->GetXaxis()->SetBinLabel(4,"fSPHighpt");  //!Global Trigger Single plus High p_T
 100   fHEventStat->GetXaxis()->SetBinLabel(5,"fSPAllpt");   //!Global Trigger Single plus All p_T
 101   fHEventStat->GetXaxis()->SetBinLabel(6,"fSMLowpt");   //!Global Trigger Single minus Low p_T
 102   fHEventStat->GetXaxis()->SetBinLabel(7,"fSMHighpt");  //!Global Trigger Single minus High p_T
 103   fHEventStat->GetXaxis()->SetBinLabel(8,"fSMAllpt");   //!Global Trigger Single minus All p_T
 104   fHEventStat->GetXaxis()->SetBinLabel(9,"fSULowpt");   //!Global Trigger Single undefined Low p_T
 105   fHEventStat->GetXaxis()->SetBinLabel(10,"fSUHighpt"); //!Global Trigger Single undefined High p_T
 106   fHEventStat->GetXaxis()->SetBinLabel(11,"fSUAllpt");  //!Global Trigger Single undefined All p_T
 107   fHEventStat->GetXaxis()->SetBinLabel(12,"fUSLowpt");  //!Global Trigger UnlikeSign pair Low p_T
 108   fHEventStat->GetXaxis()->SetBinLabel(13,"fUSHighpt"); //!Global Trigger UnlikeSign pair High p_T
 109   fHEventStat->GetXaxis()->SetBinLabel(14,"fUSAllpt");  //!Global Trigger UnlikeSign pair All p_T
 110   fHEventStat->GetXaxis()->SetBinLabel(15,"fLSLowpt");  //!Global Trigger LikeSign pair pair Low p_T
 111   fHEventStat->GetXaxis()->SetBinLabel(16,"fLSHighpt"); //!Global Trigger LikeSign pair pair High p_T
 112   fHEventStat->GetXaxis()->SetBinLabel(17,"fLSAllpt");  //!Global Trigger LikeSign pair pair All p_T
 113   fHEventStat->GetXaxis()->SetBinLabel(18,"fSPLowpt");   //!Global Trigger Single plus Low p_T
 114   fOutputList->Add(fHEventStat);
 115
 116   fHEvt = new TH2F("fHEvt", "Event-level variables; Zvtx; Cent", 30, -15, 15, 103, -2, 101);
 117   fOutputList->Add(fHEvt);
 118
 119
 120   Int_t iTrackBin[5];
 121   Double_t* trackBins[5];
 122   const char* trackAxisTitle[5];
 123
 124   // eta
 125   Double_t etaBins[fNEtaBins+1];
 126   for (Int_t i=0; i<=fNEtaBins; i++) { etaBins[i] = (Double_t)(-4.0 + 1.5/fNEtaBins*i); }
 127   iTrackBin[0] = fNEtaBins;
 128   trackBins[0] = etaBins;
 129   trackAxisTitle[0] = "#eta";
 130
 131   // p_T
 132   Double_t pTBins[fNpTBins+1];
 133   for (Int_t i=0; i<=fNpTBins; i++) { pTBins[i] = (Double_t)(10.0/fNpTBins*i); }
 134   iTrackBin[1] = fNpTBins;
 135   trackBins[1] = pTBins;
 136   trackAxisTitle[1] = "p_{T} (GeV/c)";
 137
 138   // centrality
 139   Double_t CentBins[fNCentBins+1];
 140   for (Int_t i=0; i<=fNCentBins; i++) { CentBins[i] = (Double_t)(100.0/fNCentBins*i); }
 141   iTrackBin[2] = fNCentBins;
 142   trackBins[2] = CentBins;
 143   trackAxisTitle[2] = "Cent";
 144
 145   // Z-vertex
 146   Double_t ZvtxBins[fNZvtxBins+1];
 147   for (Int_t i=0; i<=fNZvtxBins; i++) { ZvtxBins[i] = (Double_t)(-10.0 + 20.0/fNZvtxBins*i); }
 148   iTrackBin[3] = fNZvtxBins;
 149   trackBins[3] = ZvtxBins;
 150   trackAxisTitle[3] = "Zvtx";
 151
 152   // phi
 153   Double_t phiBins[fNPhiBins+1];
 154   for (Int_t i=0; i<=fNPhiBins; i++) { phiBins[i] = (Double_t)(TMath::TwoPi()/fNPhiBins * i); }
 155   iTrackBin[4] = fNPhiBins;
 156   trackBins[4] = phiBins;
 157   trackAxisTitle[4] = "#phi";
 158
 159   // THn for tracking efficiency
 160   fHMuonParGen = new THnF("fHMuonParGen", "", 5, iTrackBin, 0, 0);
 161   for (Int_t i=0; i<5; i++)
 162   {
 163     fHMuonParGen->SetBinEdges(i, trackBins[i]);
 164     fHMuonParGen->GetAxis(i)->SetTitle(trackAxisTitle[i]);
 165   }
 166   fHMuonParGen->Sumw2();
 167   fOutputList->Add(fHMuonParGen);
 168
 169   fHMuonDetGen = (THnF*) fHMuonParGen->Clone("fHMuonDetGen");
 170   fHMuonDetGen->Sumw2();
 171   fOutputList->Add(fHMuonDetGen);
 172
 173   fHMuonDetRec = (THnF*) fHMuonParGen->Clone("fHMuonDetRec");
 174   fHMuonDetRec->Sumw2();
 175   fOutputList->Add(fHMuonDetRec);
 176
 177   fHEtcDetRec = (THnF*) fHMuonParGen->Clone("fHEtcDetRec");
 178   fHEtcDetRec->Sumw2();
 179   fOutputList->Add(fHEtcDetRec);
 180
 181   if(fMDProcess)
 182   {
 183     const char* MotherSpecies[4] = {"Pion","Kaon","D", "Etc"};
 184
 185     for(Int_t i=0; i<4; i++)
 186     {
 187       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);
 188       fOutputList->Add(fHMuMotherGenPt[i]);
 189
 190       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);
 191       fOutputList->Add(fHMuMotherRecPt[i]);
 192
 193       fHMuMotherGenPhi[i] = new TH2F(Form("fHMuMotherGenPhi_%s",MotherSpecies[i]),";#phi_{gen};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
 194       fOutputList->Add(fHMuMotherGenPhi[i]);
 195
 196       fHMuMotherRecPhi[i] = new TH2F(Form("fHMuMotherRecPhi_%s",MotherSpecies[i]),";#phi_{rec};mother #phi;",100, 0, TMath::TwoPi(), 100, 0, TMath::TwoPi());
 197       fOutputList->Add(fHMuMotherRecPhi[i]);
 198
 199       fHMuMotherGenEta[i] = new TH2F(Form("fHMuMotherGenEta_%s",MotherSpecies[i]),";#eta_{gen};mother #eta;",100, -5., -1., 100, -5., -1.);
 200       fOutputList->Add(fHMuMotherGenEta[i]);
 201
 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
 205       fHMuDCA[i] =  new TH1F(Form("fHMuDCA_%s",MotherSpecies[i]), ";DCA", 100, 0, 50);
 206       fOutputList->Add(fHMuDCA[i]);
 207     }
 208   }
 209   PostData(1, fOutputList);
 210 }
 211
 212 //________________________________________________________________________
 213 void AliMuonEffMC::UserExec(Option_t *)
 214 {
 215   // Main loop, Called for each event
 216   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 217   if (!fESD)
 218   {
 219     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 220   }
 221
 222   if (!fESD && !fAOD)
 223   {
 224     AliError("Neither fESD nor fAOD available");
 225     return;
 226   }
 227
 228   fHEventStat->Fill(0.5);
 229   if(fIsMc)
 230   {
 231     fMC = MCEvent();
 232     if (!fMC) {
 233       printf("ERROR: fMC not available\n");
 234       return;
 235     }
 236   }
 237
 238   const AliESDVertex* vertex = fESD->GetPrimaryVertex();
 239   fZVertex = vertex->GetZ();
 240   if(fESD->GetCentrality()) fCentrality = fESD->GetCentrality()->GetCentralityPercentile("fCentralityEstimator");
 241
 242   // Fill Event histogram
 243   fHEvt->Fill(fZVertex, fCentrality);
 244
 245   Int_t iVerb = 0;
 246
 247    // Centrality, vertex, other event variables...
 248   if (!VertexOk(fESD)) {
 249     if (iVerb>1)
 250       AliInfo(Form("Event REJECTED. z = %.1f", fZVertex));
 251     return;
 252   }
 253
 254   if (fCentrality > 100. || fCentrality < -1.5) {
 255     if (iVerb>1)
 256       AliInfo(Form("Event REJECTED. fCentrality = %.1f", fCentrality));
 257     return;
 258   }
 259
 260   if(fCentrality < 0) fCentrality = 1.0;
 261   fHEventStat->Fill(1.5);
 262
 263   ULong64_t trigword=fESD->GetTriggerMask();
 264
 265   if (trigword & 0x01) fHEventStat->Fill(17.5);
 266   if (trigword & 0x02) fHEventStat->Fill(3.5);
 267   if (trigword & 0x04) fHEventStat->Fill(4.5);
 268   if (trigword & 0x08) fHEventStat->Fill(5.5);
 269   if (trigword & 0x010) fHEventStat->Fill(6.5);
 270   if (trigword & 0x020) fHEventStat->Fill(7.5);
 271   if (trigword & 0x040) fHEventStat->Fill(8.5);
 272   if (trigword & 0x080) fHEventStat->Fill(9.5);
 273   if (trigword & 0x100) fHEventStat->Fill(10.5);
 274   if (trigword & 0x200) fHEventStat->Fill(11.5);
 275   if (trigword & 0x400) fHEventStat->Fill(12.5);
 276   if (trigword & 0x800) fHEventStat->Fill(13.5);
 277   if (trigword & 0x1000) fHEventStat->Fill(14.5);
 278   if (trigword & 0x2000) fHEventStat->Fill(15.5);
 279   if (trigword & 0x4000) fHEventStat->Fill(16.5);
 280
 281   if(fIsMc)
 282   {
 283     for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfMuonTracks(); iTrack++)
 284     {
 285       // loop on muon tracks
 286       AliESDMuonTrack* muonTrack = fESD->GetMuonTrack(iTrack);
 287       if(muonTrack)
 288       {
 289         if(!IsGoodMUONtrack(*muonTrack)) continue;
 290         Double_t fillArrayDetRec[5] = { muonTrack->Eta(), muonTrack->Pt(), fCentrality, fZVertex, muonTrack->Phi() };
 291         fHMuonDetRec->Fill(fillArrayDetRec);
 292
 293         Int_t label = TMath::Abs(muonTrack->GetLabel());
 294         AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(label);
 295         if(TMath::Abs(McParticle->PdgCode()) != 13)
 296         {
 297           fHEtcDetRec->Fill(fillArrayDetRec);
 298           continue;
 299         }
 300         Double_t fillArrayDetGen[5] = { McParticle->Eta(), McParticle->Pt(), fCentrality, fZVertex, McParticle->Phi() };
 301         fHMuonDetGen->Fill(fillArrayDetGen);
 302
 303         if(fMDProcess)
 304         {
 305           if(McParticle->GetMother() > 0)
 306           {
 307             AliMCParticle *MotherParticle  = (AliMCParticle*)fMC->GetTrack(McParticle->GetMother());
 308             Int_t motherlabel = TMath::Abs(MotherParticle->PdgCode());
 309
 310             if(motherlabel==411 || motherlabel==413 || motherlabel==421 || motherlabel==423 || motherlabel==431 || motherlabel==433 || motherlabel==10413 || motherlabel==10411 || motherlabel==10423 || motherlabel==10421 || motherlabel==10433 || motherlabel==10431 || motherlabel==20413 || motherlabel==415 || motherlabel==20423 || motherlabel==425 || motherlabel==20433 || motherlabel==435)
 311             {
 312               fHMuMotherGenPt[2]->Fill(McParticle->Pt(), MotherParticle->Pt());
 313               fHMuMotherRecPt[2]->Fill(muonTrack->Pt(), MotherParticle->Pt());
 314               fHMuMotherGenPhi[2]->Fill(McParticle->Phi(), MotherParticle->Phi());
 315               fHMuMotherRecPhi[2]->Fill(muonTrack->Phi(), MotherParticle->Phi());
 316               fHMuMotherGenEta[2]->Fill(McParticle->Eta(), MotherParticle->Eta());
 317               fHMuMotherRecEta[2]->Fill(muonTrack->Eta(), MotherParticle->Eta());
 318               fHMuDCA[2]->Fill(muonTrack->GetDCA());
 319             }
 320             else if(motherlabel==211)
 321             {
 322               fHMuMotherGenPt[0]->Fill(McParticle->Pt(), MotherParticle->Pt());
 323               fHMuMotherRecPt[0]->Fill(muonTrack->Pt(), MotherParticle->Pt());
 324               fHMuMotherGenPhi[0]->Fill(McParticle->Phi(), MotherParticle->Phi());
 325               fHMuMotherRecPhi[0]->Fill(muonTrack->Phi(), MotherParticle->Phi());
 326               fHMuMotherGenEta[0]->Fill(McParticle->Eta(), MotherParticle->Eta());
 327               fHMuMotherRecEta[0]->Fill(muonTrack->Eta(), MotherParticle->Eta());
 328               fHMuDCA[0]->Fill(muonTrack->GetDCA());
 329             }
 330             else if(motherlabel==321)
 331             {
 332               fHMuMotherGenPt[1]->Fill(McParticle->Pt(), MotherParticle->Pt());
 333               fHMuMotherRecPt[1]->Fill(muonTrack->Pt(), MotherParticle->Pt());
 334               fHMuMotherGenPhi[1]->Fill(McParticle->Phi(), MotherParticle->Phi());
 335               fHMuMotherRecPhi[1]->Fill(muonTrack->Phi(), MotherParticle->Phi());
 336               fHMuMotherGenEta[1]->Fill(McParticle->Eta(), MotherParticle->Eta());
 337               fHMuMotherRecEta[1]->Fill(muonTrack->Eta(), MotherParticle->Eta());
 338               fHMuDCA[1]->Fill(muonTrack->GetDCA());
 339             }
 340             else
 341             {
 342               fHMuMotherGenPt[3]->Fill(McParticle->Pt(), MotherParticle->Pt());
 343               fHMuMotherRecPt[3]->Fill(muonTrack->Pt(), MotherParticle->Pt());
 344               fHMuMotherGenPhi[3]->Fill(McParticle->Phi(), MotherParticle->Phi());
 345               fHMuMotherRecPhi[3]->Fill(muonTrack->Phi(), MotherParticle->Phi());
 346               fHMuMotherGenEta[3]->Fill(McParticle->Eta(), MotherParticle->Eta());
 347               fHMuMotherRecEta[3]->Fill(muonTrack->Eta(), MotherParticle->Eta());
 348               fHMuDCA[3]->Fill(muonTrack->GetDCA());
 349             }
 350           }
 351         } // (mother hadron) : (daughter muon) QA
 352       }
 353     }
 354
 355     for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)
 356     {
 357       AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(ipart);
 358       if(TMath::Abs(McParticle->PdgCode())==13)
 359       {
 360         Double_t fillArrayParGen[5] = { McParticle->Eta(), McParticle->Pt(), fCentrality, fZVertex, McParticle->Phi() };
 361         fHMuonParGen->Fill(fillArrayParGen);
 362       }
 363     }
 364   }
 365   PostData(1, fOutputList);
 366   return;
 367 }
 368
 369 //________________________________________________________________________
 370 void AliMuonEffMC::Terminate(Option_t *)
 371 {
 372   // Draw result to the screen
 373   // Called once at the end of the query
 374
 375   fOutputList = dynamic_cast<TList*> (GetOutputData(1));
 376   if (!fOutputList) {
 377     AliError("Output list not available");
 378     return;
 379   }
 380 }
 381
 382 //________________________________________________________________________
 383 Bool_t AliMuonEffMC::VertexOk(TObject* obj) const
 384 {
 385   // Modified from AliAnalyseLeadingTrackUE::VertexSelection()
 386
 387   Int_t nContributors  = 0;
 388   Double_t zVertex     = 999;
 389   TString name("");
 390
 391   if (obj->InheritsFrom("AliESDEvent")) {
 392     AliESDEvent* esdevt = (AliESDEvent*) obj;
 393     const AliESDVertex* vtx = esdevt->GetPrimaryVertex();
 394     if (!vtx)
 395       return 0;
 396     nContributors = vtx->GetNContributors();
 397     zVertex       = vtx->GetZ();
 398     name          = vtx->GetName();
 399   }
 400   else if (obj->InheritsFrom("AliAODEvent")) {
 401     AliAODEvent* aodevt = (AliAODEvent*) obj;
 402     if (aodevt->GetNumberOfVertices() < 1)
 403       return 0;
 404     const AliAODVertex* vtx = aodevt->GetPrimaryVertex();
 405     nContributors = vtx->GetNContributors();
 406     zVertex       = vtx->GetZ();
 407     name          = vtx->GetName();
 408   }
 409
 410   // Reject if TPC-only vertex
 411   if (name.CompareTo("TPCVertex")==0)
 412     return kFALSE;
 413
 414   // Check # contributors and range...
 415   if( nContributors < 1 || TMath::Abs(zVertex) > 10 ) {
 416     return kFALSE;
 417   }
 418
 419   return kTRUE;
 420 }
 421
 422
 423 //________________________________________________________________________
 424 Bool_t AliMuonEffMC::IsGoodMUONtrack(AliESDMuonTrack &track)
 425 {
 426   // Applying track cuts for MUON tracks
 427   if(!track.ContainTrackerData()) return kFALSE;
 428   if(!track.ContainTriggerData()) return kFALSE;
 429
 430   Double_t thetaTrackAbsEnd = TMath::ATan(track.GetRAtAbsorberEnd()/505.) * TMath::RadToDeg();
 431   Double_t eta = track.Eta();
 432
 433   // Theta cut at absorber end
 434   if(thetaTrackAbsEnd <= 2. || thetaTrackAbsEnd >= 10.) return kFALSE;
 435   // Eta cut
 436   if(eta <= -4. || eta >= -2.5) return kFALSE;
 437   if(track.GetMatchTrigger() <= 0) return kFALSE;
 438
 439   return kTRUE;
 440 }