Generator for cosmic muons. (A. Dainese)
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 4 Jun 2007 14:03:30 +0000 (14:03 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 4 Jun 2007 14:03:30 +0000 (14:03 +0000)
EVGEN/AliGenCosmicsParam.cxx [new file with mode: 0644]
EVGEN/AliGenCosmicsParam.h [new file with mode: 0644]
EVGEN/EVGENLinkDef.h
EVGEN/libEVGEN.pkg

diff --git a/EVGEN/AliGenCosmicsParam.cxx b/EVGEN/AliGenCosmicsParam.cxx
new file mode 100644 (file)
index 0000000..a8f506d
--- /dev/null
@@ -0,0 +1,204 @@
+/**************************************************************************
+ * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+
+// Generator for muons according to kinematic parametrizations at ALICE
+// (not at the surface).
+// Origin: andrea.dainese@lnl.infn.it
+
+#include <TParticle.h>
+#include <TF1.h>
+
+#include "AliRun.h"
+#include "AliESDtrack.h"
+#include "AliESDVertex.h"
+#include "AliGenCosmicsParam.h"
+
+ClassImp(AliGenCosmicsParam)
+
+//-----------------------------------------------------------------------------
+AliGenCosmicsParam::AliGenCosmicsParam():
+AliGenerator(),
+fParamMI(kFALSE),
+fParamACORDE(kFALSE),
+fYOrigin(600.),
+fMaxAngleWRTVertical(-99.),
+fBkG(0.),
+fTPC(kFALSE),
+fITS(kFALSE),
+fSPDouter(kFALSE),
+fSPDinner(kFALSE)
+{
+  //
+  // Default constructor
+  //
+}
+//-----------------------------------------------------------------------------
+void AliGenCosmicsParam::Generate()
+{
+  //
+  // Generate one muon
+  //
+  
+  //
+  Float_t origin[3];
+  Float_t p[3];
+  Int_t nt;
+  Double_t ptot=0,pt=0,angleWRTVertical=0;
+  Bool_t okMom=kFALSE,okAngle=kFALSE;
+  //
+
+  // mu+ or mu-
+  Int_t ipart=13;
+  if(gRandom->Rndm()<0.5) ipart *= -1;
+
+  if(fParamACORDE) { // extracted from AliGenACORDE events
+    // sample total momentum only once (to speed up)
+    TF1 *dNdpACORDE = new TF1("dNdpACORDE","x/(1.+(x/12.8)*(x/12.8))^1.96",fPMin,fPMax);
+    ptot = (Double_t)dNdpACORDE->GetRandom();
+    delete dNdpACORDE;
+    dNdpACORDE = 0;
+  }
+
+  Int_t trials=0;
+
+  while(1) {
+    trials++;
+    // origin
+    origin[0]  = fYOrigin*TMath::Tan(fMaxAngleWRTVertical)*(-1.+2.*gRandom->Rndm());
+    origin[1]  = fYOrigin;
+    origin[2]  = fYOrigin*TMath::Tan(fMaxAngleWRTVertical)*(-1.+2.*gRandom->Rndm());
+
+    // momentum
+    while(1) {
+      okMom=kFALSE; okAngle=kFALSE;
+
+      if(fParamMI) { // parametrization by M.Ivanov of LEP cosmics data
+       Float_t pref  = 1. + gRandom->Exp(30.);
+       p[1] = -pref; 
+       p[0] = gRandom->Gaus(0.0,0.2)*pref;
+       p[2] = gRandom->Gaus(0.0,0.2)*pref;
+       if(gRandom->Rndm()>0.9) {
+         p[0] = gRandom->Gaus(0.0,0.4)*pref;
+         p[2] = gRandom->Gaus(0.0,0.4)*pref;
+       }
+       ptot=TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
+       pt=TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
+      } else if(fParamACORDE) { // extracted from AliGenACORDE events
+       Float_t theta,phi;
+       while(1) {
+         theta = gRandom->Gaus(0.5*TMath::Pi(),0.42);
+         if(TMath::Abs(theta-0.5*TMath::Pi())<fMaxAngleWRTVertical) break;
+       }
+       while(1) {
+         phi = gRandom->Gaus(-0.5*TMath::Pi(),0.42);
+         if(TMath::Abs(phi+0.5*TMath::Pi())<fMaxAngleWRTVertical) break;
+       }
+       pt = ptot*TMath::Sin(theta);
+       p[0] = pt*TMath::Cos(phi); 
+       p[1] = pt*TMath::Sin(phi); 
+       p[2] = ptot*TMath::Cos(theta);
+      } else {
+       AliFatal("Parametrization not set: use SetParamMI or SetParamACORDE");
+      }
+
+      
+      // check kinematic cuts
+      if(TestBit(kMomentumRange)) {
+       if(ptot>fPMin && ptot<fPMax) okMom=kTRUE;
+      }
+
+      angleWRTVertical=TMath::ACos(TMath::Abs(p[1])/ptot); // acos(|py|/ptot)
+      if(angleWRTVertical<fMaxAngleWRTVertical) okAngle=kTRUE;
+      
+      if(okAngle&&okMom) break;
+    }
+
+    // acceptance
+    if(!fTPC && !fITS && !fSPDinner && !fSPDouter) break;
+    if(fTPC) if(IntersectCylinder(250.,250.,ipart,origin,p)) break;
+    if(fITS) if(IntersectCylinder(50.,50.,ipart,origin,p)) break;
+    if(fSPDouter) if(IntersectCylinder(6.5,14.,ipart,origin,p)) break;
+    if(fSPDinner) if(IntersectCylinder(3.5,14.0,ipart,origin,p)) break;
+  }
+
+  Float_t polarization[3]= {0,0,0};
+  PushTrack(fTrackIt,-1,ipart,p,origin,polarization,0,kPPrimary,nt);
+
+  //printf("TRIALS %d\n",trials);
+  
+
+  return;
+}
+//-----------------------------------------------------------------------------
+void AliGenCosmicsParam::Init()
+{
+  // 
+  // Initialisation, check consistency of selected ranges
+  //
+  if(TestBit(kPtRange)) 
+    AliFatal("You cannot set the pt range for this generator! Only momentum range");
+  if(fPMin<8.) { 
+    fPMin=8.; 
+    if(TestBit(kMomentumRange)) 
+      AliWarning("Minimum momentum cannot be < 8 GeV/c"); 
+  }
+  if(fMaxAngleWRTVertical<0.) 
+    AliFatal("You must use SetMaxAngleWRTVertical() instead of SetThetaRange(), SetPhiRange()");
+
+  printf("************ AliGenCosmicsParam ****************\n");
+  printf("***** Muons generated at Y = %f cm\n",fYOrigin);
+  printf("************************************************\n");
+
+  return;
+}
+//-----------------------------------------------------------------------------
+Bool_t AliGenCosmicsParam::IntersectCylinder(Float_t r,Float_t z,Int_t pdg,
+                                            Float_t o[3],Float_t p[3]) const
+{
+  //
+  // Intersection between muon and cylinder [-z,+z] with radius r
+  //
+
+  Float_t en = TMath::Sqrt(0.105*0.105+p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
+  TParticle part(pdg,0,0,0,0,0,p[0],p[1],p[2],en,o[0],o[1],o[2],0);
+  AliESDtrack track(&part);
+  Double_t pos[3]={0.,0.,0.},sigma[3]={0.,0.,0.};
+  AliESDVertex origin(pos,sigma);
+
+  track.RelateToVertex(&origin,fBkG,10000.);
+
+  Float_t d0z0[2],covd0z0[3];
+  track.GetImpactParameters(d0z0,covd0z0);
+
+  // check rphi 
+  if(TMath::Abs(d0z0[0])>r) return kFALSE;
+  // check z
+  if(TMath::Abs(d0z0[1])>z) return kFALSE;
+
+  /*
+    if(TMath::Abs(fB)<0.01) {  // NO FIELD
+    Float_t drphi = TMath::Abs(o[1]-p[1]/p[0]*o[0])/
+    TMath::Sqrt(p[1]*p[1]/p[0]/p[0]+1.);
+    if(drphi>r) return kFALSE;
+    Float_t dz = o[2]-p[2]/p[0]*o[0]+p[2]/p[0]*
+    (p[1]*p[1]/p[0]/p[0]*o[0]-p[1]/p[0]*o[1])/(1.+p[1]*p[1]/p[0]/p[0]);
+    if(TMath::Abs(dz)>z) return kFALSE;
+    }
+  */
+
+  return kTRUE;
+}
+//-----------------------------------------------------------------------------
diff --git a/EVGEN/AliGenCosmicsParam.h b/EVGEN/AliGenCosmicsParam.h
new file mode 100644 (file)
index 0000000..8582e26
--- /dev/null
@@ -0,0 +1,53 @@
+#ifndef ALIGENCOSMICSPARAM_H
+#define ALIGENCOSMICSPARAM_H
+/* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// Generator for muons according to kinematic parametrizations at ALICE
+// (not at the surface).
+// Origin: andrea.dainese@lnl.infn.it
+
+
+#include "AliLog.h"
+#include "AliGenerator.h"
+
+class AliGenCosmicsParam : public AliGenerator
+{
+public:
+
+  AliGenCosmicsParam();
+  virtual ~AliGenCosmicsParam() {}
+  virtual void Generate();
+  virtual void Init();
+  void SetParamMI() { fParamMI=kTRUE; fParamACORDE=kFALSE; return; }
+  void SetParamACORDE() { fParamMI=kFALSE; fParamACORDE=kTRUE; return; }
+  void SetYOrigin(Float_t y=600.) { fYOrigin=y; return; }
+  void SetMaxAngleWRTVertical(Float_t max=45.) { 
+      if(max<0. || max>90.) AliFatal("angle must be in [0,pi/2]");
+      fMaxAngleWRTVertical=max; return; }
+  void SetBkG(Float_t b) { fBkG=b; return; }
+  void SetInTPC() { fTPC=kTRUE; return; }
+  void SetInITS() { fITS=kTRUE; return; }
+  void SetInSPDouter() { fSPDouter=kTRUE; return; }
+  void SetInSPDinner() { fSPDinner=kTRUE; return; }
+
+private:
+
+  Bool_t IntersectCylinder(Float_t r,Float_t z,Int_t pdg,
+                          Float_t o[3],Float_t p[3]) const;  
+
+  Bool_t fParamMI;              // parametrization from M.Ivanov
+  Bool_t fParamACORDE;          // parametrization from AliGenACORDE 
+                                // (parametrized at ALICE y)
+  Float_t fYOrigin;             // y of muon origin
+  Float_t fMaxAngleWRTVertical; // maximum angle between momentum and y axis
+  Float_t fBkG;                 // field in kGauss
+  Bool_t fTPC;                  // acceptance cuts
+  Bool_t fITS;                  // acceptance cuts
+  Bool_t fSPDouter;             // acceptance cuts
+  Bool_t fSPDinner;             // acceptance cuts
+
+  ClassDef(AliGenCosmicsParam,1) // parametrized cosmics generator
+};
+
+#endif
index 7638660..a228998 100644 (file)
@@ -56,4 +56,5 @@
 #pragma link C++ class  AliGenReaderEMD+;
 #pragma link C++ class  AliDecayerPolarized+;
 #pragma link C++ class  AliGenCorrHF+;
+#pragma link C++ class  AliGenCosmicsParam+;
 #endif
index cb1bbf1..4731810 100644 (file)
@@ -17,7 +17,7 @@ SRCS          = AliGenHIJINGpara.cxx AliGenBox.cxx AliGenFixed.cxx \
                AliGenGeVSimEventHeader.cxx\
                AliSlowNucleonModel.cxx AliSlowNucleonModelExp.cxx \
                AliGenMUONCocktail.cxx AliGenMUONCocktailpp.cxx AliGenHBTosl.cxx AliGenCocktailEventHeader.cxx \
-               AliGenReaderEMD.cxx AliDecayerPolarized.cxx AliGenCorrHF.cxx
+               AliGenReaderEMD.cxx AliDecayerPolarized.cxx AliGenCorrHF.cxx AliGenCosmicsParam.cxx
 
 # Headerfiles for this particular package (Path respect to own directory)
 HDRS= $(SRCS:.cxx=.h)