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-Nov-2006: modified by in order to loop on files
  27
  28 #if !defined(__CINT__) || defined(__MAKECINT__)
  29 #include <Riostream.h>
  30 #include <TClonesArray.h>
  31 #include <TFile.h>
  32 #include <TTree.h>
  33 #include <TMath.h>
  34 #include <TLorentzVector.h>
  35 #include <TParticle.h>
  36 #include <TSystem.h>
  37 #include <TGeoManager.h>
  38 #include "AliMUONRecoCheck.h"
  39 #include "AliMUONTrack.h"
  40 #include "AliMUONTrackParam.h"
  41 #include "AliRunLoader.h"
  42 #include "AliMUONTrackLight.h"
  43 #include "AliMUONPairLight.h"
  44 #include "AliMUONTrackExtrap.h"
  45 #include "AliESD.h"
  46 #include "AliESDVertex.h"
  47 #include "AliITSVertexerPPZ.h"
  48 #include "AliITSLoader.h"
  49 #include "AliTracker.h"
  50 #include "AliMagFMaps.h"
  51 #endif
  52
  53
  54 void DecodeRecoCocktail(char* dirname=".", char* outFileName = "MuonLight.root", char* geoFilename = "geometry.root"){
  55
  56
  57   char startingDir[200];
  58   sprintf (startingDir,"%s",gSystem->pwd());
  59   gSystem->cd(dirname);
  60
  61
  62   TClonesArray *muonArray   = new TClonesArray("AliMUONTrackLight",10);
  63   TClonesArray *dimuonArray = new TClonesArray("AliMUONPairLight",10);
  64   TTree *treeOut = new TTree("tree","tree");
  65
  66   TFile *fout = new TFile(outFileName,"recreate");
  67   fout->cd();
  68
  69   treeOut->Branch("muons",&muonArray);
  70   treeOut->Branch("dimuons",&dimuonArray);
  71
  72   TFile* esdFile = TFile::Open("AliESDs.root");
  73   if (!esdFile || !esdFile->IsOpen()) {
  74     Error("DecodeRecoCocktailNew", "opening ESD file %s failed", esdFileName);
  75     return;
  76   }
  77
  78   AliESD* esd = new AliESD();
  79   TTree* treeESD = (TTree*) esdFile->Get("esdTree");
  80   if (!treeESD) {
  81     Error("CheckESD", "no ESD tree found");
  82     return;
  83   }
  84   treeESD->SetBranchAddress("ESD", &esd);
  85
  86   // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
  87   if (!gGeoManager) {
  88     TGeoManager::Import(geoFilename);
  89     if (!gGeoManager) {
  90       Error("MUONmass_ESD", "getting geometry from file %s failed", geoFilename);
  91       return;
  92     }
  93   }
  94
  95   // set  mag field
  96   // waiting for mag field in CDB
  97   printf("Loading field map...\n");
  98   AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 1, 1., 10., AliMagFMaps::k5kG);
  99   AliTracker::SetFieldMap(field, kFALSE);
 100   // set the magnetic field for track extrapolations
 101   AliMUONTrackExtrap::SetField(AliTracker::GetFieldMap());
 102
 103   AliMUONRecoCheck *rc = new AliMUONRecoCheck("galice.root");
 104   AliRunLoader *runLoader = rc->GetRunLoader();
 105
 106   runLoader->LoadKinematics("READ");
 107   Int_t nev = runLoader->GetNumberOfEvents();
 108 //   nevent = nevent +nev;
 109 //     printf(" number of files and events = %d - %d \n",irun,nevent);
 110
 111   AliITSLoader* ITSloader =  (AliITSLoader*) runLoader->GetLoader("ITSLoader");
 112   AliITSVertexerPPZ *dovert = 0;
 113   if (ITSloader) {
 114     dovert = new AliITSVertexerPPZ("default",0,0);
 115     dovert->SetDebug(0);
 116     dovert->SetDiffPhiMax(0.05);
 117     dovert->SetWindow(3.);
 118   }
 119   AliESDVertex *vert = 0;
 120
 121   TLorentzVector v;
 122
 123   for(Int_t ievent = 0; ievent < nev; ievent++){ // loop over events
 124     AliHeader *header = runLoader->GetHeader();
 125     if (ITSloader) {
 126       vert = dovert->FindVertexForCurrentEvent(ievent);
 127     }
 128     // printf ("Event %d of %d\n",ievent,nev);
 129     muonArray->Clear();     // clean muon and dimuon arrays
 130     dimuonArray->Clear();
 131     runLoader->GetEvent(ievent);
 132     treeESD->GetEvent(ievent);
 133     rc->ResetTracks();
 134     rc->MakeTrackRef(); // make reconstructable tracks
 135
 136     TClonesArray * trackRecoArray = rc->GetTrackReco();
 137     TClonesArray * trackRefArray = rc->GetMuonTrackRef();
 138     Int_t nTrackReco = trackRecoArray->GetEntriesFast();
 139     Int_t nTracksESD = (Int_t)esd->GetNumberOfMuonTracks() ;
 140     if (nTrackReco != nTracksESD) printf ("Tracks in AliMUONTrack and in ESD do not match!\n");
 141     Int_t nreftracks = 0;
 142     for (Int_t itrRec = 0; itrRec<nTrackReco; itrRec++) {
 143       // assign parameters concerning the reconstructed tracks
 144       AliMUONTrackLight muLight;
 145       AliMUONTrack *trackReco = (AliMUONTrack *)trackRecoArray->At(itrRec);
 146       muLight.FillFromESD(esd->GetMuonTrack(itrRec));
 147       //        muLight.FillFromAliMUONTrack(trackReco);
 148
 149       // find the reference track and store further information
 150       TParticle *part = muLight.FindRefTrack(trackReco,trackRefArray,runLoader);
 151       if (part) {
 152         v.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
 153         muLight.SetPGen(v);
 154         muLight.FillMuonHistory(runLoader, part);
 155         //        muLight.PrintInfo("A");
 156         //store the referenced track in the muonArray:
 157         TClonesArray &muons = *muonArray;
 158         new (muons[nreftracks++]) AliMUONTrackLight(muLight);
 159       }
 160     }  // end reco track
 161
 162     // now loop over muon pairs to build dimuons
 163     Int_t nmuons = muonArray->GetEntriesFast();
 164     Int_t ndimuons = 0;
 165     for(Int_t itrRec1 = 0; itrRec1 < nmuons-1; itrRec1++) {
 166       AliMUONTrackLight* mu1 = (AliMUONTrackLight*) muonArray->At(itrRec1);
 167       for(Int_t itrRec2 = itrRec1+1; itrRec2 < nmuons; itrRec2++){
 168         AliMUONTrackLight* mu2 = (AliMUONTrackLight*) muonArray->At(itrRec2);
 169         AliMUONPairLight dimuLight;
 170         dimuLight.SetMuons(*mu1, *mu2);
 171         TClonesArray &dimuons = *dimuonArray;
 172         new (dimuons[ndimuons++]) AliMUONPairLight(dimuLight);
 173       }
 174     }
 175     treeOut->Fill();
 176   }
 177
 178   fout->cd();
 179   treeOut->Write();
 180   gSystem->cd(startingDir);
 181   fout->Close();
 182   delete muonArray;
 183   delete dimuonArray;
 184   delete treeOut;
 185   delete fout;
 186   delete rc;
 187   delete esd;
 188 }
 189