MUON/fastMUONSim.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 /// \ingroup macros
  19 /// \file fastMUONSim.C
  20 /// \brief An example how to do the fast simulation and storing
  21 /// the muons that survive the reconstruction
  22 ///
  23 /// \author A. De Falco, H. Woehri, INFN Cagliari, April 2007
  24 ///
  25 /// An example how to do the fast simulation and storing
  26 /// the muons that survive the reconstruction in
  27 /// AliMUONTrackLight and AliMUONPairLight objects for
  28 /// further analysis. Should be used together with the
  29 /// macro fastMUONGen.C
  30
  31 #if !defined(__CINT__) || defined(__MAKECINT__)
  32 //STEER includes
  33 #include "AliStack.h"
  34 #include "AliRun.h"
  35 #include "AliHeader.h"
  36
  37 //FASTSIM includes
  38 #include "AliFastMuonTriggerEff.h"
  39 #include "AliFastMuonTrackingAcc.h"
  40 #include "AliFastMuonTrackingEff.h"
  41 #include "AliFastMuonTrackingRes.h"
  42 #include "AliFastDetector.h"
  43
  44 //MUON includes
  45 #include "AliMUONTrackLight.h"
  46 #include "AliMUONPairLight.h"
  47
  48 //ROOT includes
  49 #include "TParticle.h"
  50 #include "TRandom.h"
  51 #include "TFile.h"
  52 #include "TClonesArray.h"
  53 #include "TTree.h"
  54
  55 #endif
  56
  57 Float_t fThetaMin = 171., fThetaMax = 178.;
  58 //====================================================
  59 void fastMUONSim(Float_t bkg=0, char* outFileName = "fastSim_pp.root"){
  60
  61   //setting up the fast simulator
  62   AliFastMuonTriggerEff *trigeff = new AliFastMuonTriggerEff();
  63   AliFastMuonTrackingAcc *acc = new AliFastMuonTrackingAcc();
  64   AliFastMuonTrackingEff *eff = new AliFastMuonTrackingEff();
  65   AliFastMuonTrackingRes *res = new AliFastMuonTrackingRes();
  66   acc->SetBackground(bkg);
  67   eff->SetBackground(bkg);
  68   res->SetBackground(bkg);
  69   acc->Init();
  70   eff->Init();
  71   res->Init();
  72   trigeff->Init();
  73
  74   trigeff->SetCut(0); // 0... trigger pt cut at 1 GeV
  75                       // 1... trigger pt cut at 2 GeV
  76
  77   //prepare the arrays to store the generated/simulated muons
  78   TClonesArray *muonArray   = new TClonesArray("AliMUONTrackLight",100);
  79   TClonesArray *dimuonArray = new TClonesArray("AliMUONPairLight",100);
  80   TTree *treeOut = new TTree("tree","tree");
  81   TFile *fout = new TFile(outFileName,"recreate");
  82   fout->cd();
  83   treeOut->Branch("muons",&muonArray);
  84   treeOut->Branch("dimuons",&dimuonArray);
  85   //
  86   // Stack of particle for each event
  87   AliStack* stack;
  88   // Creating Run Loader and opening file containing Hits
  89   // AliRunLoader * RunLoader = AliRunLoader::Open("galice.root","MUONFolder","READ");
  90   AliRunLoader *RunLoader = AliRunLoader::Open("galice.root");
  91   if (RunLoader == 0x0) {
  92     printf(">>> Error : Error Opening %s file \n", "galice.root");
  93     return;
  94   }
  95
  96   RunLoader->LoadKinematics("READ");
  97   Int_t nevents = RunLoader->GetNumberOfEvents();
  98   printf("nevents %d\n", nevents);
  99
 100   TParticle *part = 0, *partMother=0;
 101   AliMUONTrackLight *muon = 0;
 102   Double_t radeg =  180./TMath::Pi();
 103   for(int iev = 0; iev < nevents; iev++){
 104
 105     if(iev%1000==0) printf("Event %d\n", iev);
 106     RunLoader->GetEvent(iev);
 107     stack = RunLoader->Stack();
 108     //stack->DumpPStack();
 109
 110     //reset muon info
 111     Int_t nMuons = 0;
 112     muonArray->Clear();     // clean muon and dimuon arrays
 113     dimuonArray->Clear();
 114     //
 115     Int_t nMuonsStored = 0;
 116     Int_t nPart = stack->GetNtrack();
 117     for(int iPart = 0; iPart < nPart; iPart++){
 118
 119       part = stack->Particle(iPart);
 120       if(fabs(part->GetPdgCode()) == 13){ //muon found
 121
 122         nMuons++;
 123         Double_t xyz[3] = {part->Vx(),part->Vy(),part->Vz()};
 124         // do not take muons that were decayed after the absorber
 125         if(xyz[2] < -(90.+ 40.)) continue; //decay muons after 90cm + 1 int. length
 126         Int_t charge = (part->GetPdgCode() > 0) ? -1 : +1;
 127         Double_t px = part->Px();
 128         Double_t py = part->Py();
 129         Double_t pz = part->Pz();
 130         Double_t E = part->Energy();
 131         TLorentzVector pGen(px,py,pz,E);
 132         //
 133         //fast simulation code
 134         res->SetCharge(charge);
 135         eff->SetCharge(charge);
 136         acc->SetCharge(charge);
 137         Float_t p = (Float_t) pGen.P();
 138         Float_t pt = (Float_t) pGen.Pt();
 139         Float_t theta = (Float_t) radeg*pGen.Theta();
 140         Float_t phi = (Float_t) radeg*pGen.Phi();
 141
 142         //continue only if we have a muon within the MUON acceptance
 143         if(theta < fThetaMin || theta > fThetaMax) continue;
 144
 145         theta = 180. - theta; // correct by hand the 'new' coordinates frame
 146         phi   = 180. - phi; // correct by hand the 'new' coordinates frame
 147         Float_t prec = 0, thetarec = 0, phirec = 0;
 148         res->Evaluate(p, theta ,phi, prec, thetarec, phirec);
 149         Float_t precx = prec * TMath::Sin(thetarec/radeg) * TMath::Cos(phirec/radeg);
 150         Float_t precy = prec * TMath::Sin(thetarec/radeg) * TMath::Sin(phirec/radeg);
 151         Float_t precz = prec * TMath::Cos(thetarec/radeg);
 152         precz = -precz; // correct by hand the 'new' coordinates frame
 153         precx   = -precx; // correct by hand the 'new' coordinates frame
 154
 155         Float_t trkeff = eff->Evaluate(charge, pt, theta, phi);
 156         Float_t accep  = acc->Evaluate(charge, pt, theta, phi);
 157         Float_t trgeff = trigeff->Evaluate(charge,pt,theta,phi);
 158         if(trkeff > 1.) printf("tracking efficiency > 1\n");
 159         if(accep > 1.) printf("acceptance efficiency > 1\n");
 160         if(trgeff > 1.) printf("trigger efficiency > 1\n");
 161         if (gRandom->Rndm() > trkeff || gRandom->Rndm() > accep) continue;
 162
 163         //only if we have a muon in the acceptance, store it as an
 164         //AliMUONTrackLight object and use it to form dimuons...
 165         muon = new AliMUONTrackLight();
 166         muon->SetCharge(charge);
 167         muon->SetTrackPDGCode(part->GetPdgCode()); //must be set, otherwise -999
 168         muon->SetTrackPythiaLine(iPart);//must be set, otherwise -999
 169         muon->SetPGen(pGen);
 170         muon->FillMuonHistory(stack, part);
 171         // set vertex to mother's vx,vy,vz to have the primary vertex
 172
 173         Int_t idMother = -1;
 174         Int_t id2 = 0;
 175         while (idMother < 0) idMother = muon->GetParentPythiaLine(id2++);
 176         partMother = stack->Particle(idMother);
 177         Double_t xyzMother[3] = {partMother->Vx(),
 178                                  partMother->Vy(),
 179                                  partMother->Vz()};
 180         muon->SetVertex(xyzMother);
 181         if (gRandom->Rndm() > trgeff) muon->SetTriggered(kFALSE);
 182         else muon->SetTriggered(kTRUE);
 183         muon->SetPxPyPz(precx,precy,precz);
 184
 185         //store the referenced track in the muonArray:
 186         TClonesArray &muons = *muonArray;
 187         new (muons[nMuonsStored++]) AliMUONTrackLight(*muon);
 188         delete muon;
 189       }
 190
 191     }//part
 192
 193     Int_t nmuons = muonArray->GetEntriesFast();
 194     Int_t ndimuons = 0;
 195     for(Int_t nMu1 = 0; nMu1 < nmuons-1; nMu1++) {
 196       AliMUONTrackLight* mu1 = (AliMUONTrackLight*) muonArray->At(nMu1);
 197       for(Int_t nMu2 = nMu1+1; nMu2 < nmuons; nMu2++){
 198         AliMUONTrackLight* mu2 = (AliMUONTrackLight*) muonArray->At(nMu2);
 199         AliMUONPairLight *dimuLight = new AliMUONPairLight();
 200         dimuLight->SetMuons(*mu1, *mu2);
 201         TClonesArray &dimuons = *dimuonArray;
 202         new (dimuons[ndimuons++]) AliMUONPairLight(*dimuLight);
 203       }
 204     }// end dimuons
 205     treeOut->Fill();
 206     stack->Reset();
 207   }//end of events
 208
 209   RunLoader->UnloadKinematics();
 210   fout->cd();
 211   treeOut->Write();
 212 }