PWG/muon/AliAnalysisTaskAODvsESD.cxx

   1
   2 /* $Id$ */
   3
   4 #include "TChain.h"
   5 #include "TTree.h"
   6 #include "TNtuple.h"
   7 #include "TBranch.h"
   8 #include "TH1F.h"
   9 #include "TCanvas.h"
  10 #include <TLorentzVector.h>
  11 #include <TVector3.h>
  12
  13 #include "AliAnalysisTask.h"
  14 #include "AliAnalysisManager.h"
  15 #include "AliESDInputHandler.h"
  16 #include "AliESDEvent.h"
  17 #include "AliESDMuonTrack.h"
  18 #include "AliAODInputHandler.h"
  19 #include "AliAODHandler.h"
  20 #include "AliAODEvent.h"
  21 #include "AliAnalysisTaskAODvsESD.h"
  22
  23 ClassImp(AliAnalysisTaskAODvsESD)
  24
  25 //________________________________________________________________________
  26 AliAnalysisTaskAODvsESD::AliAnalysisTaskAODvsESD(const char *name)
  27 : AliAnalysisTask(name, ""), fESD(0), fAOD(0), fList(0), fMuonNtuple(0), fMuonNtupleAOD(0), fInvMass(0), fInvMassAOD(0)
  28 {
  29   // Constructor
  30   // Define input and output slots here
  31   // Input slot #0 works with a TChain
  32   DefineInput(0, TChain::Class());
  33   // Output slot #0 writes NTuple/histos into a TList
  34   DefineOutput(0, TList::Class());
  35 }
  36
  37 //________________________________________________________________________
  38 void AliAnalysisTaskAODvsESD::ConnectInputData(Option_t *)
  39 {
  40   // Connect ESD here
  41   TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
  42   if (!esdTree) {
  43     Printf("ERROR: Could not read chain from input slot 0");
  44   } else {
  45     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  46     if (!esdH) {
  47       Printf("ERROR: Could not get ESDInputHandler");
  48     } else
  49       fESD = esdH->GetEvent();
  50   }
  51
  52   // Connect AOD here
  53   AliAODHandler *aodH = dynamic_cast<AliAODHandler*> (AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
  54   if (!aodH) {
  55     Printf("ERROR: Could not get AODHandler");
  56   } else
  57     fAOD = aodH->GetAOD();
  58 }
  59
  60 //________________________________________________________________________
  61 void AliAnalysisTaskAODvsESD::CreateOutputObjects()
  62 {
  63   // This method has to be called INSIDE the user redefined CreateOutputObjects
  64   // method, before creating each object corresponding to the output containers
  65   // that are to be written to a file. This need to be done in general for the big output
  66   // objects that may not fit memory during processing.
  67   OpenFile(0);
  68
  69   // Create Ntuples
  70   // For single muons: thetaX, thetaY, ptInv, eta, phi, theta, px, py, pz,ch
  71   fMuonNtuple = new TNtuple("fMuonNtuple","Muon information","thX:thY:ptI:eta:phi:theta:px:py:pz");
  72   fMuonNtupleAOD = new TNtuple("fMuonNtupleAOD","Muon information","etaAOD:phiAOD:thetaAOD:pxAOD:pyAOD:pzAOD");
  73   // Create histos for inv mass
  74   fInvMass = new TH1F("fInvMass","Inv. mass from ESDs",140,0,7);
  75   fInvMassAOD = new TH1F("fInvMassAOD","Inv. mass from AOD",140,0,7);
  76
  77   // Add Ntuples to the list
  78   fList = new TList();
  79   fList->Add(fMuonNtuple);
  80   fList->Add(fMuonNtupleAOD);
  81   fList->Add(fInvMass);
  82   fList->Add(fInvMassAOD);
  83 }
  84
  85 //________________________________________________________________________
  86 void AliAnalysisTaskAODvsESD::Exec(Option_t *)
  87 {
  88   // Main loop, called for each event
  89   if (!fESD) {
  90     Printf("ERROR: fESD not available");
  91     return;
  92   }
  93   if (!fAOD) {
  94     Printf("ERROR: fAOD not available");
  95     return;
  96   }
  97
  98   // ESDs
  99   // for muons
 100   Float_t muonMass = 0.105658369;
 101   Float_t thetaX1=0;    Float_t thetaY1=0;
 102   Float_t ptInv1=0;     Float_t pYZ1=0;
 103   Float_t pxM1=0;       Float_t pyM1=0;        Float_t pzM1=0;
 104   Float_t etaM1=0;      Float_t phiM1=0;       Float_t thetaM1=0;
 105   Int_t chargeM1=0;      Float_t energyM1=0;    TLorentzVector lvM1;
 106   Float_t thetaX2=0;    Float_t thetaY2=0;
 107   Float_t ptInv2=0;     Float_t pYZ2=0;
 108   Float_t pxM2=0;       Float_t pyM2=0;        Float_t pzM2=0;
 109   Int_t chargeM2=0;      Float_t energyM2=0;    TLorentzVector lvM2;
 110   //for J/psi
 111   Float_t invM=0;
 112   TLorentzVector lvJpsi;
 113
 114   int nmt = fESD->GetNumberOfMuonTracks();
 115   // first loop over muon tracks in event
 116   for (Int_t mTrack = 0; mTrack < nmt; mTrack++) {
 117     thetaX1=0; thetaY1=0; ptInv1=0; pYZ1=0; chargeM1=0;
 118     pxM1=0;    pyM1=0;    pzM1=0;   energyM1=0;
 119     AliESDMuonTrack* muonTrack1 = fESD->GetMuonTrack(mTrack);
 120     thetaX1 = muonTrack1->GetThetaX();
 121     thetaY1 = muonTrack1->GetThetaY();
 122     ptInv1 = TMath::Abs(muonTrack1->GetInverseBendingMomentum());
 123     pYZ1 = 1/ptInv1;
 124     pzM1  = - pYZ1 / TMath::Sqrt(1.0 + TMath::Tan(thetaY1)*TMath::Tan(thetaY1));
 125     pxM1  = pzM1 * TMath::Tan(thetaX1);
 126     pyM1  = pzM1 * TMath::Tan(thetaY1);
 127     energyM1 = TMath::Sqrt(muonMass*muonMass + pxM1*pxM1 + pyM1*pyM1 + pzM1*pzM1);
 128     lvM1.SetPxPyPzE(pxM1,pyM1,pzM1,energyM1);
 129     chargeM1 = Int_t(TMath::Sign(1.,muonTrack1->GetInverseBendingMomentum()));
 130     etaM1 = muonTrack1->Eta();
 131     phiM1 = muonTrack1->Phi();
 132     thetaM1 = muonTrack1->Theta();
 133     fMuonNtuple->Fill(thetaX1,thetaY1,ptInv1,etaM1,phiM1,thetaM1,pxM1,pyM1,pzM1);
 134     // second loop over muon tracks in event
 135     for (Int_t mTrack2 = mTrack+1; mTrack2 < nmt; mTrack2++) {
 136       thetaX2=0; thetaY2=0; ptInv2=0; pYZ2=0; chargeM2=0;
 137       pxM2=0; pyM2=0; pzM2=0; energyM2=0;
 138       AliESDMuonTrack* muonTrack2 = fESD->GetMuonTrack(mTrack2);
 139       thetaX2 = muonTrack2->GetThetaX();
 140       thetaY2 = muonTrack2->GetThetaY();
 141       ptInv2 = TMath::Abs(muonTrack2->GetInverseBendingMomentum());
 142       pYZ2 = 1/ptInv2;
 143       pzM2  = - pYZ2 / TMath::Sqrt(1.0 + TMath::Tan(thetaY2)*TMath::Tan(thetaY2));
 144       pxM2  = pzM2 * TMath::Tan(thetaX2);
 145       pyM2  = pzM2 * TMath::Tan(thetaY2);
 146       chargeM2 = Int_t(TMath::Sign(1.,muonTrack2->GetInverseBendingMomentum()));
 147       energyM2 = TMath::Sqrt(muonMass*muonMass + pxM2*pxM2 + pyM2*pyM2 + pzM2*pzM2);
 148       // if muons have opposite charge
 149       if(chargeM1*chargeM2 == -1){
 150         lvM2.SetPxPyPzE(pxM2, pyM2, pzM2,energyM2);
 151         lvJpsi = lvM1 + lvM2;
 152         invM = lvJpsi.M();
 153         fInvMass->Fill(invM);
 154       } // end if muons with opposite charge
 155     } // end second loop over muon tracks in event
 156   } // end first loop over muon tracks in event
 157
 158
 159   // Created AOD
 160   Float_t pxAodM1=0;   Float_t pyAodM1=0;     Float_t pzAodM1=0;
 161   Float_t etaAodM1=0;  Float_t phiAodM1=0;    Float_t thetaAodM1=0;
 162   Int_t chargeAodM1=0;  Float_t energyAodM1=0; TLorentzVector lvAodM1;
 163   Float_t pxAodM2=0;   Float_t pyAodM2=0;     Float_t pzAodM2=0;
 164   Float_t etaAodM2=0;  Float_t phiAodM2=0;    Float_t thetaAodM2=0;
 165   Int_t chargeAodM2=0;  Float_t energyAodM2=0; TLorentzVector lvAodM2;
 166   //for J/psi
 167   Float_t invMAOD=0;
 168   TLorentzVector lvJpsiAOD;
 169
 170   int nmtAOD = fAOD->GetNumberOfTracks();
 171   // first loop over tracks
 172   for (Int_t mTrack = 0; mTrack < nmtAOD; mTrack++) {
 173     pxM1=0; pyM1=0; pzM1=0;
 174     AliAODTrack* muonTrack1 = fAOD->GetTrack(mTrack);
 175     if(muonTrack1->IsMuonTrack()){
 176       etaAodM1 = muonTrack1->Eta();
 177       phiAodM1 = muonTrack1->Phi();
 178       thetaAodM1 = muonTrack1->Theta();
 179       pxAodM1 = muonTrack1->Px();
 180       pyAodM1 = muonTrack1->Py();
 181       pzAodM1 = muonTrack1->Pz();
 182       chargeAodM1 = muonTrack1->Charge();
 183       energyAodM1 = TMath::Sqrt(muonMass*muonMass + pxAodM1*pxAodM1 + pyAodM1*pyAodM1 + pzAodM1*pzAodM1);
 184       lvAodM1.SetPxPyPzE(pxAodM1, pyAodM1, pzAodM1,energyAodM1);
 185       fMuonNtupleAOD->Fill(etaAodM1,phiAodM1,thetaAodM1,pxAodM1,pyAodM1,pzAodM1);
 186     }
 187     // second loop over tracks in event
 188     for (Int_t mTrack2 = mTrack+1; mTrack2 < nmtAOD; mTrack2++) {
 189       chargeAodM2=0;
 190       pxAodM2=0; pyAodM2=0; pzAodM2=0; energyAodM2=0;
 191       AliAODTrack* muonTrack2 = fAOD->GetTrack(mTrack2);
 192       if(muonTrack2->IsMuonTrack()){
 193         etaAodM2 = muonTrack2->Eta();
 194         phiAodM2 = muonTrack2->Phi();
 195         thetaAodM2 = muonTrack2->Theta();
 196         pxAodM2 = muonTrack2->Px();
 197         pyAodM2 = muonTrack2->Py();
 198         pzAodM2 = muonTrack2->Pz();
 199         energyAodM2 = TMath::Sqrt(muonMass*muonMass + pxAodM2*pxAodM2 + pyAodM2*pyAodM2 + pzAodM2*pzAodM2);
 200         chargeAodM2 = muonTrack2->Charge();
 201         // if muons of opposite charge
 202         if(chargeAodM1*chargeAodM2 == -1){
 203           lvAodM2.SetPxPyPzE(pxAodM2, pyAodM2, pzAodM2,energyAodM2);
 204           lvJpsiAOD = lvAodM1 + lvAodM2;
 205           invMAOD = lvJpsiAOD.M();
 206           fInvMassAOD->Fill(invMAOD);
 207         }// end if muons with opposite charge
 208       }// end if muon track
 209     }// end second loop over tracks in event
 210   }// end first loop over tracks in event
 211
 212   // Post final data. Write histo list to a file with option "RECREATE"
 213   PostData(0,fList);
 214 }
 215
 216 //________________________________________________________________________
 217 void AliAnalysisTaskAODvsESD::Terminate(const Option_t*)
 218 {
 219   // check if major differences between the two Ntuples (more comparisons can be added)
 220   int n1 = fMuonNtuple->GetEntries();
 221   int n2 = fMuonNtupleAOD->GetEntries();
 222
 223   if(n1!=n2){
 224     printf("ERROR: Different number of entries in single muon Ntuples\n");
 225     return;
 226   }
 227   else
 228     printf("Same number of entries in single muon Ntuples\n");
 229
 230 //   TCanvas* cv1 = new TCanvas("cvn1","cvn1",500,350);
 231 //   cv1->cd(1);
 232 //   fInvMass->SetMarkerStyle(29);
 233 //   fInvMass->SetMarkerColor(4);
 234 //   fInvMass->Draw();
 235 //   fInvMassAOD->SetMarkerStyle(28);
 236 //   fInvMassAOD->SetMarkerColor(2);
 237 //   fInvMassAOD->Draw("same");
 238 }
 239