PWG/muon/ReadAOD_MCBranch.C

   1 void ReadAOD_MCBranch(char *filein="AliAODs.root"){
   2
   3 //=============================================================================
   4 // Macro to read AOD+MC Branch:
   5 // in this example the macro reads generated and reconstructed JPsi
   6 // Plots refer to all MC generated particles, to MC particles which have
   7 // been reconstructed and reconstructed tracks
   8 //
   9 // Each aod event has 3 MC particles (J/Psi, mu1, mu2)
  10 // To access MC muons:   mctrackall->IsPhysicalPrimary()
  11 // To access JPsi muons: mctrackall->IsPrimary() && !mctrackall->IsPhysicalPrimary()
  12
  13 // AliStack::IsPhysicalPrimary()
  14 // Test if a particle is a physical primary according to the following definition:
  15 // Particles produced in the collision including products of strong and
  16 // electromagnetic decay and excluding feed-down from weak decays of strange
  17 // particles.
  18 //=============================================================================
  19
  20   gStyle->SetFrameFillColor(0);
  21   gStyle->SetCanvasColor(0);
  22   gStyle->SetPalette(1);
  23
  24 // Open AOD file
  25   TFile *file = new TFile(filein);
  26   TTree *aodTree = file->Get("aodTree");
  27
  28 // Read AOD event
  29   AliAODEvent *event = new AliAODEvent();
  30   event->ReadFromTree(aodTree);
  31
  32 // Read MC info
  33   TClonesArray *mcarray = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
  34   if(!mcarray) {
  35     printf("No MC info in AOD!\n");
  36     return;
  37   }
  38
  39 // Book Histos
  40   TH2D *hpt_mc_rec = new TH2D("hpt_mc_rec","Pt correlations MC-REC",50,0.,10.,50,0.,10.);
  41   TH1D *hpt_rec = new TH1D("hpt_rec","Pt REC tracks",20,0.,10.);
  42   TH1D *hpt_mc = new TH1D("hpt_mc","PT MC particles corresponding to rec tracks",20,0.,10.);
  43   TH1D *hpt_allmc = new TH1D("hpt_allmc","PT all MC particles",20,0.,10.);
  44   TH1D *hpt_primarymc = new TH1D("hpt_primarymc","PT MC primaries (J/Psi and muons)",20,0.,10.);
  45   TH1D *hpt_physprimarymc = new TH1D("hpt_physprimarymc","PT Muons MC",20,0.,10.);
  46   TH1D *hpt_notphysprimarymc = new TH1D("hpt_notphysprimarymc","Pt JPsi MC",20,0.,10.);
  47   TH1D *hpt_diff = new TH1D("hpt_diff","Pt: Rec - MC",100,-5.,5.);
  48
  49   TH1D *hCharge_primarymc = new TH1D("hCharge_primarymc","Charge MC primary particle",10,-5.,5.);
  50   TH1D *hCharge_physprimmc = new TH1D("hCharge_physprimmc","Charge Muons MC",10,-5.,5.);
  51   TH1D *hCharge_notphysprimmc = new TH1D("hCharge_notphysprimmc","Charge JPsi MC",10,-5.,5.);
  52   TH1D *hCharge_rec = new TH1D("hCharge_rec","hCharge_rec",10,-5.,5.);
  53   TH1D *hCharge_mc = new TH1D("hCharge_mc","hCharge_mc",10,-5.,5.);
  54
  55   TH1D *hMass_mc = new TH1D("hMass_mc","Mass MC particles",100,0.,5..);
  56   TH1D *hMass_rec = new TH1D("hMass_rec","Mass rec tracks",100,0.,5..);
  57   TH1D *hMass_allmc = new TH1D("hMass_allmc","Mass all MC particles",100,0.,5..);
  58   TH1D *hMass_primarymc = new TH1D("hMass_primarymc","Mass primary particles MC",100,0.,5..);
  59   TH1D *hMass_physprimarymc = new TH1D("hMass_physprimarymc","Mass muons MC",100,0.,5..);
  60   TH1D *hMass_notphysprimarymc = new TH1D("hMass_notphysprimarymc","Mass JPsi MC",100,0.,5..);
  61   TH1D *hMassJPsi_mc = new TH1D("hMassJPsi_mc","Mass JPsi MC",100,0.,5..);
  62
  63 // Loop on AOD events
  64   Int_t nev=0;
  65
  66   while(aodTree->GetEvent(nev++)){
  67
  68 // Loop on all MC tracks
  69     for(int ii=0;ii<mcarray->GetEntries();ii++){
  70       AliAODMCParticle *mctrackall = (AliAODMCParticle*) mcarray->At(ii);
  71       // all particles
  72       hpt_allmc->Fill(mctrackall->Pt());
  73       hMass_allmc->Fill(mctrackall->M());
  74
  75       // Muons
  76       if(mctrackall->IsPhysicalPrimary()) {
  77         hCharge_physprimmc->Fill(mctrackall->Charge());
  78         hpt_physprimarymc->Fill(mctrackall->Pt());
  79         hMass_physprimarymc->Fill(mctrackall->M());
  80       }
  81
  82       if(mctrackall->IsPrimary()) {
  83         hpt_primarymc->Fill(mctrackall->Pt());
  84         hMass_primarymc->Fill(mctrackall->M());
  85         hCharge_primarymc->Fill(mctrackall->Charge());
  86       }
  87
  88       // JPsi
  89       if(mctrackall->IsPrimary() && !mctrackall->IsPhysicalPrimary()) {
  90         hCharge_notphysprimmc->Fill(mctrackall->Charge());
  91         hMass_notphysprimarymc->Fill(mctrackall->M());
  92         hpt_notphysprimarymc->Fill(mctrackall->Pt());
  93       }
  94     }
  95
  96 // Loop on reconstructed tracks
  97     for(int i=0;i<event->GetNumberOfTracks();i++){
  98       AliAODTrack *track = dynamic_cast<AliAODTrack*>(event->GetTrack(i));
  99       if(!track) {
 100         printf("Error: no track! \n");
 101         continue;
 102       }
 103       hpt_rec->Fill(track->Pt());
 104       hMass_rec->Fill(track->M());
 105       hCharge_rec->Fill(track->Charge());
 106
 107 // Read MC particle corresponding to reconstructed tracks
 108       AliAODMCParticle *mctrack = (AliAODMCParticle*) mcarray->At(track->GetLabel());
 109       if(!mctrack){
 110         printf("Error: no MC track! \n");
 111         continue;
 112       }
 113       //printf("nev=%d i=%d mc=%f rec=%f\n",nev-1,i,mctrack->Pt(),track->Pt());
 114       hpt_mc->Fill(mctrack->Pt());
 115       hpt_mc_rec->Fill(track->Pt(),mctrack->Pt());
 116       hpt_diff->Fill(track->Pt()-mctrack->Pt());
 117       hMass_mc->Fill(mctrack->M());
 118       hCharge_mc->Fill(mctrack->Charge());
 119    }
 120  }
 121
 122 // Plots
 123  TCanvas *c = new TCanvas("c","PT spectra",20,20,600,600);
 124  c->Divide(3,3);
 125  c->cd(1);
 126  hpt_primarymc->Draw();
 127  hpt_primarymc->GetXaxis()->SetTitle("Pt MC primary particles");
 128  c->cd(2);
 129  hpt_mc->GetXaxis()->SetTitle("Pt MC particles");
 130  hpt_mc->Draw();
 131  hpt_mc->SetLineColor(4);
 132  c->cd(3);
 133  hpt_rec->GetXaxis()->SetTitle("Rec Pt track");
 134  hpt_rec->Draw();
 135  hpt_rec->SetLineColor(2);
 136  c->cd(4);
 137  hpt_physprimarymc->Draw();
 138  hpt_physprimarymc->GetXaxis()->SetTitle("Pt MC muon track");
 139  c->cd(5);
 140  hpt_notphysprimarymc->Draw();
 141  hpt_notphysprimarymc->GetXaxis()->SetTitle("Pt MC JPsi track");
 142  c->cd(7);
 143  hpt_mc_rec->Draw("colz");
 144  hpt_mc_rec->GetXaxis()->SetTitle("Pt rec track");
 145  hpt_mc_rec->GetYaxis()->SetTitle("Pt MC track");
 146  c->cd(8);
 147  hpt_diff->Draw();
 148  hpt_diff->GetXaxis()->SetTitle("Rec-MC Pt");
 149
 150  TCanvas *c1 = new TCanvas("c1","Mass distributions",40,40,600,600);
 151  c1->Divide(2,2);
 152  c1->cd(1);
 153  hMass_primarymc->Draw();
 154  hMass_primarymc->GetXaxis()->SetTitle("Mass MC primary particles");
 155  c1->cd(2);
 156  hMass_physprimarymc->Draw();
 157  hMass_physprimarymc->GetXaxis()->SetTitle("Mass MC muon particles");
 158  c1->cd(3);
 159  hMass_notphysprimarymc->Draw();
 160  hMass_notphysprimarymc->GetXaxis()->SetTitle("Mass MC J/Psi");
 161  c1->cd(4);
 162  hMass_rec->Draw();
 163  hMass_rec->GetXaxis()->SetTitle("Mass Rec tracks");
 164
 165  TCanvas *c2= new TCanvas("c2","Charge distributions",60,60,600,600);
 166  c2->Divide(2,2);
 167  c2->cd(1);
 168  hCharge_primarymc->Draw();
 169  hCharge_primarymc->GetXaxis()->SetTitle("Charge MC primary particles");
 170  c2->cd(2);
 171  hCharge_physprimmc->Draw();
 172  hCharge_physprimmc->GetXaxis()->SetTitle("Charge muon MC");;
 173  c2->cd(3);
 174  hCharge_notphysprimmc->Draw();
 175  hCharge_notphysprimmc->GetXaxis()->SetTitle("Charge J/Psi MC");;
 176  c2->cd(4);
 177  hCharge_rec->Draw();
 178  hCharge_rec->GetXaxis()->SetTitle("Charge rec tracks");
 179 }