MUON/DecodeRecoCocktail.C

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 // A. De Falco, H. Woehri, INFN Cagliari, July 2006
  19 // This macro reads the generation/reconstruction files in the
  20 // input directory (dirname),
  21 // builds a tree that contains, for each event, an array of muons and dimuons
  22 // (TClonesArrays of AliMUONTrackLight and AliMUONPairLight objects)
  23 // and writes the tree in an output file (outFileName)
  24 // Note that if the path for the output file is not explicitly specified,
  25 // it will be written in the directory containing the generation/reconstruction
  26
  27 #include <iostream>
  28 #include <TClonesArray.h>
  29 #include <TFile.h>
  30 #include <TTree.h>
  31 #include <TMath.h>
  32 #include <TLorentzVector.h>
  33 #include <TParticle.h>
  34 #include <TSystem.h>
  35 #include <AliMUONRecoCheck.h>
  36 #include <AliMUONTrack.h>
  37 #include <AliMUONTrackParam.h>
  38 #include <AliRunLoader.h>
  39 #include "AliMUONTrackLight.h"
  40 #include "AliMUONPairLight.h"
  41
  42 void DecodeRecoCocktail(char* dirname=".", char* outFileName = "MuonLight.root"){
  43   const char *startingDir = gSystem->pwd();
  44   gSystem->cd(dirname);
  45   TFile *fout = new TFile(outFileName,"recreate");
  46
  47   AliMUONRecoCheck *rc = new AliMUONRecoCheck("galice.root");
  48   AliRunLoader *runLoader = rc->GetRunLoader();
  49
  50   TClonesArray *muonArray   = new TClonesArray("AliMUONTrackLight",10);
  51   TClonesArray *dimuonArray = new TClonesArray("AliMUONPairLight",10);
  52   TTree *treeOut = new TTree("tree","tree");
  53
  54   fout->cd();
  55   treeOut->Branch("muons",&muonArray);
  56   treeOut->Branch("dimuons",&dimuonArray);
  57
  58   runLoader->LoadKinematics("READ");;
  59   Int_t nev = runLoader->GetNumberOfEvents();
  60
  61   TLorentzVector v;
  62   for(Int_t ievent = 0; ievent < nev; ievent++){ // loop over events
  63     printf ("Event %d of %d\n",ievent,nev);
  64     muonArray->Clear();     // clean muon and dimuon arrays
  65     dimuonArray->Clear();
  66     runLoader->GetEvent(ievent);
  67     rc->ResetTracks();
  68     rc->MakeTrackRef(); // make reconstructable tracks
  69
  70     TClonesArray * trackRecoArray = rc->GetTrackReco();
  71     TClonesArray * trackRefArray = rc->GetMuonTrackRef();
  72     Int_t nTrackReco = trackRecoArray->GetEntriesFast();
  73     Int_t nreftracks = 0;
  74     for (Int_t itrRec = 0; itrRec<nTrackReco; itrRec++) {
  75       // assign parameters concerning the reconstructed tracks
  76       AliMUONTrack *trackReco = (AliMUONTrack *)trackRecoArray->At(itrRec);
  77       AliMUONTrackLight *muLight = new AliMUONTrackLight();
  78       AliMUONTrackParam *trPar = trackReco->GetTrackParamAtVertex();
  79       muLight->SetCharge(Int_t(TMath::Sign(1.,trPar->GetInverseBendingMomentum())));
  80       muLight->SetPxPyPz(trPar->Px(),trPar->Py(), trPar->Pz());
  81       muLight->SetTriggered(trackReco->GetMatchTrigger());
  82       Double_t xyz[3] = { trPar->GetNonBendingCoor(),
  83                           trPar->GetBendingCoor(),
  84                           trPar->GetZ()};
  85       muLight->SetVertex(xyz);
  86       // find the reference track and store further information
  87       TParticle *part = muLight->FindRefTrack(trackReco,trackRefArray,runLoader);
  88       if (part) {
  89         v.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
  90         muLight->SetPGen(v);
  91         muLight->FillMuonHistory(runLoader, part);
  92         muLight->PrintInfo("A");
  93         //store the referenced track in the muonArray:
  94         TClonesArray &muons = *muonArray;
  95         new (muons[nreftracks++]) AliMUONTrackLight(*muLight);
  96       }
  97     }
  98
  99     // now loop over muon pairs to build dimuons
 100     Int_t nmuons = muonArray->GetEntriesFast();
 101     Int_t ndimuons = 0;
 102     for(Int_t itrRec1 = 0; itrRec1 < nmuons-1; itrRec1++) {
 103       AliMUONTrackLight* mu1 = (AliMUONTrackLight*) muonArray->At(itrRec1);
 104       for(Int_t itrRec2 = itrRec1+1; itrRec2 < nmuons; itrRec2++){
 105         AliMUONTrackLight* mu2 = (AliMUONTrackLight*) muonArray->At(itrRec2);
 106         AliMUONPairLight *dimuLight = new AliMUONPairLight();
 107         dimuLight->SetMuons(*mu1, *mu2);
 108         //      if(dimuLight->GetCreationProcess() == 2){
 109         //        printf ("#dimuon = %d (%d, %d) \n", ndimuons, itrRec1, itrRec2);
 110         // dimuLight->PrintInfo("A");
 111         //      }
 112         //store the referenced track in the dimuonArray:
 113         TClonesArray &dimuons = *dimuonArray;
 114         new (dimuons[ndimuons++]) AliMUONPairLight(*dimuLight);
 115       }
 116     }
 117     Int_t ndimu2 = dimuonArray->GetEntriesFast();
 118     //     printf ("dimuonArray has %d entries\n",ndimu2);
 119     //    dimuonArray->Dump();
 120     //     printf ("filling tree\n");
 121     //     // fill the tree
 122     treeOut->Fill();
 123     //     printf ("done\n");
 124
 125   }
 126   fout->cd();
 127   treeOut->Write();
 128   gSystem->cd(startingDir);
 129 }