Move HMPID QA code to PWG1.
[u/mrichter/AliRoot.git] / HMPID / AliHMPIDv1.cxx
index 58537bb..d0db224 100644 (file)
@@ -18,6 +18,7 @@
 #include "AliHMPIDParam.h"  //StepManager()
 #include "AliHMPIDHit.h"    //Hits2SDigs(),StepManager()
 #include "AliHMPIDDigit.h"  //Digits2Raw(), Raw2SDigits()
+#include "AliHMPIDRawStream.h"  //Digits2Raw(), Raw2SDigits()
 #include "AliRawReader.h"  //Raw2SDigits()
 #include <TVirtualMC.h>    //StepManager() for gMC
 #include <TPDGCode.h>      //StepHistory() 
 #include <AliMC.h>            //StepManager()      
 #include <AliRun.h>           //CreateMaterials()    
 #include <AliMagF.h>          //CreateMaterials()
-#include <TGeoManager.h>      //CreateGeometry()
+//#include <TGeoManager.h>      //CreateGeometry()
+#include <AliCDBEntry.h>      //CreateMaterials()
+#include <AliCDBManager.h>    //CreateMaterials()
 #include <TF1.h>              //DefineOpticalProperties()
 #include <TF2.h>              //DefineOpticalProperties()
+#include <TGeoGlobalMagField.h>
 #include <TLorentzVector.h>   //IsLostByFresnel() 
-#include <AliCDBManager.h>    //CreateMaterials()
-#include <AliCDBEntry.h>      //CreateMaterials()
  
 ClassImp(AliHMPIDv1)    
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -65,8 +67,8 @@ void AliHMPIDv1::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;                    //mas step allowed [cm]
@@ -206,7 +208,7 @@ void AliHMPIDv1::DefineOpticalProperties()
   Float_t                                                    aQeAll [kNbins], aQePc [kNbins];
   Double_t dReflMet[kNbins], dQePc[kNbins];
 
-  TF2 *pRaIF=new TF2("HidxRad","sqrt(1+0.554*(1239.84/x)^2/((1239.84/x)^2-5796)-0.0005*(y-20))"                                       ,emin,emax,0,50); //DiMauro mail temp 0-50 degrees C
+  TF2 *pRaIF=new TF2("HidxRad","sqrt(1+0.554*(1239.84/x)^2/((1239.84/x)^2-5769)-0.0005*(y-20))"                                       ,emin,emax,0,50); //DiMauro mail temp 0-50 degrees C
   TF1 *pWiIF=new TF1("HidxWin","sqrt(1+46.411/(10.666*10.666-x*x)+228.71/(18.125*18.125-x*x))"                                        ,emin,emax);      //SiO2 idx TDR p.35
   TF1 *pGaIF=new TF1("HidxGap","1+0.12489e-6/(2.62e-4 - x*x/1239.84/1239.84)"                                                         ,emin,emax);      //?????? from where  
 
@@ -281,7 +283,7 @@ void AliHMPIDv1::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)
@@ -388,7 +390,9 @@ void AliHMPIDv1::Digits2Raw()
   }
   treeD->GetEntry(0);
   
-  AliHMPIDDigit::WriteRaw(DigLst());
+  //AliHMPIDDigit::WriteRaw(DigLst());
+   AliHMPIDRawStream *pRS=0x0;
+   pRS->WriteRaw(DigLst());
     
   GetLoader()->UnloadDigits();
   AliDebug(1,"Stop.");      
@@ -419,7 +423,7 @@ Float_t AliHMPIDv1::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;
@@ -463,18 +467,20 @@ Bool_t AliHMPIDv1::Raw2SDigits(AliRawReader *pRR)
 // Interface methode ivoked from AliSimulation to create a list of sdigits from raw digits. Events loop is done in AliSimulation
 // Arguments: pRR- raw reader 
 //   Returns: kTRUE on success (currently ignored in AliSimulation::ConvertRaw2SDigits())      
-  AliHMPIDDigit sdi; //tmp sdigit, raw digit will be converted to it
+  //AliHMPIDDigit sdi; //tmp sdigit, raw digit will be converted to it
   
   if(!GetLoader()->TreeS()) {MakeTree("S");  MakeBranch("S");}
     
   TClonesArray *pSdiLst=SdiLst(); Int_t iSdiCnt=0; //tmp list of sdigits for all chambers
-  pRR->Select("HMPID",0,13);//select all HMPID DDL files
-  UInt_t w32=0;
-  while(pRR->ReadNextInt(w32)){//raw records loop (in selected DDL files)
-    UInt_t ddl=pRR->GetDDLID(); //returns 0,1,2 ... 13
-    if(!sdi.Raw(ddl,w32,pRR)) continue;  
-    new((*pSdiLst)[iSdiCnt++]) AliHMPIDDigit(sdi); //add this digit to the tmp list
-  }//raw records loop
+  AliHMPIDRawStream stream(pRR);
+  while(stream.Next())
+  {
+    for(Int_t iPad=0;iPad<stream.GetNPads();iPad++) {
+      AliHMPIDDigit sdi(stream.GetPadArray()[iPad],stream.GetChargeArray()[iPad]);
+      new((*pSdiLst)[iSdiCnt++]) AliHMPIDDigit(sdi); //add this digit to the tmp list
+    }
+  }
+  
   GetLoader()->TreeS()->Fill(); GetLoader()->WriteSDigits("OVERWRITE");//write out sdigits
   SdiReset();
   return kTRUE;
@@ -503,12 +509,13 @@ void AliHMPIDv1::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
@@ -547,8 +554,9 @@ void AliHMPIDv1::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
+      Float_t  hitTime=(Float_t)gMC->TrackTime();                                                //hit formation time
       Float_t xl,yl;   AliHMPIDParam::Instance()->Mars2Lors(copy,x,xl,yl);                       //take LORS position
-      new((*fHits)[fNhits++])AliHMPIDHit(copy,etot,pid,tid,xl,yl,x);                             //HIT for photon, position at P, etot will be set to Q
+      new((*fHits)[fNhits++])AliHMPIDHit(copy,etot,pid,tid,xl,yl,hitTime,x);                     //HIT for photon, position at P, etot will be set to Q
       GenFee(etot);                                                                              //generate feedback photons etot is modified in hit ctor to Q of hit
     }//photon hit PC and DE >0 
   }//photon hit PC
@@ -566,11 +574,12 @@ void AliHMPIDv1::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]);                                                                  //>
       Float_t xl,yl;AliHMPIDParam::Instance()->Mars2Lors(copy,out,xl,yl);                         //take LORS position
-      new((*fHits)[fNhits++])AliHMPIDHit(copy,eloss,pid,tid,xl,yl,out);                           //HIT for MIP, position near anod plane, eloss will be set to Q 
+      new((*fHits)[fNhits++])AliHMPIDHit(copy,eloss,pid,tid,xl,yl,hitTime,out);                   //HIT for MIP, position near anod plane, eloss will be set to Q 
       GenFee(eloss);                                                                              //generate feedback photons 
     }else                                                                                         //just going inside
       eloss          += gMC->Edep();                                                              //collect this step eloss