PWG3/muondep/AliAnalysisTaskRecoCheck.cxx

   1 #include "TChain.h"
   2 #include "TClonesArray.h"
   3 #include "TLorentzVector.h"
   4 #include "TIterator.h"
   5
   6 #include "AliStack.h"
   7 #include "AliMagF.h"
   8 #include "AliTracker.h"
   9
  10 #include "AliAnalysisTask.h"
  11 #include "AliAnalysisManager.h"
  12
  13 #include "AliESDEvent.h"
  14 #include "AliESDInputHandler.h"
  15 #include "AliESDMuonTrack.h"
  16
  17 #include "AliMCEventHandler.h"
  18 #include "AliMCEvent.h"
  19
  20 #include "AliAnalysisTaskRecoCheck.h"
  21
  22 #include "AliMUONTrackLight.h"
  23 #include "AliMUONPairLight.h"
  24 #include "AliMUONVTrackStore.h"
  25 #include "AliMUONTrack.h"
  26 #include "AliMUONRecoCheck.h"
  27 #include "AliMUONTrackExtrap.h"
  28
  29 // analysis task for decoding reconstructed tracks and kinematics (AliMUONRecoCheck)
  30 // Authors: Bogdan Vulpescu
  31
  32 ClassImp(AliAnalysisTaskRecoCheck)
  33
  34 //________________________________________________________________________
  35 AliAnalysisTaskRecoCheck::AliAnalysisTaskRecoCheck(const char *name)
  36   : AliAnalysisTask(name, ""), fESDEvent(0), fTree(0), fArray1Mu(0), fArray2Mu(0), fL3Current(30000.0)
  37 {
  38   // Constructor
  39
  40   // Define input and output slots here
  41   // Input slot #0 works with a TChain
  42   DefineInput(0, TChain::Class());
  43   // Output slot #0 writes into a TTree container
  44   DefineOutput(0, TTree::Class());
  45
  46 }
  47
  48 //________________________________________________________________________
  49 void AliAnalysisTaskRecoCheck::ConnectInputData(Option_t *)
  50 {
  51   // Connect ESD or AOD here
  52   // Called once
  53
  54   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  55   if (!tree) {
  56     Printf("ERROR: Could not read chain from input slot 0");
  57   } else {
  58     tree->SetBranchStatus("*", kFALSE);
  59     tree->SetBranchStatus("MuonTracks.*", kTRUE);
  60     tree->SetBranchStatus("AliESDHeader.*", kTRUE);
  61
  62     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  63
  64     if (!esdH) {
  65       Printf("ERROR: Could not get ESDInputHandler");
  66     } else
  67       fESDEvent = esdH->GetEvent();
  68   }
  69
  70   // calculate the filed map in the L3 magnet using the current value
  71   AliMUONTrackExtrap::SetField();
  72
  73 }
  74
  75 //________________________________________________________________________
  76 void AliAnalysisTaskRecoCheck::CreateOutputObjects()
  77 {
  78   // Create output tree with branches of arrays
  79   // Called once
  80
  81   fArray1Mu = new TClonesArray("AliMUONTrackLight",0);
  82   fArray2Mu = new TClonesArray("AliMUONPairLight",0);
  83   fTree = new TTree("tree","tree",0);
  84   fTree->Branch("muons",  &fArray1Mu);
  85   fTree->Branch("dimuons",&fArray2Mu);
  86   fTree->SetDirectory(0);
  87
  88 }
  89
  90 //________________________________________________________________________
  91 void AliAnalysisTaskRecoCheck::Exec(Option_t *)
  92 {
  93   // Main loop
  94   // Called for each event
  95
  96   // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
  97   // This handler can return the current MC event
  98
  99   Int_t nTracks = fESDEvent->GetNumberOfMuonTracks();
 100   if (nTracks == 0) return;
 101
 102   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 103   if (!eventHandler) {
 104     Printf("ERROR: Could not retrieve MC event handler");
 105     return;
 106   }
 107
 108   AliMCEvent* mcEvent = eventHandler->MCEvent();
 109   if (!mcEvent) {
 110      Printf("ERROR: Could not retrieve MC event");
 111      return;
 112   }
 113
 114   Int_t eventNumber = fESDEvent->GetEventNumberInFile();
 115
 116   fArray1Mu->Clear();
 117   fArray2Mu->Clear();
 118
 119   TLorentzVector v;
 120
 121   AliMUONRecoCheck *rc = new AliMUONRecoCheck(fESDEvent, eventHandler);
 122
 123   AliMUONVTrackStore* trackRefs = rc->TrackRefs(eventNumber);
 124   AliMUONVTrackStore* recoTracks = rc->ReconstructedTracks(eventNumber);
 125   AliStack* pstack = mcEvent->Stack();
 126
 127   TIter next(recoTracks->CreateIterator());
 128   AliMUONTrack* trackReco = NULL;
 129
 130   Int_t nTrackReco = recoTracks->GetSize();
 131   Int_t nTracksESD = fESDEvent->GetNumberOfMuonTracks();
 132   //printf("Tracks in recoTracks (%d) and in ESD (%d).\n", nTrackReco, nTracksESD);
 133
 134   if (nTrackReco != nTracksESD) printf ("Tracks in recoTracks (%d) and in ESD (%d) do not match!\n", nTrackReco, nTracksESD);
 135
 136   Int_t nreftracks = 0;
 137   Int_t itrRec = 0;
 138   while ( (trackReco = static_cast<AliMUONTrack*>(next())) != NULL ) {
 139     // assign parameters concerning the reconstructed tracks
 140     AliMUONTrackLight muLight;
 141
 142     muLight.FillFromESD(fESDEvent->GetMuonTrack(itrRec));
 143
 144     // find the reference track and store further information
 145     TParticle *part = muLight.FindRefTrack(trackReco, trackRefs, pstack);
 146     if (part) {
 147       v.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
 148       muLight.SetPGen(v);
 149       muLight.FillMuonHistory(pstack, part);
 150       //muLight.PrintInfo("A");
 151       //store the referenced track in the muonArray:
 152       TClonesArray &muons = *fArray1Mu;
 153       new (muons[nreftracks++]) AliMUONTrackLight(muLight);
 154     }
 155
 156     itrRec++;
 157   }  // end reco track
 158
 159   // now loop over muon pairs to build dimuons
 160   Int_t nmuons = fArray1Mu->GetEntriesFast();
 161   Int_t ndimuons = 0;
 162   for(Int_t itrRec1 = 0; itrRec1 < nmuons-1; itrRec1++) {
 163     AliMUONTrackLight* mu1 = (AliMUONTrackLight*) fArray1Mu->At(itrRec1);
 164     for(Int_t itrRec2 = itrRec1+1; itrRec2 < nmuons; itrRec2++){
 165       AliMUONTrackLight* mu2 = (AliMUONTrackLight*) fArray1Mu->At(itrRec2);
 166       AliMUONPairLight dimuLight;
 167       dimuLight.SetMuons(*mu1, *mu2);
 168       //dimuLight.PrintInfo("A");
 169       TClonesArray &dimuons = *fArray2Mu;
 170       new (dimuons[ndimuons++]) AliMUONPairLight(dimuLight);
 171     }
 172   }
 173
 174   delete rc;
 175
 176   // Post output data.
 177   if (nreftracks != 0) {
 178     fTree->Fill();
 179   }
 180
 181   PostData(0, fTree);
 182
 183
 184 }
 185
 186 //________________________________________________________________________
 187 void AliAnalysisTaskRecoCheck::Terminate(Option_t *)
 188 {
 189   // the end
 190
 191 }