]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HMPID/AliHMPIDv3.cxx
minors
[u/mrichter/AliRoot.git] / HMPID / AliHMPIDv3.cxx
index b37bf7472997a1ddac5dd00844de4b230e9eb45e..fd84ee67581b666c947d1f52599f87eab6bc4d77 100644 (file)
 #include <AliMC.h>              //StepManager()      
 #include <AliRun.h>             //CreateMaterials()    
 #include <AliMagF.h>            //CreateMaterials()
+#include "AliGeomManager.h"   //AddAlignableVolumes()
+#include <AliCDBEntry.h>        //CreateMaterials()
+#include <AliCDBManager.h>      //CreateMaterials()
 #include <TF1.h>                //DefineOpticalProperties()
 #include <TF2.h>                //DefineOpticalProperties()
-#include <TLorentzVector.h>     //IsLostByFresnel() 
-#include <AliCDBManager.h>      //CreateMaterials()
-#include <AliCDBEntry.h>        //CreateMaterials()
+#include <TGeoCompositeShape.h> //CradleBaseVolume()
+#include <TGeoGlobalMagField.h>
 #include <TGeoPhysicalNode.h>   //AddAlignableVolumes()
 #include <TGeoXtru.h>           //CradleBaseVolume()
-#include <TGeoCompositeShape.h> //CradleBaseVolume()
+#include <TLorentzVector.h>     //IsLostByFresnel() 
 #include <TString.h>            //StepManager()
-#include "AliGeomManager.h"   //AddAlignableVolumes()
+#include <TTree.h>
 
 ClassImp(AliHMPIDv3)    
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -60,7 +62,7 @@ void AliHMPIDv3::AddAlignableVolumes()const
   for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++) {
     modUID = AliGeomManager::LayerToVolUID(idHMPID,modnum++);
     if(!gGeoManager->SetAlignableEntry(Form("/HMPID/Chamber%i",iCh),Form("ALIC_1/Hmp%i_0",iCh),modUID))
-           AliFatal("AliHMPIDv3::Unable to set alignable entry!!");  //aligment without AliCluster3D
+           AliError("AliHMPIDv3::Unable to set alignable entry!!");  //aligment without AliCluster3D
     //Get Tracking To Local matricies for alignment with AliCluster3D
     TGeoPNEntry *eCh = gGeoManager->GetAlignableEntryByUID(modUID);
     TGeoHMatrix *globMatrix = eCh->GetGlobalOrig();
@@ -99,8 +101,8 @@ void AliHMPIDv3::CreateMaterials()
 
     Int_t   matId=0;                           //tmp material id number
     Int_t   unsens =  0, sens=1;               //sensitive or unsensitive medium
-    Int_t   itgfld = gAlice->Field()->Integ(); //type of field intergration 0 no field -1 user in guswim 1 Runge Kutta 2 helix 3 const field along z
-    Float_t maxfld = gAlice->Field()->Max();   //max field value
+    Int_t   itgfld = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Integ(); //type of field intergration 0 no field -1 user in guswim 1 Runge Kutta 2 helix 3 const field along z
+    Float_t maxfld = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Max();   //max field value
     Float_t tmaxfd = -10.0;                    //max deflection angle due to magnetic field in one step
     Float_t deemax = - 0.2;                    //max fractional energy loss in one step   
     Float_t stemax = - 0.1;                    //max step allowed [cm]
@@ -131,13 +133,11 @@ void AliHMPIDv3::CreateMaterials()
     AliMaterial(++matId,"W"   ,aW   ,zW   ,dW   ,radW   ,absW   );  AliMedium(kW   ,"W"   , matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
     AliMaterial(++matId,"Al"  ,aAl  ,zAl  ,dAl  ,radAl  ,absAl  );  AliMedium(kAl  ,"Al"  , matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
     AliMaterial(++matId,"Ar"  ,aAr  ,zAr  ,dAr  ,radAr  ,absAr  );  AliMedium(kAr  ,"Ar"  , matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
-
-    InitProperties();
         
 }//void AliHMPID::CreateMaterials()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDv3::InitProperties()
-{
+//void AliHMPIDv3::InitProperties()
+//{
 /*
 * HMPID
 * ====
@@ -155,7 +155,7 @@ HMPID  9  1.e-5 1.e-5 1.e-4  -1.   1.e-4 -1.   -1.    5.e-5 5.e-5 -1.  -1   -1
 HMPID  6  1.e-5 1.e-5 1.e-4  -1.   1.e-4 -1.   -1.    5.e-5 5.e-5 -1.  -1   -1   -1   -1   1    -1   1    -1   -1   -1   -1 
 * PCB backplane   (>  50 keV delta-electrons)
 HMPID 12  1.e-5 1.e-5 1.e-4  -1.   1.e-4 -1.   -1.    5.e-5 5.e-5 -1.  -1   -1   -1   -1   1    -1   1    -1   -1   -1   -1 
-*/
+
     Int_t *idtmed = fIdtmed->GetArray();
     Int_t imed;
     
@@ -236,7 +236,7 @@ HMPID 12  1.e-5 1.e-5 1.e-4  -1.   1.e-4 -1.   -1.    5.e-5 5.e-5 -1.  -1   -1
     gMC->Gstpar(idtmed[imed], "DRAY",1);    
     gMC->Gstpar(idtmed[imed], "LOSS",1);    
     
-}
+}*/
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 void AliHMPIDv3::CreateGeometry()
 {
@@ -508,14 +508,18 @@ void AliHMPIDv3::Init()
 // Arguments: none
 //   Returns: none      
   AliDebug(1,"Start v2 HMPID.");    
+  //InitProperties();
   AliDebug(1,"Stop v2 HMPID.");    
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 void AliHMPIDv3::DefineOpticalProperties() 
 {
+  AliDebug(1,"");    
+
 // Optical properties definition.
-  const Int_t kNbins=30;       //number of photon energy points
-  Float_t emin=5.5,emax=8.5;         //Photon energy range,[eV]
+  const Int_t kNbins=30;                 //number of photon energy points
+  Float_t emin=5.5,emax=8.5;             //Photon energy range,[eV]
+  Float_t deltaE = (emax - emin)/kNbins;
   Float_t aEckov [kNbins]; 
   Double_t dEckov [kNbins]; 
   Float_t aAbsRad[kNbins], aAbsWin[kNbins], aAbsGap[kNbins], aAbsMet[kNbins];
@@ -538,7 +542,7 @@ void AliHMPIDv3::DefineOpticalProperties()
   Bool_t isFlatIdx=title.Contains("FlatIdx"); 
   
   for(Int_t i=0;i<kNbins;i++){
-    Float_t eV=emin+0.1*i;  //Ckov energy in eV
+    Float_t eV=emin+deltaE*i;  //Ckov energy in eV
     aEckov [i] =1e-9*eV;    //Ckov energy in GeV
     dEckov [i] = aEckov[i];
     aAbsRad[i]=pRaAF->Eval(eV); (isFlatIdx)? aIdxRad[i]=1.292: aIdxRad[i]=pRaIF->Eval(eV,20);     
@@ -547,7 +551,7 @@ void AliHMPIDv3::DefineOpticalProperties()
     aQeAll[i] =1;                     //QE for all other materials except for PC must be 1.  
     aAbsMet[i] =0.0001;                aIdxMet[i]=0;                                             //metal ref idx must be 0 in order to reflect photon
                                        aIdxPc [i]=1;           aQePc [i]=pQeF->Eval(eV);         //PC ref idx must be 1 in order to apply photon to QE conversion 
-    dQePc [i]=pQeF->Eval(eV);
+    dQePc [i]= pQeF->Eval(eV);
     dReflMet[i] = 0.;     // no reflection on the surface of the pc (?)                                       
   }
   gMC->SetCerenkov((*fIdtmed)[kC6F14]    , kNbins, aEckov, aAbsRad  , aQeAll , aIdxRad );    
@@ -559,10 +563,11 @@ void AliHMPIDv3::DefineOpticalProperties()
   gMC->SetCerenkov((*fIdtmed)[kAl]       , kNbins, aEckov, aAbsMet  , aQeAll , aIdxMet );    
 
   // Define a skin surface for the photocatode to enable 'detection' in G4
-  gMC->DefineOpSurface("surfPc", kGlisur /*kUnified*/,kDielectric_metal,kPolished, 0.);
-  gMC->SetMaterialProperty("surfPc", "EFFICIENCY", kNbins, dEckov, dQePc);
-  gMC->SetMaterialProperty("surfPc", "REFLECTIVITY", kNbins, dEckov, dReflMet);
-  gMC->SetSkinSurface("skinPc", "Rpc", "surfPc");
+  for(Int_t i=0; i<7; i++){
+  gMC->DefineOpSurface(Form("surfPc%i",i), kGlisur /*kUnified*/,kDielectric_metal,kPolished, 0.);
+  gMC->SetMaterialProperty(Form("surfPc%i",i), "EFFICIENCY", kNbins, dEckov, dQePc);
+  gMC->SetMaterialProperty(Form("surfPc%i",i), "REFLECTIVITY", kNbins, dEckov, dReflMet);
+  gMC->SetSkinSurface(Form("skinPc%i",i), Form("Hpad%i",i),Form("surfPc%i",i)); }
 
   delete pRaAF;delete pWiAF;delete pGaAF; delete pRaIF; delete pWiIF; delete pGaIF; delete pQeF;
 }
@@ -601,7 +606,7 @@ void AliHMPIDv3::GenFee(Float_t qtot)
     gMC->GetRandom()->RndmArray(2,ranf);    //Sample direction
     cthf=ranf[0]*2-1.0;
     if(cthf<0) continue;
-    sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
+    sthf = TMath::Sqrt((1. - cthf) * (1. + cthf));
     phif = ranf[1] * 2 * TMath::Pi();
     
     if(Double_t randomNumber=gMC->GetRandom()->Rndm()<=0.57)
@@ -741,7 +746,7 @@ Float_t AliHMPIDv3::Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
     //FORMULAE FROM HANDBOOK OF OPTICS, 33.23 OR
     //W.R. HUNTER, J.O.S.A. 54 (1964),15 , J.O.S.A. 55(1965),1197
 
-    Float_t sinin=TMath::Sqrt(1-pdoti*pdoti);
+    Float_t sinin=TMath::Sqrt((1.-pdoti)*(1.+pdoti));
     Float_t tanin=sinin/pdoti;
 
     Float_t c1=cn*cn-ck*ck-sinin*sinin;
@@ -829,12 +834,13 @@ void AliHMPIDv3::StepHistory()
   }
 
   TString flag="fanny combination";
-  if(gMC->IsTrackAlive())
-      if(gMC->IsTrackEntering())      flag="enters to";
-      else if(gMC->IsTrackExiting())  flag="exits from";
-      else if(gMC->IsTrackInside())   flag="inside";
-  else
-      if(gMC->IsTrackStop())          flag="stoped in";        
+  if(gMC->IsTrackAlive()) {
+    if(gMC->IsTrackEntering())      flag="enters to";
+    else if(gMC->IsTrackExiting())  flag="exits from";
+    else if(gMC->IsTrackInside())   flag="inside";
+  } else {
+    if(gMC->IsTrackStop())          flag="stopped in";
+  }
   
   Int_t vid=0,copy=0;
   TString path=gMC->CurrentVolName(); path.Prepend("-");path.Prepend(gMC->CurrentVolOffName(1));//current volume and his mother are always there
@@ -888,9 +894,10 @@ void AliHMPIDv3::StepManager()
       Int_t   pid=     gMC->TrackPid();                                                          //take PID
       Float_t etot=    gMC->Etot();                                                              //total hpoton energy, [GeV] 
       Double_t x[3];   gMC->TrackPosition(x[0],x[1],x[2]);                                       //take MARS position at entrance to PC
-    TString tmpname = volname; tmpname.Remove(0,4); Int_t idch = tmpname.Atoi();                 //retrieve the chamber number
-    Float_t xl,yl;   AliHMPIDParam::Instance()->Mars2Lors(idch,x,xl,yl);                         //take LORS position 
-      new((*fHits)[fNhits++])AliHMPIDHit(idch,etot,pid,tid,xl,yl,x);                             //HIT for photon, position at P, etot will be set to Q
+      Float_t hitTime= (Float_t)gMC->TrackTime();                                                //hit formation time       
+      TString tmpname = volname; tmpname.Remove(0,4); Int_t idch = tmpname.Atoi();               //retrieve the chamber number
+      Float_t xl,yl;   AliHMPIDParam::Instance()->Mars2Lors(idch,x,xl,yl);                       //take LORS position 
+      new((*fHits)[fNhits++])AliHMPIDHit(idch,etot,pid,tid,xl,yl,hitTime,x);                             //HIT for photon, position at P, etot will be set to Q
       if(fDoFeed) GenFee(etot);                                                                  //generate feedback photons etot is modified in hit ctor to Q of hit
     }//photon hit PC and DE >0 
   }//photon hit PC
@@ -913,13 +920,14 @@ void AliHMPIDv3::StepManager()
       Int_t tid=          gMC->GetStack()->GetCurrentTrackNumber();                               //take TID
       Int_t pid=          gMC->TrackPid();                                                        //take PID
       Double_t out[3];    gMC->TrackPosition(out[0],out[1],out[2]);                               //take MARS position at exit
+      Float_t hitTime= (Float_t)gMC->TrackTime();                                                         //hit formation time       
       out[0]=0.5*(out[0]+in[0]);                                                                  //
       out[1]=0.5*(out[1]+in[1]);                                                                  //take hit position at the anod plane
       out[2]=0.5*(out[2]+in[2]);
       TString tmpname = volname;  tmpname.Remove(0,4);  Int_t idch = tmpname.Atoi();              //retrieve the chamber number
       Float_t xl,yl;AliHMPIDParam::Instance()->Mars2Lors(idch,out,xl,yl);                         //take LORS position
       if(eloss>0) {
-        new((*fHits)[fNhits++])AliHMPIDHit(idch,eloss,pid,tid,xl,yl,out);                           //HIT for MIP, position near anod plane, eloss will be set to Q 
+        new((*fHits)[fNhits++])AliHMPIDHit(idch,eloss,pid,tid,xl,yl,hitTime,out);                           //HIT for MIP, position near anod plane, eloss will be set to Q 
         if(fDoFeed) GenFee(eloss);                                                                  //generate feedback photons 
       }
     }else                                                                                         //just going inside
@@ -1230,7 +1238,7 @@ return cradle;
 }//CreateCradle()
 
 
-TGeoVolume * AliHMPIDv3::CradleBaseVolume(TGeoMedium *med, Double_t l[7],char *name)
+TGeoVolume * AliHMPIDv3::CradleBaseVolume(TGeoMedium *med, Double_t l[7],const char *name)
 {
 /*
 The trapezoid is build in the xy plane