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 AliESDs.root failed");
  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", "galice_sim.root");
 104   AliRunLoader *runLoader = rc->GetRunLoader();
 105   AliRunLoader *runLoaderSim = rc->GetRunLoaderSim();
 106
 107   runLoaderSim->LoadKinematics("READ");
 108   Int_t nev = runLoaderSim->GetNumberOfEvents();
 109 //   nevent = nevent +nev;
 110 //     printf(" number of files and events = %d - %d \n",irun,nevent);
 111
 112   AliITSLoader* ITSloader =  (AliITSLoader*) runLoaderSim->GetLoader("ITSLoader");
 113   AliITSVertexerPPZ *dovert = 0;
 114   if (ITSloader) {
 115     dovert = new AliITSVertexerPPZ("default",0,0);
 116     // dovert->SetDebug(0);
 117     dovert->SetDiffPhiMax(0.05);
 118     dovert->SetWindow(3.);
 119   }
 120   AliESDVertex *vert = 0;
 121
 122   TLorentzVector v;
 123
 124   for(Int_t ievent = 0; ievent < nev; ievent++){ // loop over events
 125     runLoaderSim->GetHeader();
 126     if (ITSloader) {
 127       vert = dovert->FindVertexForCurrentEvent(ievent);
 128     }
 129     //printf ("Event %d of %d\n",ievent,nev);
 130     muonArray->Clear();     // clean muon and dimuon arrays
 131     dimuonArray->Clear();
 132     runLoader->GetEvent(ievent);
 133     runLoaderSim->GetEvent(ievent);
 134     treeESD->GetEvent(ievent);
 135     rc->ResetTracks();
 136     rc->MakeTrackRef(); // make reconstructable tracks
 137
 138     TClonesArray * trackRecoArray = rc->GetTrackReco();
 139     TClonesArray * trackRefArray = rc->GetMuonTrackRef();
 140     Int_t nTrackReco = trackRecoArray->GetEntriesFast();
 141     Int_t nTracksESD = (Int_t)esd->GetNumberOfMuonTracks() ;
 142     if (nTrackReco != nTracksESD) printf ("Tracks in AliMUONTrack and in ESD do not match!\n");
 143     Int_t nreftracks = 0;
 144     for (Int_t itrRec = 0; itrRec<nTrackReco; itrRec++) {
 145       // assign parameters concerning the reconstructed tracks
 146       AliMUONTrackLight muLight;
 147       AliMUONTrack *trackReco = (AliMUONTrack *)trackRecoArray->At(itrRec);
 148
 149       muLight.FillFromESD(esd->GetMuonTrack(itrRec));
 150       //        muLight.FillFromAliMUONTrack(trackReco);
 151
 152       // find the reference track and store further information
 153       TParticle *part = muLight.FindRefTrack(trackReco,trackRefArray,runLoaderSim);
 154       if (part) {
 155         v.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->Energy());
 156         muLight.SetPGen(v);
 157         muLight.FillMuonHistory(runLoaderSim, part);
 158         //        muLight.PrintInfo("A");
 159         //store the referenced track in the muonArray:
 160         TClonesArray &muons = *muonArray;
 161         new (muons[nreftracks++]) AliMUONTrackLight(muLight);
 162       }
 163     }  // end reco track
 164
 165     // now loop over muon pairs to build dimuons
 166     Int_t nmuons = muonArray->GetEntriesFast();
 167     Int_t ndimuons = 0;
 168     for(Int_t itrRec1 = 0; itrRec1 < nmuons-1; itrRec1++) {
 169       AliMUONTrackLight* mu1 = (AliMUONTrackLight*) muonArray->At(itrRec1);
 170       for(Int_t itrRec2 = itrRec1+1; itrRec2 < nmuons; itrRec2++){
 171         AliMUONTrackLight* mu2 = (AliMUONTrackLight*) muonArray->At(itrRec2);
 172         AliMUONPairLight dimuLight;
 173         dimuLight.SetMuons(*mu1, *mu2);
 174         TClonesArray &dimuons = *dimuonArray;
 175         new (dimuons[ndimuons++]) AliMUONPairLight(dimuLight);
 176       }
 177     }
 178     treeOut->Fill();
 179   }
 180
 181   fout->cd();
 182   treeOut->Write();
 183   gSystem->cd(startingDir);
 184   fout->Close();
 185   delete muonArray;
 186   delete dimuonArray;
 187   delete treeOut;
 188   delete fout;
 189   delete rc;
 190   delete esd;
 191 }
 192