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