// GetLoader()->UnloadHits();
// return 0;
// }
-//__________________________________________________________________________________________________
-void AliHMPID::HitPrint(Int_t iEvtN)const
-{
-//Prints a list of HMPID hits for a given event. Default is event number 0.
- if(GetLoader()->GetRunLoader()->GetEvent(iEvtN)) return;
- AliInfo(Form("List of HMPID hits for event %i",iEvtN));
- if(GetLoader()->LoadHits()) return;
-
- Int_t iTotalHits=0;
- for(Int_t iPrimN=0;iPrimN<GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop
- GetLoader()->TreeH()->GetEntry(iPrimN);
- Hits()->Print();
- iTotalHits+=Hits()->GetEntries();
- }
- GetLoader()->UnloadHits();
- AliInfo(Form("totally %i hits",iTotalHits));
-}
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPID::SdiPrint(Int_t iEvt)const
-{
-//prints a list of HMPID sdigits for a given event
- if(GetLoader()->GetRunLoader()->GetEvent(iEvt)) return;
- Info("PrintSDigits","List of HMPID sdigits for event %i",iEvt);
- if(GetLoader()->LoadSDigits()) return;
-
- GetLoader()->TreeS()->GetEntry(0);
- SdiLst()->Print();
- GetLoader()->UnloadSDigits();
- Printf("totally %i sdigits",SdiLst()->GetEntries());
-}
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPID::DigPrint(Int_t iEvt)const
-{
-//prints a list of HMPID digits for a given event
- if(GetLoader()->GetRunLoader()->GetEvent(iEvt)) return;
- Printf("List of HMPID digits for event %i",iEvt);
- if(GetLoader()->LoadDigits()) return;
-
- GetLoader()->TreeD()->GetEntry(0);
- DigLst()->Print();
- Int_t totDigs=0;
- for(Int_t i=0;i<7;i++) {totDigs+=DigLst(i)->GetEntries();}
- Printf("totally %i Digits",totDigs);
- GetLoader()->UnloadDigits();
-}
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPID::CluPrint(Int_t iEvtN)const
-{
-//prints a list of HMPID clusters for a given event
- Printf("List of HMPID clusters for event %i",iEvtN);
- if(GetLoader()->GetRunLoader()->GetEvent(iEvtN)) return;
- if(GetLoader()->LoadRecPoints()) return;
-
- Int_t iCluCnt=0;
- GetLoader()->TreeR()->GetEntry(0);
- for(Int_t iCh=0;iCh<7;iCh++){
- TClonesArray *pCluLst=(TClonesArray*)fClu->At(iCh); iCluCnt+=pCluLst->GetEntries(); pCluLst->Print();
- }
- GetLoader()->UnloadRecPoints();
- Printf("totally %i clusters for event %i",iCluCnt,iEvtN);
-}
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void SetTreeAddress ( ); //from AliModule invoked from AliRun::GetEvent(), AliLoader::SetTAddrInDet()
virtual void StepManager ( )=0; //from AliModule invoked from AliMC
//private part +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
- void HitCreate( ) {if(fHits)return; fHits=new TClonesArray("AliHMPIDHit"); fNhits=0; }//create hits list
- void HitPrint (Int_t evt)const; //print hits list
+ void HitCreate( ) {if(fHits)return; fHits=new TClonesArray("AliHMPIDHit"); fNhits=0; }//create hits list
TClonesArray* SdiLst ( )const{return fSdi; }//get sdigits list
- void SdiCreate( ) {if(fSdi)return; fSdi=new TClonesArray("AliHMPIDDigit"); }//create sdigits list
+ void SdiCreate( ) {if(fSdi)return; fSdi=new TClonesArray("AliHMPIDDigit"); }//create sdigits list
void SdiReset ( ) {if(fSdi) fSdi ->Clear(); }//clean sdigits list
- void SdiPrint (Int_t evt)const; //print sdigits
TObjArray* DigLst ( )const{return fDig; }//get digits list for all chambers
TClonesArray* DigLst (Int_t c )const{return fDig ? (TClonesArray *)fDig->At(c):0; }//get digits list for chamber
if (fDig) return; //PH do not recreate existing containers
fDig=new TObjArray(7);for(Int_t i=0;i<7;i++)fDig->AddAt(new TClonesArray("AliHMPIDDigit"),i); }//create digits list
void DigReset ( ) {if(fDig)for(int i=0;i<7;i++)fDig->At(i)->Clear(); }//clean digits list
- void DigPrint (Int_t evt)const; //print digits
TObjArray* CluLst ( )const{return fClu; }//get clusters list for all chambers
TClonesArray* CluLst (Int_t c )const{return fClu ? (TClonesArray *)fClu->At(c):0; }//get clusters list for chamber
if (fClu) return; //PH do not recreate existing containers
fClu=new TObjArray(7); for(Int_t i=0;i<7;i++)fClu->AddAt(new TClonesArray("AliHMPIDCluster"),i); }//create clusters list
void CluReset ( ) {if(fClu)for(int i=0;i<7;i++)fClu->At(i)->Clear(); }//clean clusters list
- void CluPrint (Int_t evt)const; //print clusters list
protected:
TClonesArray *fSdi; //! list of sdigits
TObjArray *fDig; //! each chamber holds it's one list of digits
TObjArray *fClu; //! each chamber holds it's one list of clusters
- private:
+private:
AliHMPID(const AliHMPID &rich );
AliHMPID& operator=(const AliHMPID&);
// Correction of cluster x position due to sinoid, see HMPID TDR page 30
// Arguments: none
// Returns: none
- AliHMPIDDigit dig;dig.Manual1(Ch(),fX,fY); //tmp digit to get it center
- Float_t x=fX-dig.LorsX();
+ Int_t pc,px,py;
+ AliHMPIDDigit::Lors2Pad(fX,fY,pc,px,py); //tmp digit to get it center
+ Float_t x=fX-AliHMPIDDigit::LorsX(pc,px); //diff between cluster x and center of the pad contaning this cluster
fX+=3.31267e-2*TMath::Sin(2*TMath::Pi()/0.8*x)-2.66575e-3*TMath::Sin(4*TMath::Pi()/0.8*x)+2.80553e-3*TMath::Sin(6*TMath::Pi()/0.8*x)+0.0070;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
}
Double_t ratio=0;
if(Q()>0&&QRaw()>0) ratio = Q()/QRaw()*100;
- Printf("%sCLU:(%7.3f,%7.3f) Qfitted=%8.3f QRaw=%8.3f(%3.0f%%) ch=%i Npads=%2i DimBox %i NlocMax %i Chi2 %f %s",
+ Printf("%sCLU: ch=%i (%7.3f,%7.3f) Q=%8.3f Qraw=%8.3f(%3.0f%%) Size=%2i DimBox=%i LocMax=%i Chi2=%7.3f %s",
opt,X(),Y(),Q(),QRaw(),ratio,Ch(),Size(),fBox,fNlocMax,fChi2,status);
if(fDigs) fDigs->Print();
}//Print()
// * provided "as is" without express or implied warranty. *
// **************************************************************************
-#include "AliHMPIDDigit.h" //class header
-#include <TClonesArray.h> //Hit2Sdi()
-#include <TBox.h> //Draw()
+#include "AliHMPIDDigit.h" //class header
+#include <TClonesArray.h> //WriteRaw()
+#include <TBox.h> //Draw()
+#include <AliRawDataHeader.h> //WriteRaw()
+#include <AliDAQ.h> //WriteRaw()
+#include <Riostream.h> //WriteRaw()
ClassImp(AliHMPIDDigit)
-const Float_t AliHMPIDDigit::fMinPcX[]={0,SizePcX() + SizeDead(),0,SizePcX() + SizeDead(), 0,SizePcX() + SizeDead()};
-const Float_t AliHMPIDDigit::fMinPcY[]={0, 0,SizePcY()+SizeDead(),SizePcY() + SizeDead(),2*(SizePcY()+SizeDead()),2*(SizePcY()+SizeDead())};
+const Float_t AliHMPIDDigit::fMinPcX[]={ 0.00 , 66.60 , 0.00 , 66.60 , 0.00 , 66.60};
+const Float_t AliHMPIDDigit::fMaxPcX[]={64.00 , 130.60 , 64.00 , 130.60 , 64.00 , 130.60};
-const Float_t AliHMPIDDigit::fMaxPcX[]={SizePcX(),SizeAllX(),SizePcX(),SizeAllX(),SizePcX(),SizeAllX()};
-const Float_t AliHMPIDDigit::fMaxPcY[]={SizePcY(),SizePcY(),SizeAllY() - SizePcY(),SizeAllY() - SizePcY(),SizeAllY(),SizeAllY()};
+const Float_t AliHMPIDDigit::fMinPcY[]={ 0.00 , 42.92 , 85.84, 0.00 , 42.92 , 85.84};
+const Float_t AliHMPIDDigit::fMaxPcY[]={40.32 , 83.24 , 126.16, 40.32 , 83.24 , 126.16};
/*
+
+ y6 ---------- ---------- 126.16
+ | | | |
+ | 4 | | 5 |
+ y5 ---------- ---------- 85.84
+
+ y4 ---------- ---------- 83.24
+ | | | |
+ | 2 | | 3 | view from electronics side
+ y3 ---------- ---------- 42.92
+
+ y2 ---------- ---------- 40.32
+ | | | |
+ | 0 | | 1 |
+ y1 ---------- ---------- 0
+ x1 x2 x3 x4
+
+ 0 64.00 66.60 130.60
Preface: all geometrical information (like left-right sides) is reported as seen from electronic side.
pad->Draw();
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDDigit::Hit2Sdi(AliHMPIDHit *pHit,TClonesArray *pSdiLst)
-{
-// Creates a list of sdigits out of provided hit
-// Arguments: pHit- hit
-// Returns: none
-
- Int_t iSdiCnt=pSdiLst->GetEntries(); //list of sdigits contains sdigits from previous ivocations of Hit2Sdi, do not override them
- AliHMPIDDigit dig;
- for(Int_t i=0;i<9;i++){ //affected pads loop
- dig.Set(pHit,i); //c,q,tid,x,y create tmp sdigit for pad i around hit position
- if(dig.PadPcX()==-1) continue;
- new((*pSdiLst)[iSdiCnt++]) AliHMPIDDigit(dig);
- }
-}//Hit2Sdi()
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDDigit::Print(Option_t*)const
+void AliHMPIDDigit::Print(Option_t *opt)const
{
// Print current digit
// Arguments: option string not used
// Returns: none
- UInt_t w32; Raw(w32);
- Printf("DIG:(%7.3f,%7.3f) Q=%8.3f (ch=%1i,pc=%1i,x=%2i,y=%2i) TID=(%5i,%5i,%5i) ddl=%i raw=0x%x (r=%2i,d=%2i,a=%2i) %s",
- LorsX(),LorsY(),Q(), A2C(fPad),A2P(fPad),A2X(fPad),A2Y(fPad),
- fTracks[0],fTracks[1],fTracks[2],DdlIdx(),w32,Row(),Dilogic(),Addr(),
- (IsOverTh(Q()))?"":"!!!");
-}
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDDigit::PrintSize()
-{
-// Print all segmentaion related sizes
-// Arguments: none
-// Returns: none
- Printf("-->pad =(%6.2f,%6.2f) cm dead zone %.2f cm\n"
- "-->PC =(%6.2f,%6.2f) cm (%3i,%3i) pads\n"
- "-->all PCs=(%6.2f,%6.2f) cm (%3i,%3i) pads",
- SizePadX(),SizePadY(),SizeDead(),
- SizePcX() ,SizePcY() ,kPadPcX ,kPadPcY,
- SizeAllX(),SizeAllY(),kPadAllX,kPadAllY);
+ UInt_t w32; Int_t ddl,r,d,a;
+ Raw(w32,ddl,r,d,a);
+ Printf("%sDIG:(ch=%1i,pc=%1i,x=%2i,y=%2i) (%7.3f,%7.3f) Q=%8.3f TID=(%5i,%5i,%5i) raw=0x%x (ddl=%2i,r=%2i,d=%2i,a=%2i) %s",
+ opt, A2C(fPad),A2P(fPad),A2X(fPad),A2Y(fPad),LorsX(),LorsY(), Q(), fTracks[0],fTracks[1],fTracks[2],w32,ddl,r,d,a, (IsOverTh(Q()))?"":"below thr");
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDDigit::Test()
+void AliHMPIDDigit::WriteRaw(TObjArray *pDigAll)
{
- AliHMPIDDigit d1,d2; Int_t ddl;UInt_t w32;
- for(Int_t i=0;i<10000;i++){
- Int_t ch=Int_t(gRandom->Rndm()*7);
- Int_t pc=Int_t(gRandom->Rndm()*6);
- Int_t px=Int_t(gRandom->Rndm()*80);
- Int_t py=Int_t(gRandom->Rndm()*48);
- d1.Manual2(ch,pc,px,py);
- ddl=d1.Raw(w32); d2.Raw(ddl,w32);
- if(d1.Compare(&d2)) Printf("Problem!!!");
- }
- Printf("OK");
-}//Test()
+// Write a list of digits for a given chamber in raw data stream
+// Arguments: pDigAll- list of digits
+// Returns: none
+ ofstream ddlL,ddlR; //output streams, 2 per chamber
+ Int_t cntL=0,cntR=0; //data words counters for DDLs
+ AliRawDataHeader header; header.SetAttribute(0); //empty DDL header
+
+ UInt_t w32=0; Int_t ddl,r,d,a; //32 bits data word
+ for(Int_t iCh=kMinCh;iCh<=kMaxCh;iCh++){//chambers loop
+ ddlL.open(AliDAQ::DdlFileName("HMPID",2*iCh));
+ ddlR.open(AliDAQ::DdlFileName("HMPID",2*iCh+1)); //open both DDL of this chamber in parallel
+ ddlL.write((char*)&header,sizeof(header)); //write dummy header as place holder, actual
+ ddlR.write((char*)&header,sizeof(header)); //will be rewritten later when total size of DDL is known
+
+ TClonesArray *pDigCh=(TClonesArray *)pDigAll->At(iCh); //list of digits for current chamber
+ for(Int_t iDig=0;iDig<pDigCh->GetEntriesFast();iDig++){//digits loop
+ AliHMPIDDigit *pDig=(AliHMPIDDigit*)pDigCh->At(iDig);
+ pDig->Raw(w32,ddl,r,d,a);
+ if(ddl%2){
+ ddlL.write((char*)&w32,sizeof(w32)); cntL++;
+ }else{
+ ddlR.write((char*)&w32,sizeof(w32)); cntR++;
+ }
+ }//digits loop
+
+ header.fSize=sizeof(header)+cntL*sizeof(w32); ddlL.seekp(0); ddlL.write((char*)&header,sizeof(header)); ddlL.close(); //rewrite header with size set to
+ header.fSize=sizeof(header)+cntR*sizeof(w32); ddlR.seekp(0); ddlR.write((char*)&header,sizeof(header)); ddlR.close(); //number of bytes and close file
+ }//chambers loop
+}//WriteRaw()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
#include <TMath.h> //Mathieson()
#include <TRandom.h> //IsOverTh()
#include <AliBitPacking.h> //Raw()
-#include "AliHMPIDHit.h" //Set()
+
class TClonesArray; //Hit2Sdi()
class AliHMPIDDigit :public AliDigit //TObject-AliDigit-AliHMPIDDigit
{
public:
- enum EAbsPad {kChAbs=100000000,kPcAbs=1000000,kPadAbsX=1000,kPadAbsY=1}; //absolute pad number structure
- enum ERawData{kNddls=14}; //RAW data structure
- enum EPadData{kPcX=2,kPcY=3,kPadPcX=80,kPadPcY=48,kPadAllX=kPadPcX*kPcX,kPadAllY=kPadPcY*kPcY,kPcAll=kPcX*kPcY,kPadAll=kPadAllX*kPadAllY}; //Segmentation structure
+ enum EChamberData{kMinCh=0,kMaxCh=6,kMinPc=0,kMaxPc=5};
+ enum EPadxData{kPadPcX=80,kMinPx=0,kMaxPx=kPadPcX-1,kMaxPadChX=159}; //Segmentation structure along x
+ enum EPadyData{kPadPcY=48,kMinPy=0,kMaxPy=kPadPcY-1,kMaxPadChY=143}; //Segmentation structure along y
//ctor&dtor
- AliHMPIDDigit( ):AliDigit( ),fPad(Abs(-1,-1,-1,-1)),fQ(0) {} //default ctor
- AliHMPIDDigit(Int_t pad,Int_t q,Int_t *t ):AliDigit(t),fPad(pad ),fQ(q ) {} //ctor used in digitizer
- virtual ~AliHMPIDDigit( ) {} //dtor
+ AliHMPIDDigit( ):AliDigit( ),fPad(Abs(-1,-1,-1,-1)),fQ(-1) {} //default ctor
+ AliHMPIDDigit(Int_t pad,Int_t q,Int_t *t):AliDigit(t),fPad(pad ),fQ(q ) {} //digit ctor
+ AliHMPIDDigit(const AliHMPIDDigit &d ):AliDigit(d),fPad(d.fPad),fQ(d.fQ) {} //copy ctor
+ virtual ~AliHMPIDDigit() {} //dtor
//framework part
Bool_t IsSortable ( )const{return kTRUE;} //provision to use TObject::Sort()
inline Int_t Compare (const TObject *pObj )const; //provision to use TObject::Sort()
void Draw (Option_t *opt="" ); //TObject::Draw() overloaded
void Print (Option_t *opt="" )const; //TObject::Print() overloaded
//private part
- static Int_t Abs (Int_t c,Int_t s,Int_t x,Int_t y) {return c*kChAbs+s*kPcAbs+x*kPadAbsX+y*kPadAbsY; } //(ch,pc,padx,pady)-> abs pad
- static Int_t A2C (Int_t pad ) {return pad/kChAbs; } //abs pad -> chamber
- static Int_t A2P (Int_t pad ) {return pad%kChAbs/kPcAbs; } //abs pad -> pc
- static Int_t A2X (Int_t pad ) {return pad%kPcAbs/kPadAbsX; } //abs pad -> pad X
- static Int_t A2Y (Int_t pad ) {return pad%kPadAbsX; } //abs pad -> pad Y
- Int_t Addr ( )const{Int_t map[6]={5,3,1,0,2,4};return map[A2Y(fPad)%6]+6*(A2X(fPad)%8);}//ADDRESS 0..47
- void AddTidOffset(Int_t offset ) {for (Int_t i=0; i<3; i++) if (fTracks[i]>0) fTracks[i]+=offset;}; //needed for merging
- Int_t Ch ( )const{return A2C(fPad); } //chamber number
- Int_t Dilogic ( )const{return 1+PadPcX()/8; } //DILOGIC# 1..10
- Int_t DdlIdx ( )const{return 2*Ch()+Pc()%2; } //DDL# 0..13
- Int_t DdlId ( )const{return (6<<8)+DdlIdx(); } //DDL ID 0x600..0x60d
- static void Hit2Sdi (AliHMPIDHit *pHit,TClonesArray*); //hit -> 9 sdigits
- static Bool_t IsOverTh (Float_t q ) {return q >= 4; } //is digit over threshold????
- static Bool_t IsInside (Float_t x,Float_t y ) {return x>0&&y>0&&x<SizeAllX()&&y<SizeAllY(); } //is point inside chamber boundary?
- inline static Bool_t IsInDead (Float_t x,Float_t y ); //is point in dead area?
- Float_t LorsX ( )const{return (PadPcX()+0.5)*SizePadX()+(Pc()%2)*(SizePcX()+SizeDead());} //center of the pad x, [cm]
+ static Int_t Abs (Int_t c,Int_t s,Int_t x,Int_t y) {return c*100000000+s*1000000+x*1000+y; } //(ch,pc,padx,pady)-> abs pad
+ static Int_t A2C (Int_t pad ) {return pad/100000000; } //abs pad -> chamber
+ static Int_t A2P (Int_t pad ) {return pad%100000000/1000000; } //abs pad -> pc
+ static Int_t A2X (Int_t pad ) {return pad%1000000/1000; } //abs pad -> pad X
+ static Int_t A2Y (Int_t pad ) {return pad%1000; } //abs pad -> pad Y
+ void AddTidOffset(Int_t offset ) {for (Int_t i=0; i<3; i++) if (fTracks[i]>0) fTracks[i]+=offset; } //needed for merging
+ Int_t Ch ( )const{return A2C(fPad); } //chamber number
+ static Bool_t IsOverTh (Float_t q ) {return q >= 4; } //is digit over threshold????
+ static Bool_t IsInside (Float_t x,Float_t y ) {return x>0&&y>0&&x<SizeAllX()&&y<SizeAllY(); } //is point inside chamber boundary?
+ Float_t LorsX ( )const{return LorsX(A2P(fPad),A2X(fPad)); } //center of the pad x, [cm]
static Float_t LorsX (Int_t pc,Int_t padx ) {return (padx +0.5)*SizePadX()+(pc %2)*(SizePcX()+SizeDead());} //center of the pad x, [cm]
- Float_t LorsY ( )const{return (PadPcY()+0.5)*SizePadY()+(Pc()/2)*(SizePcY()+SizeDead());} //center of the pad y, [cm]
+ Float_t LorsY ( )const{return LorsY(A2P(fPad),A2Y(fPad)); } //center of the pad y, [cm]
static Float_t LorsY (Int_t pc,Int_t pady ) {return (pady +0.5)*SizePadY()+(pc /2)*(SizePcY()+SizeDead());} //center of the pad y, [cm]
- void Manual1 (Int_t c,Float_t x,Float_t y ) {AliHMPIDHit h(c,200e-9,2212,3,x,y); Set(&h,0);} //manual from hit
- void Manual2 (Int_t c,Int_t p,Int_t x,Int_t y,Float_t q=0) {fPad=Abs(c,p,x,y);fQ=q;} //manual creation
inline Float_t Mathieson (Float_t x,Float_t y )const; //Mathieson distribution
Int_t PadPcX ( )const{return A2X(fPad);} //pad pc x # 0..79
Int_t PadPcY ( )const{return A2Y(fPad);} //pad pc y # 0..47
Int_t PadChY ( )const{return (Pc()/2)*kPadPcY+PadPcY();} //pad ch y # 0..143
Int_t Pad ( )const{return fPad;} //absolute id of this pad
Int_t Pc ( )const{return A2P(fPad);} //PC position number
- static void PrintSize ( ); //print all segmentation sizes
Float_t Q ( )const{return fQ;} //charge, [QDC]
- inline Int_t Raw ( UInt_t &raw32 )const; //digit->(ddl,raw32)
- inline void Raw (Int_t l,UInt_t raw32 ); //(ddl,raw32)->digit
- static Int_t Raw2Ch (UInt_t l ) {return l/2;} //ch=f(ddl)
- static Int_t Raw2Pc (UInt_t l,UInt_t r ) {r=(r-1)/8;return (l%2)?5-2*r:2*r;} //pc=f(ddl,r)
- static Int_t Raw2X ( UInt_t d,UInt_t a ) { return (d-1)*8+a/6;} //padx=f(d,a)
- static Int_t Raw2Y (UInt_t l,UInt_t r,UInt_t a ) {Int_t a2y[6]={3,2,4,1,5,0};r=(l%2)?(24-r):r-1;return 6*(r%8)+a2y[a%6];}//pady=f(ddl,r,a)
- Int_t Row ( )const{Int_t r=1+Pc()/2*8+PadPcY()/6; return (Pc()%2)? 25-r:r;} //row r=1..24
- inline Bool_t Set (AliHMPIDHit *pHit,Int_t pad ); //sdigit from hit in given pad
+ inline void Raw (UInt_t &w32,Int_t &ddl,Int_t &r,Int_t &d,Int_t &a)const; //digit->(w32,ddl,r,d,a)
+ inline void Raw (UInt_t w32,Int_t ddl ); //(w32,ddl)->digit
+ inline Bool_t Set (Int_t c,Int_t p,Int_t x,Int_t y,Float_t q=0,Int_t tid=0); //manual creation
+ static void WriteRaw (TObjArray *pDigLst ); //write as raw stream
- static Float_t SizeAllX ( ) {return SizePadX()*kPadAllX+SizeDead();} //all PCs size x, [cm]
- static Float_t SizeAllY ( ) {return SizePadY()*kPadAllY+2*SizeDead();} //all PCs size y, [cm]
- static Float_t SizeArea ( ) {return SizePcX()*SizePcY()*kPcAll;} //sence area, [cm^2]
+ static Float_t SizeAllX ( ) {return fMaxPcX[5];} //all PCs size x, [cm]
+ static Float_t SizeAllY ( ) {return fMaxPcY[5];} //all PCs size y, [cm]
+ static Float_t SizeArea ( ) {return SizePcX()*SizePcY()*(kMaxPc-kMinPc+1);} //sence area, [cm^2]
static Float_t SizeDead ( ) {return 2.6;} //dead zone size x, [cm]
static Float_t SizeGap ( ) {return 8; }
static Float_t SizePadX ( ) {return 0.8;} //pad size x, [cm]
static Float_t SizePadY ( ) {return 0.84;} //pad size y, [cm]
- static Float_t SizePcX ( ) {return SizePadX()*kPadPcX;} //PC size x, [cm]
- static Float_t SizePcY ( ) {return SizePadY()*kPadPcY;} //PC size y, [cm]
+ static Float_t SizePcX ( ) {return fMaxPcX[0];} //PC size x, [cm]
+ static Float_t SizePcY ( ) {return fMaxPcY[0];} //PC size y, [cm]
static Float_t SizeWin ( ) {return 0.5;} //Quartz window width
static Float_t SizeRad ( ) {return 1.5;} //Rad width
- static void Test ( ); //Test conversions
static const Float_t fMinPcX[6];
static const Float_t fMinPcY[6];
static const Float_t fMaxPcX[6];
static const Float_t fMaxPcY[6];
+ inline static Bool_t IsInDead(Float_t x,Float_t y ); //is point in dead area?
+ inline static void Lors2Pad(Float_t x,Float_t y,Int_t &pc,Int_t &px,Int_t &py); //(x,y)->(pc,px,py)
protected: //AliDigit has fTracks[3]
Int_t fPad; //absolute pad number
typedef AliHMPIDDigit AliRICHDigit; // for backward compatibility
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Bool_t AliHMPIDDigit::Set(AliHMPIDHit *pHit,Int_t pad)
+void AliHMPIDDigit::Lors2Pad(Float_t x,Float_t y,Int_t &pc,Int_t &px,Int_t &py)
{
-// Creates digit
-// Arguments: pHit- pointer to the hit
-// pad - for which pad to create
-// Returns: none
-
- fPad=Abs(-1,-1,-1,-1); fQ=0; //reset
- Int_t pc,px,py;
- Float_t x=pHit->LorsX(),y=pHit->LorsY();
-
+// Check the pad of given position
+// Arguments: x,y- position [cm] in LORS; pc,px,py- pad where to store the result
+// Returns: none
+ pc=px=py=-1;
if (x>= 0 && x<= SizePcX() ) {pc=0; px=Int_t( x / SizePadX());}//PC 0 or 2 or 4
else if(x>=SizePcX()+SizeDead() && x<= SizeAllX() ) {pc=1; px=Int_t((x- SizePcX()- SizeDead()) / SizePadX());}//PC 2 or 4 or 6
- else return kFALSE;
+ else return;
if (y>= 0 && y<= SizePcY() ) { py=Int_t( y / SizePadY());}//PC 0 or 1
else if(y>=SizePcY()+SizeDead() && y<=2*SizePcY()+SizeDead() ) {pc+=2;py=Int_t((y- SizePcY()- SizeDead()) / SizePadY());}//PC 2 or 3
else if(y>=SizeAllY()-SizePcY() && y<= SizeAllY() ) {pc+=4;py=Int_t((y-2*SizePcY()-2*SizeDead()) / SizePadY());}//PC 4 or 5
- else return kFALSE;
-
- switch(pad){
- case 8: px--;py++;break; case 1:py++;break; case 2:px++; py++;break;
-
- case 7: px--; break; case 0: break; case 3:px++; break;
-
- case 6: px--;py--;break; case 5:py--;break; case 4:px++; py--;break;
- }
- if(px<0 || px>=kPadPcX) return kFALSE;
- if(py<0 || py>=kPadPcY) return kFALSE;
- fPad=Abs(pHit->Ch(),pc,px,py);
- fQ=pHit->Q()*Mathieson(x,y);
- fTracks[0]=pHit->Tid();
- return kTRUE;
-}//Set()
+ else return;
+}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Int_t AliHMPIDDigit::Compare(const TObject *pObj) const
{
return 4*k4*(TMath::ATan(ux2)-TMath::ATan(ux1))*k4*(TMath::ATan(uy2)-TMath::ATan(uy1));
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Int_t AliHMPIDDigit::Raw(UInt_t &w32)const
+void AliHMPIDDigit::Raw(UInt_t &w32,Int_t &ddl,Int_t &r,Int_t &d,Int_t &a)const
{
// Convert digit structure to raw word format
-// Arguments: 32 bits raw word to fill
-// Returns: DDL ID where to write this digit
- w32=0;
- AliBitPacking::PackWord((UInt_t)fQ ,w32, 0,11); // 0000 0rrr rrdd ddaa aaaa qqqq qqqq qqqq Qdc bits (00..11) counts (0..4095)
- AliBitPacking::PackWord( Addr() ,w32,12,17); // 3322 2222 2222 1111 1111 1000 0000 0000 DILOGIC address bits (12..17) counts (0..47)
- AliBitPacking::PackWord( Dilogic(),w32,18,21); // 1098 7654 3210 9876 5432 1098 7654 3210 DILOGIC number bits (18..21) counts (1..10)
- AliBitPacking::PackWord( Row() ,w32,22,26); // Row number bits (22..26) counts (1..24)
- return DdlIdx(); //ddl 0..13 where to write this digit
+// Arguments: w32,ddl,r,d,a where to write the results
+// Returns: none
+ Int_t y2a[6]={5,3,1,0,2,4};
+
+ ddl=2*Ch()+Pc()%2; //DDL# 0..13
+ Int_t tmp=1+Pc()/2*8+PadPcY()/6; r=(Pc()%2)? 25-tmp:tmp; //row r=1..24
+ d=1+PadPcX()/8; //DILOGIC# 1..10
+ a=y2a[PadPcY()%6]+6*(PadPcX()%8); //ADDRESS 0..47
+
+ w32=0;
+ AliBitPacking::PackWord((UInt_t)fQ,w32, 0,11); // 0000 0rrr rrdd ddaa aaaa qqqq qqqq qqqq Qdc bits (00..11) counts (0..4095)
+ AliBitPacking::PackWord( a ,w32,12,17); // 3322 2222 2222 1111 1111 1000 0000 0000 DILOGIC address bits (12..17) counts (0..47)
+ AliBitPacking::PackWord( d ,w32,18,21); // 1098 7654 3210 9876 5432 1098 7654 3210 DILOGIC number bits (18..21) counts (1..10)
+ AliBitPacking::PackWord( r ,w32,22,26); // Row number bits (22..26) counts (1..24)
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDDigit::Raw(Int_t ddl,UInt_t w32)
+void AliHMPIDDigit::Raw(UInt_t w32,Int_t ddl)
{
// Converts a given raw data word to a digit
// Arguments: w32 - 32 bits raw data word
// ddl - DDL idx 0 1 2 3 4 ... 13
// Returns: none
- fQ = AliBitPacking::UnpackWord(w32, 0,11); // 0000 0rrr rrdd ddaa aaaa qqqq qqqq qqqq Qdc bits (00..11) counts (0..4095)
- UInt_t a = AliBitPacking::UnpackWord(w32,12,17); // 3322 2222 2222 1111 1111 1000 0000 0000 DILOGIC address bits (12..17) counts (0..47)
- UInt_t d = AliBitPacking::UnpackWord(w32,18,21); // 1098 7654 3210 9876 5432 1098 7654 3210 DILOGIC number bits (18..21) counts (1..10)
- UInt_t r = AliBitPacking::UnpackWord(w32,22,26); // Row number bits (22..26) counts (1..24)
- fPad=Abs(Raw2Ch(ddl),Raw2Pc(ddl,r),Raw2X(d,a),Raw2Y(ddl,r,a));
+ Int_t a2y[6]={3,2,4,1,5,0};//pady for a given address (for single DILOGIC chip)
+ Int_t r = AliBitPacking::UnpackWord(w32,22,26); assert(1<=r&&r<=24); // Row number (1..24)
+ Int_t d = AliBitPacking::UnpackWord(w32,18,21); assert(1<=d&&d<=10); // 3322 2222 2222 1111 1111 1000 0000 0000 DILOGIC number (1..10)
+ Int_t a = AliBitPacking::UnpackWord(w32,12,17); assert(0<=a&&a<=47); // 1098 7654 3210 9876 5432 1098 7654 3210 DILOGIC address (0..47)
+ Int_t q = AliBitPacking::UnpackWord(w32, 0,11); assert(0<=q&&q<=4095); // 0000 0rrr rrdd ddaa aaaa qqqq qqqq qqqq Qdc (0..4095)
+ Int_t ch=ddl/2;
+ Int_t tmp=(r-1)/8; Int_t pc=(ddl%2)? 5-2*tmp:2*tmp;
+ Int_t px=(d-1)*8+a/6;
+ tmp=(ddl%2)?(24-r):r-1; Int_t py=6*(tmp%8)+a2y[a%6];
+ fPad=Abs(ch,pc,px,py);fQ=q;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Bool_t AliHMPIDDigit::Set(Int_t ch,Int_t pc,Int_t px,Int_t py,Float_t qdc,Int_t tid)
+{
+// Manual creation of digit
+// Arguments: ch,pc,px,py,qdc,tid
+// Returns: none
+ if(px<kMinPx || px>kMaxPx) return kTRUE;
+ if(py<kMinPy || py>kMaxPy) return kTRUE;
+
+ fPad=Abs(ch,pc,px,py);fQ=qdc;fTracks[0]=tid;
+ return kFALSE;
+}
#endif
// **************************************************************************
#include "AliHMPIDHit.h" //class header
-#include "AliHMPIDDigit.h"
-#include <TPDGCode.h>
-#include <TMarker.h>
+#include <TPDGCode.h> //Draw() Print()
+#include <TMarker.h> //Draw()
+#include <TClonesArray.h> //Hit2Sdi()
ClassImp(AliHMPIDHit)
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
case 50000051: iMark=27; break;
default: iMark=26; break;
}
- TMarker *pMark=new TMarker(LorsX(),LorsY(),iMark); pMark->SetMarkerColor(kRed); pMark->Draw();
+ TMarker *pMark=new TMarker(fLx,fLy,iMark); pMark->SetMarkerColor(kRed); pMark->Draw();
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDHit::Print(Option_t*)const
+void AliHMPIDHit::Hit2Sdi(TClonesArray *pSdiLst,Int_t iHow)const
+{
+// Adds sdigits of this hit to the list
+// Arguments: pSdiLst- sigits list where to add new sdgits
+// iHow- how many pads to check
+// Returns: none
+ Int_t pc,px,py;
+ AliHMPIDDigit::Lors2Pad(fLx,fLy,pc,px,py); if(pc<0) return; //check if the hit in dead zone. Should never happen during trasport!
+
+ AliHMPIDDigit dig;
+ Int_t iSdiCnt=pSdiLst->GetEntries(); //list of sdigits contains sdigits from previous ivocations of Hit2Sdi, do not override them
+
+ for(Int_t i=-iHow;i<=iHow;i++){ //horizontal loop
+ for(Int_t j=-iHow;j<=iHow;j++){ //vertical loop
+ if(dig.Set(fCh,pc,px+i,py+j,fQ*dig.Mathieson(fLx,fLy),fTrack)) continue;
+ new((*pSdiLst)[iSdiCnt++]) AliHMPIDDigit(dig);
+ }
+ }
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void AliHMPIDHit::Print(Option_t *opt)const
{
//Print hit
char *sPart=Form("pid=%i",Pid());
case 50000051: sPart="feed";break;
}
- Printf("HIT:(%7.3f,%7.3f) Q=%8.3f ch=%i TID= %5i, MARS=(%7.2f,%7.2f,%7.2f) %s %s",
- LorsX(),LorsY(), Q(), Ch(), Tid(), X(), Y(), Z(), sPart,
+ Printf("%sHIT: ch=%i (%7.3f,%7.3f) Q=%8.3f TID= %5i, MARS=(%7.2f,%7.2f,%7.2f) %s %s",
+ opt, Ch(), fLx,fLy, fQ, fTrack, X(), Y(), Z(), sPart,
(AliHMPIDDigit::IsInDead(LorsX(),LorsY()))? "IN DEAD ZONE":"");
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include <AliHit.h> //base class
#include <TVector3.h> //ctor
#include <TRandom.h> //QdcTot()
-#include <TMath.h>
+#include "AliHMPIDDigit.h" //QdcTot()
class AliHMPIDHit : public AliHit // TObject-AliHit-AliHMPIDHit
{
public:
- AliHMPIDHit( ):AliHit( ),fCh(-1),fPid(-1),fQ(0),fLx(0),fLy(0){} //default ctor
- AliHMPIDHit(Int_t c,Float_t e,Int_t id,Int_t tn,Float_t x,Float_t y,const TVector3 &p):AliHit(0,tn),fCh(c ),fPid(id),fQ(0),fLx(x),fLy(y){QdcTot(e);fX=p.X();fY=p.Y();fZ=p.Z();}
- AliHMPIDHit(Int_t c,Float_t e,Int_t id,Int_t tn,Float_t x,Float_t y ): fCh(c ),fPid(id),fQ(0),fLx(x),fLy(y){QdcTot(e);fTrack=tn;}//manual ctor
+ AliHMPIDHit( ):AliHit( ),fCh(-1),fPid(-1 ),fQ(-1),fLx(0),fLy(0){} //default ctor
+ AliHMPIDHit(Int_t c,Float_t &e,Int_t pid,Int_t tid,Float_t x,Float_t y,const TVector3 &p):AliHit(0,tid),fCh(c ),fPid(pid),fQ(0 ),fLx(x),fLy(y){QdcTot(e);fX=p.X();fY=p.Y();fZ=p.Z();}
+ AliHMPIDHit(Int_t c,Float_t &e,Int_t pid,Int_t tid,Float_t x,Float_t y ):AliHit( ),fCh(c ),fPid(pid),fQ(0 ),fLx(x),fLy(y){QdcTot(e);fTrack=tid;}//manual ctor
+ AliHMPIDHit(const AliHMPIDHit &h):AliHit(h),fCh(h.fCh),fPid(h.fPid),fQ(h.fQ),fLx(h.fLx),fLy(h.fLy) {}//copy ctor
virtual ~AliHMPIDHit() {}
//framework part
- void Print(Option_t *opt="")const; //from TObject to print current status
- void Draw (Option_t *opt=""); //from TObject to Draw this hit in current canvas
+ void Print(Option_t *opt="")const; //from TObject to print current status
+ void Draw (Option_t *opt=""); //from TObject to Draw this hit
//private part
- Int_t Ch ( )const{return fCh; } //Chamber
- //clm: hit LorsX and Y were changed for simulation hits
- Float_t LorsX ( )const{return fLx; } //hit X position in LORS, [cm]
- Float_t LorsY ( )const{return fLy; } //hit Y position in LORS, [cm]
- Int_t Pid ( )const{return fPid; } //PID
- Float_t Q ( )const{return fQ; } //Eloss for MIP hit or Etot for photon hit, [GeV]
- inline Float_t QdcTot(Float_t e); //calculate total charge of the hit
- Int_t Tid ( )const{return fTrack; } //TID
+ Int_t Ch ( )const{return fCh; } //Chamber
+ void Hit2Sdi(TClonesArray *pSdiLst,Int_t n=1)const; //add sdigits of this hit to the list
+ Float_t LorsX ( )const{return fLx; } //hit X position in LORS, [cm]
+ Float_t LorsY ( )const{return fLy; } //hit Y position in LORS, [cm]
+ Int_t Pid ( )const{return fPid; } //PID
+ Float_t Q ( )const{return fQ; } //total charge, [QDC]
+ inline Float_t QdcTot (Float_t &e ); //calculate total charge of the hit
+ Int_t Tid ( )const{return fTrack; } //TID
-protected: //AliHit has fTid,fX,fY,fZ
+protected: //AliHit has fTrack,fX,fY,fZ
Int_t fCh; //Chamber
Int_t fPid; //PID
Float_t fQ; //total charge [QDC]
ClassDef(AliHMPIDHit,5) //HMPID hit class
};//class AliHMPIDhit
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Float_t AliHMPIDHit::QdcTot(Float_t e)
+Float_t AliHMPIDHit::QdcTot(Float_t &e)
{
// Samples total charge of the hit
-// Arguments: e- hit energy [GeV]; for mip Eloss for photon Etot
+// Arguments: e- hit energy [GeV] for mip Eloss for photon Etot
// Returns: total QDC
- Float_t x=(LorsX() > 66.6)? LorsX()-66.6:LorsX(); //sagita is for PC (0-64) and not for chamber
+ Int_t pc,px,py;
+ AliHMPIDDigit::Lors2Pad(fLx,fLy,pc,px,py);
+ Float_t y=AliHMPIDDigit::LorsY(pc,py);
+ fLy=((y-fLy)>0)?y-0.2:y+0.2; //shift to the nearest anod wire
+
+ Float_t x=(fLx > 66.6)? fLx-66.6:fLx; //sagita is for PC (0-64) and not for chamber
Float_t qdcEle=34.06311+0.2337070*x+5.807476e-3*x*x-2.956471e-04*x*x*x+2.310001e-06*x*x*x*x; //reparametrised from DiMauro
Int_t iNele=Int_t(e/26e-9); if(iNele<1) iNele=1; //number of electrons created by hit
- for(Int_t i=1;i<=iNele;i++) fQ-=qdcEle*TMath::Log(gRandom->Rndm()+1e-6); //1e-6 is a protection against 0 from rndm
+ fQ=0;
+ for(Int_t i=1;i<=iNele;i++){
+ Double_t rnd=gRandom->Rndm(); if(rnd==0) rnd=1e-12; //1e-12 is a protection against 0 from rndm
+ fQ-=qdcEle*TMath::Log(rnd);
+ }
+ e=fQ;
return fQ;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
+
typedef AliHMPIDHit AliRICHHit; // for backward compatibility
#endif
// Note that TGeoManager should be already initialized from geometry.root file
fX=0.5*AliHMPIDDigit::SizeAllX();
fY=0.5*AliHMPIDDigit::SizeAllY();
- for(Int_t i=0;i<7;i++) fM[i]=(TGeoHMatrix*)gGeoManager->GetVolume("ALIC")->GetNode(Form("HMPID_%i",i))->GetMatrix();
+ for(Int_t i=0;i<7;i++)
+ if(gGeoManager)
+ fM[i]=(TGeoHMatrix*)gGeoManager->GetVolume("ALIC")->GetNode(Form("HMPID_%i",i))->GetMatrix();
+ else{
+ fM[i]=new TGeoHMatrix;
+ IdealPosition(i,fM[i]);
+ }
fgInstance=this;
}//ctor
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
for(Int_t i=0;i<7;i++) fM[i]->Print(opt);
}//Print()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void AliHMPIDParam::IdealPosition(Int_t iCh, TGeoHMatrix *pMatrix)
+{
+// Construct ideal position matrix for a given chamber
+// Arguments: iCh- chamber ID; pMatrix- pointer to precreated unity matrix where to store the results
+// Returns: none
+ const Double_t kAngHor=19.5; // horizontal angle between chambers 19.5 grad
+ const Double_t kAngVer=20; // vertical angle between chambers 20 grad
+ const Double_t kAngCom=30; // common HMPID rotation with respect to x axis 30 grad
+ const Double_t trans[3]={490,0,0}; //center of the chamber is on window-gap surface
+ pMatrix->RotateY(90); //rotate around y since initial position is in XY plane -> now in YZ plane
+ pMatrix->SetTranslation(trans); //now plane in YZ is shifted along x
+ switch(iCh){
+ case 0: pMatrix->RotateY(kAngHor); pMatrix->RotateZ(-kAngVer); break; //right and down
+ case 1: pMatrix->RotateZ(-kAngVer); break; //down
+ case 2: pMatrix->RotateY(kAngHor); break; //right
+ case 3: break; //no rotation
+ case 4: pMatrix->RotateY(-kAngHor); break; //left
+ case 5: pMatrix->RotateZ(kAngVer); break; //up
+ case 6: pMatrix->RotateY(-kAngHor); pMatrix->RotateZ(kAngVer); break; //left and up
+ }
+ pMatrix->RotateZ(kAngCom); //apply common rotation in XY plane
+
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Int_t AliHMPIDParam::Stack(Int_t evt,Int_t tid)
{
// Prints some useful info from stack
Double_t MeanIdxWin () {return 1.57819;}//???????????
static Int_t Stack(Int_t evt=-1,Int_t tid=-1); //Print stack info for event and tid
static Int_t StackCount(Int_t pid,Int_t evt); //Counts stack particles of given sort in given event
-//trasformation methodes
+ static void IdealPosition(Int_t iCh,TGeoHMatrix *m); //ideal position of given chamber
+ //trasformation methodes
void Lors2Mars (Int_t c,Float_t x,Float_t y,Double_t *m,Int_t pl=kPc)const{Double_t z=0; switch(pl){case kPc:z=8.0;break; case kAnod:z=7.806;break; case kRad:z=-1.25; break;} Double_t l[3]={x-fX,y-fY,z}; fM[c]->LocalToMaster(l,m); }
TVector3 Lors2Mars (Int_t c,Float_t x,Float_t y, Int_t pl=kPc)const{Double_t m[3];Lors2Mars(c,x,y,m,pl); return TVector3(m); }//MRS->LRS
void Mars2Lors (Int_t c,Double_t *m,Float_t &x,Float_t &y )const{Double_t l[3];fM[c]->MasterToLocal(m,l);x=l[0]+fX;y=l[1]+fY;}//MRS->LRS
// Return pointer to the AliHMPIDParam singleton.
// Arguments: none
// Returns: pointer to the instance of AliHMPIDParam or 0 if no geometry
- if(!fgInstance)
- if(gGeoManager) new AliHMPIDParam;
- else Printf("AliHMPIDParam> Error:: No geometry defined!");
+ if(!fgInstance) new AliHMPIDParam;
return fgInstance;
}//Instance()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
#endif
// pCluAll - list of clusters for all chambers
// isTryUnfold - flag to choose between CoG and Mathieson fitting
// Returns: none
- TMatrixF padMap(AliHMPIDDigit::kPadAllX,AliHMPIDDigit::kPadAllY); //pads map for single chamber 0..159 x 0..143
+ TMatrixF padMap(AliHMPIDDigit::kMaxPadChX,AliHMPIDDigit::kMaxPadChY); //pads map for single chamber 0..159 x 0..143
for(Int_t iCh=0;iCh<7;iCh++){ //chambers loop
TClonesArray *pDigCur=(TClonesArray*)pDigAll->At(iCh); //get list of digits for current chamber
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
- dig.Raw(ddl,w32);
+ dig.Raw(w32,ddl);
AliDebug(1,Form("Ch=%i DDL=%i raw=0x%x digit=(%3i,%3i,%3i,%3i) Q=%5.2f",iCh,ddl,w32,dig.Ch(),dig.Pc(),dig.PadPcX(),dig.PadPcY(),dig.Q()));
new((*((TClonesArray*)digLst.At(iCh)))[iDigCnt[iCh]++]) AliHMPIDDigit(dig); //add this digit to the tmp list
}//raw records loop
pChain->GetListOfFiles()->Print();
- TVirtualProof *pProof=TProof::Open("kir@lxb6046.cern.ch");
+ TProof *pProof=TProof::Open("kir@lxb6046.cern.ch");
pProof->UploadPackage("ESD.par");
pProof->EnablePackage("ESD");
pDigQ->Fill(pDig->Q());
pDigM->Fill(pDig->LorsX(),pDig->LorsY());
}//digits loop
- Int_t qmax=0;
+ Float_t qmax=0;
for(Int_t iClu=0;iClu<pCluCh->GetEntriesFast();iClu++){//clusters loop
AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluCh->UncheckedAt(iClu);
pCluQ->Fill(pClu->Q());
#include "AliHMPIDv1.h" //class header
-#include "AliHMPIDParam.h" //CreateMaterials()
+#include "AliHMPIDParam.h" //StepManager()
#include "AliHMPIDHit.h" //Hits2SDigs(),StepManager()
-#include "AliHMPIDDigit.h" //CreateMaterials()
+#include "AliHMPIDDigit.h" //Digits2Raw(), Raw2SDigits()
#include "AliRawReader.h" //Raw2SDigits()
-#include <TParticle.h> //Hits2SDigits()
-#include <TRandom.h>
#include <TVirtualMC.h> //StepManager() for gMC
#include <TPDGCode.h> //StepHistory()
#include <AliStack.h> //StepManager(),Hits2SDigits()
#include <AliConst.h>
#include <AliPDG.h>
#include <AliMC.h> //StepManager()
-#include <AliRawDataHeader.h> //Digits2Raw()
-#include <AliDAQ.h> //Digits2Raw()
#include <AliRun.h> //CreateMaterials()
#include <AliMagF.h> //CreateMaterials()
#include <TGeoManager.h> //CreateGeometry()
-#include <TMultiGraph.h> //Optics()
-#include <TGraph.h> //Optics()
-#include <TLegend.h> //Optics()
-#include <TCanvas.h> //Optics()
-#include <TF2.h> //CreateMaterials()
+#include <TF1.h> //DefineOpticalProperties()
+#include <TF2.h> //DefineOpticalProperties()
+#include <TLorentzVector.h> //IsLostByFresnel()
#include <AliCDBManager.h> //CreateMaterials()
#include <AliCDBEntry.h> //CreateMaterials()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void AliHMPIDv1::AddAlignableVolumes()const
{
-// Associates the symbolic volume name with the corresponding volume path. Interface methode from AliModule ivoked from AliMC
+// Associates the symbolic volume name with the corresponding volume path. Interface method from AliModule invoked from AliMC
// Arguments: none
// Returns: none
- for(Int_t i=0;i<7;i++)
+ for(Int_t i=AliHMPIDDigit::kMinCh;i<=AliHMPIDDigit::kMaxCh;i++)
gGeoManager->SetAlignableEntry(Form("/HMPID/Chamber%i",i),Form("ALIC_1/HMPID_%i",i));
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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);
-
DefineOpticalProperties();
-
- AliDebug(1,"Stop v1 HMPID.");
-
- TString ttl=GetTitle(); if(!ttl.Contains("ShowOptics")) return; //user didn't aks to plot optical curves
-
- const Double_t kWidth=0.25,kHeight=0.2;
- const Int_t kRadM=24 , kRadC=kRed;
- const Int_t kWinM=26 , kWinC=kBlue;
- const Int_t kGapM=25 , kGapC=kGreen;
- const Int_t kPcM = 2 , kPcC =kMagenta;
- const Int_t kNbins=30; //number of photon energy points
-
- Float_t aEckov [kNbins];
- Float_t aAbsRad[kNbins], aAbsWin[kNbins], aAbsGap[kNbins];
- Float_t aIdxRad[kNbins], aIdxWin[kNbins], aIdxGap[kNbins];
- Float_t aQePc [kNbins];
- Float_t aTraRad[kNbins],aTraWin[kNbins],aTraGap[kNbins],aTraTot[kNbins];
- for(Int_t i=0;i<kNbins;i++){//calculate probability for photon to survive during transversing a volume of material with absorption length
- aTraRad[i]=TMath::Exp(-AliHMPIDDigit::SizeRad()/ (aAbsRad[i]+0.0001)); //radiator
- aTraWin[i]=TMath::Exp(-AliHMPIDDigit::SizeWin()/ (aAbsWin[i] +0.0001)); //window
- aTraGap[i]=TMath::Exp(-AliHMPIDDigit::SizeGap()/ (aAbsGap[i] +0.0001)); //from window to PC
- aTraTot[i]=aTraRad[i]*aTraWin[i]*aTraGap[i]*aQePc[i];
- }
-
- TGraph *pRaAG=new TGraph(kNbins,aEckov,aAbsRad);pRaAG->SetMarkerStyle(kRadM);pRaAG->SetMarkerColor(kRadC);
- TGraph *pRaIG=new TGraph(kNbins,aEckov,aIdxRad);pRaIG->SetMarkerStyle(kRadM);pRaIG->SetMarkerColor(kRadC);
- TGraph *pRaTG=new TGraph(kNbins,aEckov,aTraRad);pRaTG->SetMarkerStyle(kRadM);pRaTG->SetMarkerColor(kRadC);
-
- TGraph *pWiAG=new TGraph(kNbins,aEckov,aAbsWin);pWiAG->SetMarkerStyle(kWinM);pWiAG->SetMarkerColor(kWinC);
- TGraph *pWiIG=new TGraph(kNbins,aEckov,aIdxWin);pWiIG->SetMarkerStyle(kWinM);pWiIG->SetMarkerColor(kWinC);
- TGraph *pWiTG=new TGraph(kNbins,aEckov,aTraWin);pWiTG->SetMarkerStyle(kWinM);pWiTG->SetMarkerColor(kWinC);
-
- TGraph *pGaAG=new TGraph(kNbins,aEckov,aAbsGap);pGaAG->SetMarkerStyle(kGapM);pGaAG->SetMarkerColor(kGapC);
- TGraph *pGaIG=new TGraph(kNbins,aEckov,aIdxGap);pGaIG->SetMarkerStyle(kGapM);pGaIG->SetMarkerColor(kGapC);
- TGraph *pGaTG=new TGraph(kNbins,aEckov,aTraGap);pGaTG->SetMarkerStyle(kGapM);pGaTG->SetMarkerColor(kGapC);
-
- TGraph *pQeG =new TGraph(kNbins,aEckov,aQePc); pQeG ->SetMarkerStyle(kPcM );pQeG->SetMarkerColor(kPcC);
- TGraph *pToG =new TGraph(kNbins,aEckov,aTraTot);pToG ->SetMarkerStyle(30) ;pToG->SetMarkerColor(kYellow);
-
- TMultiGraph *pIdxMG=new TMultiGraph("idx","Ref index;E_{#check{C}} [GeV]");
- TMultiGraph *pAbsMG=new TMultiGraph("abs","Absorption [cm];E_{#check{C}} [GeV]");
- TMultiGraph *pTraMG=new TMultiGraph("tra","Transmission;E_{#check{C}} [GeV]"); TLegend *pTraLe=new TLegend(0.2,0.4,0.2+kWidth,0.4+kHeight);
- pAbsMG->Add(pRaAG); pIdxMG->Add(pRaIG); pTraMG->Add(pRaTG); pTraLe->AddEntry(pRaTG, "Rad", "p");
- pAbsMG->Add(pWiAG); pIdxMG->Add(pWiIG); pTraMG->Add(pWiTG); pTraLe->AddEntry(pWiTG, "Win", "p");
- pAbsMG->Add(pGaAG); pIdxMG->Add(pGaIG); pTraMG->Add(pGaTG); pTraLe->AddEntry(pGaTG, "Gap", "p");
- pTraMG->Add(pToG); pTraLe->AddEntry(pToG, "Tot", "p");
- pTraMG->Add(pQeG); pTraLe->AddEntry(pQeG, "QE" , "p");
- TCanvas *pC=new TCanvas("c1","HMPID optics to check",1100,900); pC->Divide(2,2);
- pC->cd(1); pIdxMG->Draw("AP");
- pC->cd(2); gPad->SetLogy(); pAbsMG->Draw("AP");
- pC->cd(3); pTraLe->Draw();
- pC->cd(4); pTraMG->Draw("AP");
}//void AliHMPID::CreateMaterials()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void AliHMPIDv1::CreateGeometry()
TGeoVolume *pRich=gGeoManager->MakeBox("HMPID",gGeoManager->GetMedium("HMPID_CH4"),dx=(6*mm+1681*mm+6*mm)/2, //main HMPID volume
dy=(6*mm+1466*mm+6*mm)/2,
dz=(80*mm+40*mm)*2/2); //x,y taken from 2033P1 z from p84 TDR
- const Double_t kAngHor=19.5; // horizontal angle between chambers 19.5 grad
- const Double_t kAngVer=20; // vertical angle between chambers 20 grad
- const Double_t kAngCom=30; // common HMPID rotation with respect to x axis 30 grad
- const Double_t trans[3]={490,0,0}; //center of the chamber is on window-gap surface
- for(Int_t iCh=0;iCh<7;iCh++){//place 7 chambers
+ for(Int_t iCh=AliHMPIDDigit::kMinCh;iCh<=AliHMPIDDigit::kMaxCh;iCh++){//place 7 chambers
TGeoHMatrix *pMatrix=new TGeoHMatrix;
- pMatrix->RotateY(90); //rotate around y since initial position is in XY plane -> now in YZ plane
- pMatrix->SetTranslation(trans); //now plane in YZ is shifted along x
- switch(iCh){
- case 0: pMatrix->RotateY(kAngHor); pMatrix->RotateZ(-kAngVer); break; //right and down
- case 1: pMatrix->RotateZ(-kAngVer); break; //down
- case 2: pMatrix->RotateY(kAngHor); break; //right
- case 3: break; //no rotation
- case 4: pMatrix->RotateY(-kAngHor); break; //left
- case 5: pMatrix->RotateZ(kAngVer); break; //up
- case 6: pMatrix->RotateY(-kAngHor); pMatrix->RotateZ(kAngVer); break; //left and up
- }
- pMatrix->RotateZ(kAngCom); //apply common rotation in XY plane
+ AliHMPIDParam::IdealPosition(iCh,pMatrix);
gGeoManager->GetVolume("ALIC")->AddNode(pRich,iCh,pMatrix);
}
return kFALSE;
}//IsLostByFresnel()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPIDv1::GenFee(Int_t iCh,Float_t eloss)
+void AliHMPIDv1::GenFee(Float_t qtot)
{
// Generate FeedBack photons for the current particle. To be invoked from StepManager().
-// eloss=0 means photon so only pulse height distribution is to be analysed. This one is done in AliHMPIDParam::TotQdc()
+// eloss=0 means photon so only pulse height distribution is to be analysed.
TLorentzVector x4;
gMC->TrackPosition(x4);
- Int_t iNphotons=gMC->GetRandom()->Poisson(0.02*200); eloss++; iCh++; //??????????????????????
+ Int_t iNphotons=gMC->GetRandom()->Poisson(0.02*qtot); //# of feedback photons is proportional to the charge of hit
AliDebug(1,Form("N photons=%i",iNphotons));
Int_t j;
Float_t cthf, phif, enfp = 0, sthf, e1[3], e2[3], e3[3], vmod, uswop,dir[3], phi,pol[3], mom[4];
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void AliHMPIDv1::Hits2SDigits()
{
-// Interface methode ivoked from AliSimulation to create a list of sdigits corresponding to list of hits. Every hit generates one or more sdigits.
+// Interface method ivoked from AliSimulation to create a list of sdigits corresponding to list of hits. Every hit generates one or more sdigits.
// Arguments: none
// Returns: none
AliDebug(1,"Start.");
if(!GetLoader()->TreeH()) {GetLoader()->LoadHits(); }
if(!GetLoader()->TreeS()) {GetLoader()->MakeTree("S"); MakeBranch("S");}//to
- for(Int_t iPrimN=0;iPrimN<GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop
- GetLoader()->TreeH()->GetEntry(iPrimN);
+ for(Int_t iEnt=0;iEnt<GetLoader()->TreeH()->GetEntries();iEnt++){//prims loop
+ GetLoader()->TreeH()->GetEntry(iEnt);
Hit2Sdi(Hits(),SdiLst());
}//prims loop
GetLoader()->TreeS()->Fill();
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void AliHMPIDv1::Hit2Sdi(TClonesArray *pHitLst,TClonesArray *pSdiLst)
{
-// Converts list of hits to list of sdigits. For each hit in a loop the following steps are done:
-// - calcultion of the total charge induced by the hit
-// - determination of the pad contaning the hit and shifting hit y position to the nearest anod wire y
-// - defining a set of pads affected (up to 9 including the hitted pad)
-// - calculating charge induced to all those pads using integrated Mathieson distribution and creating sdigit
+// Converts list of hits to list of sdigits.
// Arguments: pHitLst - list of hits provided not empty
// pSDigLst - list of sdigits where to store the results
// Returns: none
for(Int_t iHit=0;iHit<pHitLst->GetEntries();iHit++){ //hits loop
- AliHMPIDHit *pHit=(AliHMPIDHit*)pHitLst->At(iHit); //get pointer to current hit
- AliHMPIDDigit::Hit2Sdi(pHit,pSdiLst); //convert this hit to list of sdigits
+ AliHMPIDHit *pHit=(AliHMPIDHit*)pHitLst->At(iHit); //get pointer to current hit
+ pHit->Hit2Sdi(pSdiLst); //convert this hit to list of sdigits
}//hits loop loop
-}//Hits2SDigs() for TVector2
+}//Hits2Sdi()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void AliHMPIDv1::Digits2Raw()
{
-// Creates raw data files in DDL format. Invoked by AliSimulation where loop over events is done
+// Interface method invoked by AliSimulation to create raw data streams from digits. Events loop is done in AliSimulation
// Arguments: none
// Returns: none
AliDebug(1,"Start.");
}
treeD->GetEntry(0);
- ofstream file[AliHMPIDDigit::kNddls]; //output streams
- Int_t cnt[AliHMPIDDigit::kNddls]; //data words counters for DDLs
- AliRawDataHeader header; //empty DDL header
- UInt_t w32=0; //32 bits data word
-
- for(Int_t i=0;i<AliHMPIDDigit::kNddls;i++){
- file[i].open(AliDAQ::DdlFileName(GetName(),i)); //open all 14 DDL in parallel
- file[i].write((char*)&header,sizeof(header)); //write dummy header as place holder, actual will be written later when total size of DDL is known
- cnt[i]=0; //reset counters
- }
-
- for(Int_t iCh=0;iCh<7;iCh++)
- for(Int_t iDig=0;iDig<DigLst(iCh)->GetEntriesFast();iDig++){//digits loop for a given chamber
- AliHMPIDDigit *pDig=(AliHMPIDDigit*)DigLst(iCh)->At(iDig);
- Int_t ddl=pDig->Raw(w32); //ddl is 0..13
- file[ddl].write((char*)&w32,sizeof(w32)); cnt[ddl]++;//write formated digit to the propriate file (as decided in Dig2Raw) and increment corresponding counter
- }//digits
+ AliHMPIDDigit::WriteRaw(DigLst());
-
- for(Int_t i=0;i<AliHMPIDDigit::kNddls;i++){
- header.fSize=sizeof(header)+cnt[i]*sizeof(w32); //now calculate total number of bytes for each DDL file
- header.SetAttribute(0);
- file[i].seekp(0); file[i].write((char*)&header,sizeof(header));//rewrite DDL header with fSize field properly set
- file[i].close(); //close DDL file
- }
GetLoader()->UnloadDigits();
AliDebug(1,"Stop.");
}//Digits2Raw()
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 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 PC
- GenFee(copy); //generate feedback photons
+ 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
+ GenFee(etot); //generate feedback photons etot is modified in hit ctor to Q of hit
}//photon hit PC and DE >0
}//photon hit PC
//Treat charged particles
- static Double_t dEdX; //need to store mip parameters between different steps
+ static Float_t eloss; //need to store mip parameters between different steps
static Double_t in[3];
if(gMC->TrackCharge() && gMC->CurrentVolID(copy)==fIdAmpGap){ //charged particle in amplification gap (fIdAmpGap)
if(gMC->IsTrackEntering()||gMC->IsNewTrack()) { //entering or newly created
- dEdX=0; //reset dEdX collector
+ eloss=0; //reset Eloss collector
gMC->TrackPosition(in[0],in[1],in[2]); //take position at the entrance
}else if(gMC->IsTrackExiting()||gMC->IsTrackStop()||gMC->IsTrackDisappeared()){ //exiting or disappeared
- dEdX +=gMC->Edep(); //take into account last step dEdX
+ eloss +=gMC->Edep(); //take into account last step Eloss
gMC->CurrentVolOffID(1,copy); //take current chamber since geometry tree is HMPID-Rgap
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
out[0]=0.5*(out[0]+in[0]); out[1]=0.5*(out[1]+in[1]); out[1]=0.5*(out[1]+in[1]); //take hit position at the anod plane
- Float_t xl,yl;AliHMPIDParam::Instance()->Mars2Lors(copy,out,xl,yl); //take LORS position
- new((*fHits)[fNhits++])AliHMPIDHit(copy,dEdX,pid,tid,xl,yl,out); //HIT for MIP, position near anod plane
- GenFee(copy,dEdX); //generate feedback photons
+ 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
+ GenFee(eloss); //generate feedback photons
}else //just going inside
- dEdX += gMC->Edep(); //collect this step dEdX
+ eloss += gMC->Edep(); //collect this step eloss
}//MIP in GAP
}//StepManager()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Bool_t Raw2SDigits (AliRawReader * ); //from AliMOdule invoked from AliSimulation
void StepManager ( ); //from AliModule invoked from AliMC::Stepping()
//private part++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
- void GenFee (Int_t iCh,Float_t eloss=0 ); //generates feedback photons; eloss=0 for photon
+ void GenFee (Float_t qtot ); //generates feedback photons
static Float_t Fresnel (Float_t geV,Float_t p, Bool_t pl ); //deals with Fresnel absorption on PC
static void Hit2Sdi (TClonesArray *pH,TClonesArray *pS);
Bool_t IsLostByFresnel ( ); //checks if the photon lost on Fresnel reflection
#pragma link off all classes;
#pragma link off all functions;
-#pragma link C++ typedef AliRICHHelix;
+
#pragma link C++ typedef AliRICHDigit;
#pragma link C++ typedef AliRICHHit;
#pragma link C++ typedef AliRICHCluster;
#pragma link C++ typedef AliRICHParam;
#pragma link C++ typedef AliRICH;
-#pragma link C++ class AliHMPIDHelix+;
+
#pragma link C++ class AliHMPIDDigit+;
#pragma link C++ class AliHMPIDHit+;
#pragma link C++ class AliHMPIDCluster+;
AliHMPID *h; AliLoader *hl; AliHMPIDParam *hp;
Bool_t isGeomType=kFALSE;
-Int_t nEvent=0;
-
+Int_t nCurEvt=0;
+Int_t nMaxEvt=0;
+TControlBar *pMenu=0;
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void GetParam()
{
{
TString status="Status: ";
if(gSystem->IsFileInIncludePath("galice.root")){
- status+="galice.root found";
+ status+="galice.root: ";
al=AliRunLoader::Open(); //try to open galice.root from current dir
if(gAlice) delete gAlice; //in case we execute this in aliroot delete default AliRun object
al->LoadgAlice(); a=al->GetAliRun(); //take new AliRun object from galice.root
hl=al->GetDetectorLoader("HMPID"); h=(AliHMPID*)a->GetDetector("HMPID"); //get HMPID object from galice.root
- status+=Form(" with %i event(s)",al->GetNumberOfEvents()); status+=(h)? " with HMPID": " without HMPID";
+ status+=(h)? "HMPID": "PROBLEM PROBLEM PROBLEM- no HMPID";
+ nMaxEvt=al->GetNumberOfEvents()-1;
+ status+=Form(" Event(s) 0-%i",nMaxEvt);
}else
- status+="No galice.root";
+ status+="PROBLEM PROBLEM PROBLEM no galice.root";
+ status+=Form(" curent event %i",nCurEvt);
GetParam();
- TControlBar *pMenu = new TControlBar("horizontal",status.Data(),0,0);
+ pMenu = new TControlBar("horizontal",status.Data(),0,0);
pMenu->AddButton(" ","","");
pMenu->AddButton(" General ","General()" ,"general items which do not depend on any files");
pMenu->AddButton(" ","" ,"");
void General()
{
TControlBar *pMenu = new TControlBar("vertical","General purpose",100,50);
- pMenu->AddButton("Debug ON","don();" ,"Switch debug on-off" );
- pMenu->AddButton("Debug OFF","doff();" ,"Switch debug on-off" );
+ pMenu->AddButton("Debug ON","don();" ,"Switch debug on-off" );
+ pMenu->AddButton("Debug OFF","doff();" ,"Switch debug on-off" );
pMenu->AddButton("Geo GUI","geo();" ,"Shows geometry" );
pMenu->AddButton("Browser","new TBrowser;" ,"Start ROOT TBrowser" );
pMenu->Show();
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void SimData()
{
- TControlBar *pMenu = new TControlBar("vertical","Sim data",310,50);
- pMenu->AddButton("Display ","hed();" ,"Display Fast");
- pMenu->AddButton("HITS QA" ,"hqa()" ,"QA plots for hits: hqa()");
- pMenu->AddButton("Print stack" ,"stack();" ,"To print hits: hp(evt)");
- pMenu->AddButton("Print hits" ,"hp();" ,"To print hits: hp(evt)");
- pMenu->AddButton("Print sdigits" ,"sp();" ,"To print sdigits: sp(evt)");
- pMenu->AddButton("Print digits" ,"dp();" ,"To print digits: dp(evt)");
- pMenu->AddButton("Print clusters" ,"cp();" ,"To print clusters: cp(evt)");
- pMenu->Show();
+ TControlBar *pSim = new TControlBar("vertical","Sim data",310,50);
+ pSim->AddButton("Display ","hed();" ,"Display Fast");
+ pSim->AddButton("HITS QA" ,"hqa()" ,"QA plots for hits: hqa()");
+ pSim->AddButton("Print stack" ,"stack();" ,"To print hits: hp(evt)");
+ pSim->AddButton("Print hits" ,"hp(nCurEvt);" ,"To print hits: hp(evt)");
+ pSim->AddButton("Print sdigits" ,"sp(nCurEvt);" ,"To print sdigits: sp(evt)");
+ pSim->AddButton("Print digits" ,"dp(nCurEvt);" ,"To print digits: dp(evt)");
+ pSim->AddButton("Print clusters" ,"cp(nCurEvt);" ,"To print clusters: cp(evt)");
+ pSim->Show();
}//SimData()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void RawData()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void Test()
{
- TControlBar *pMenu = new TControlBar("vertical","Test",820,50);
- pMenu->AddButton("TEST Display " ,"sed();" ,"Display Fast");
- pMenu->AddButton("Hits->Digits" ,"thd();" ,"test hits->sdigits->digits" );
- pMenu->AddButton("Segmentation" ,"ts()" ,"test segmentation methods" );
- pMenu->AddButton("Test response" ,"AliHMPIDParam::TestResp();","Test AliHMPIDParam response methods" );
- pMenu->AddButton("Print map" ,"PrintMap();" ,"Test AliHMPIDParam transformation methods" );
- pMenu->AddButton("Test Recon" ,"rec();" ,"Test AliHMPIDRecon" );
- pMenu->Show();
+ TControlBar *pTst = new TControlBar("vertical","Test",625,50);
+ pTst->AddButton("TEST Display " ,"sed();" ,"Display Fast");
+ pTst->AddButton("Test all" ,"tst();" ,"test hits->sdigits->digits" );
+ pTst->AddButton("Segmentation" ,"ts()" ,"test segmentation methods" );
+ pTst->AddButton("Test response" ,"AliHMPIDParam::TestResp();","Test AliHMPIDParam response methods" );
+ pTst->AddButton("Print map" ,"PrintMap();" ,"Test AliHMPIDParam transformation methods" );
+ pTst->AddButton("Test Recon" ,"rec();" ,"Test AliHMPIDRecon" );
+ pTst->Show();
}//Test()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void du ( ) {h->Dump ( );} //utility display
-void hp ( ) {h->HitPrint (nEvent);} //print hits for requested event
-void hq ( ) {h->HitQA ( );} //hits QA plots for all events
-void sp ( ) {h->SdiPrint (nEvent);} //print sdigits for requested event
-void sq ( ) {h->SdiPrint (nEvent);} //print sdigits for requested event
-void dp ( ) {h->DigPrint (nEvent);} //print digits for requested event
-void dq ( ) {AliHMPIDReconstructor::DigQA (al );} //digits QA plots for all events
-void cp ( ) {h->CluPrint (nEvent); } //print clusters for requested event
-void cq ( ) {AliHMPIDReconstructor::CluQA (al );} //clusters QA plots for all events
-
-void ep ( ) {AliHMPIDTracker::EsdQA(1); }
-void eq ( ) {AliHMPIDTracker::EsdQA(); }
-void mp (Double_t probCut=0.7 ) {AliHMPIDTracker::MatrixPrint(probCut);}
-
-void PrevEvent() {nEvent--;if(nEvent<0)nEvent=0;Printf("Current event is %i",nEvent);}
-void NextEvent() {nEvent++;Printf("Current event is %i",nEvent);}
+void PrevEvent() {nCurEvt--;if(nCurEvt<0 )nCurEvt=0 ;pMenu->SetTitle(Form("Event(s): 0-%i Current event %i",nMaxEvt,nCurEvt));}
+void NextEvent() {nCurEvt++;if(nCurEvt>=nMaxEvt)nCurEvt=nMaxEvt;pMenu->SetTitle(Form("Event(s): 0-%i Current event %i",nMaxEvt,nCurEvt));}
void stack( ) {AliHMPIDParam::Stack();}
void tid (Int_t tid,Int_t evt=0) {AliHMPIDParam::Stack(evt,tid);}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void PrintMap()
{
Printf("\n\n\n");
for(int ch=6;ch>=0;ch--){
- AliHMPIDDigit dL(ch,0,1,0,0),dR(ch,0,1,67,0);
+ AliHMPIDDigit dL,dR; dL.Manual2(ch,2,0 ,24);
+ dR.Manual2(ch,3,79,24);
TVector3 lt=rp->Lors2Mars(ch,0,y); TVector3 rt=rp->Lors2Mars(ch,x,y);
TVector3 ce=rp->Lors2Mars(ch,x/2,y/2);
TVector3 lb=rp->Lors2Mars(ch,0,0); TVector3 rb=rp->Lors2Mars(ch,x,0);
Printf("(ch=%i,pc=%i,x=%2i,y=%2i) ddl=%i raw=0x%h (r=%2i,d=%2i,a=%2i)",
ch, pc, px, py, ddl, raw, r, d, a); dd.Print();
ddl=7;raw=0x592e000;r=22;d=4;a=46;ch=3;pc=1;
-
-
}//PrintMap()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void HitQA(Double_t cut,Double_t cutele,Double_t cutR)
+void DrawSeg()
{
-// Provides a set of control plots intended primarily for charged particle flux analisys
-// Arguments: cut (GeV) - cut on momentum of any charged particles but electrons,
-// cetele (GeV) - the same for electrons-positrons
-// cutR (cm) - cut on production vertex radius (cylindrical system)
- gBenchmark->Start("HitsAna");
-
- Double_t cutPantiproton =cut;
- Double_t cutPkaonminus =cut;
- Double_t cutPpionminus =cut;
- Double_t cutPmuonminus =cut;
- Double_t cutPpositron =cutele;
-
- Double_t cutPelectron =cutele;
- Double_t cutPmuonplus =cut;
- Double_t cutPpionplus =cut;
- Double_t cutPkaonplus =cut;
- Double_t cutPproton =cut;
-
-
- TH2F *pEleHitRZ =new TH2F("EleHitRZ" ,Form("e^{+} e^{-} hit %s;z[cm];R[cm]" ,GetName()) , 400,-300,300 ,400,-500,500); //R-z
- TH2F *pEleHitRP =new TH2F("EleHitRP" ,Form("e^{+} e^{-} hit %s;p[GeV];R[cm]",GetName()) ,1000,-1 ,1 ,400, 0,550); //R-p
- TH1F *pEleAllP =new TH1F("EleAllP" , "e^{+} e^{-} all;p[GeV]" ,1000,-1 ,1 );
- TH1F *pEleHitP =new TH1F("EleHitP" ,Form("e^{+} e^{-} hit %s;p[GeV]" ,GetName()) ,1000,-1 ,1 );
- TH1F *pMuoHitP =new TH1F("MuoHitP" ,Form("#mu^{-} #mu^{+} hit %s;p[GeV]" ,GetName()) ,1000,-4 ,4 );
- TH1F *pPioHitP =new TH1F("PioHitP" ,Form("#pi^{-} #pi^{+} hit %s;p[GeV]" ,GetName()) ,1000,-4 ,4 );
- TH1F *pKaoHitP =new TH1F("KaoHitP" ,Form("K^{-} K^{+} hit %s;p[GeV]" ,GetName()) ,1000,-4 ,4 );
- TH1F *pProHitP =new TH1F("ProHitP" ,Form("p^{-} p^{+} hit %s;p[GeV]" ,GetName()) ,1000,-4 ,4 );
- TH2F *pFlux =new TH2F("flux" ,Form("%s flux with Rvertex<%.1fcm" ,GetName(),cutR),10 ,-5 ,5 , 10,0 ,10); //special text hist
- TH2F *pVertex =new TH2F("vertex" ,Form("%s 2D vertex of HMPID hit;x;y" ,GetName()) ,120 ,0 ,600 ,120,0 ,600); //special text hist
- TH1F *pRho =new TH1F("rho" ,Form("%s r of HMPID hit" ,GetName()) ,600 ,0 ,600); //special text hist
- pFlux->SetStats(0);
- pFlux->GetXaxis()->SetBinLabel(1 ,Form("p^{-}>%.3fGeV/c" ,cutPantiproton));
- pFlux->GetXaxis()->SetBinLabel(2 ,Form("K^{-}>%.3fGeV/c" ,cutPkaonminus ));
- pFlux->GetXaxis()->SetBinLabel(3 ,Form("#pi^{-}>%.3fGeV/c" ,cutPpionminus ));
- pFlux->GetXaxis()->SetBinLabel(4 ,Form("#mu^{-}>%.3fGeV/c" ,cutPmuonminus ));
- pFlux->GetXaxis()->SetBinLabel(5 ,Form("e^{+}>%.3fGeV/c" ,cutPpositron ));
-
- pFlux->GetXaxis()->SetBinLabel(6 ,Form("e^{-}>%.3fGeV/c" ,cutPelectron ));
- pFlux->GetXaxis()->SetBinLabel(7 ,Form("#mu^{+}>%.3fGeV/c" ,cutPmuonplus ));
- pFlux->GetXaxis()->SetBinLabel(8 ,Form("#pi^{+}>%.3fGeV/c" ,cutPpionplus ));
- pFlux->GetXaxis()->SetBinLabel(9 ,Form("K^{+}>%.3fGeV/c" ,cutPkaonplus ));
- pFlux->GetXaxis()->SetBinLabel(10,Form("p^{+}>%.3fGeV/c" ,cutPproton ));
-
- pFlux->GetYaxis()->SetBinLabel(1,"sum");
- pFlux->GetYaxis()->SetBinLabel(2,"ch1");
- pFlux->GetYaxis()->SetBinLabel(3,"ch2");
- pFlux->GetYaxis()->SetBinLabel(4,"ch3");
- pFlux->GetYaxis()->SetBinLabel(5,"ch4");
- pFlux->GetYaxis()->SetBinLabel(6,"ch5");
- pFlux->GetYaxis()->SetBinLabel(7,"ch6");
- pFlux->GetYaxis()->SetBinLabel(8,"ch7");
- pFlux->GetYaxis()->SetBinLabel(9,"prim");
- pFlux->GetYaxis()->SetBinLabel(10,"tot");
-
-//end of hists definition
-
- Int_t iNevents=fLoader->GetRunLoader()->GetAliRun()->GetEventsPerRun(),iCntPrimParts=0,iCntTotParts=0;
-//load all needed trees
- fLoader->LoadHits();
- fLoader->GetRunLoader()->LoadHeader();
- fLoader->GetRunLoader()->LoadKinematics();
-
- for(Int_t iEvtN=0;iEvtN < iNevents;iEvtN++){//events loop
- fLoader->GetRunLoader()->GetEvent(iEvtN);
- AliInfo(Form(" %i event processes",fLoader->GetRunLoader()->GetEventNumber()));
- AliStack *pStack= fLoader->GetRunLoader()->Stack();
-
- for(Int_t iParticleN=0;iParticleN<pStack->GetNtrack();iParticleN++){//stack loop
- TParticle *pPart=pStack->Particle(iParticleN);
-
- if(iParticleN%10000==0) AliInfo(Form(" %i particles read",iParticleN));
-
- switch(pPart->GetPdgCode()){
- case kProtonBar: pFlux->Fill(-4.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-4.5,8); break;
- case kKMinus: pFlux->Fill(-3.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-3.5,8); break;
- case kPiMinus: pFlux->Fill(-2.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-2.5,8); break;
- case kMuonMinus: pFlux->Fill(-1.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-1.5,8); break;
- case kPositron: pFlux->Fill(-0.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-0.5,8); pEleAllP->Fill(-pPart->P()); break;
-
- case kElectron: pFlux->Fill( 0.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 0.5,8); pEleAllP->Fill( pPart->P()); break;
- case kMuonPlus: pFlux->Fill( 1.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 1.5,8); break;
- case kPiPlus: pFlux->Fill( 2.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 2.5,8); break;
- case kKPlus: pFlux->Fill( 3.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 3.5,8); break;
- case kProton: pFlux->Fill( 4.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 4.5,8); break;
- }//switch
- }//stack loop
-//now hits analiser
- for(Int_t iEntryN=0;iEntryN < fLoader->TreeH()->GetEntries();iEntryN++){//TreeH loop
- fLoader->TreeH()->GetEntry(iEntryN); //get current entry (prim)
- for(Int_t iHitN=0;iHitN < Hits()->GetEntries();iHitN++){//hits loop
- AliHMPIDHit *pHit = (AliHMPIDHit*)Hits()->At(iHitN); //get current hit
- TParticle *pPart=pStack->Particle(pHit->GetTrack()); //get stack particle which produced the current hit
-
- if(pPart->GetPDG()->Charge()!=0&&pPart->Rho()>0.1) pVertex->Fill(pPart->Vx(),pPart->Vy()); //safe margin for sec.
- if(pPart->GetPDG()->Charge()!=0) pRho->Fill(pPart->Rho()); //safe margin for sec.
- if(pPart->R()>cutR) continue; //cut on production radius (cylindrical system)
-
- switch(pPart->GetPdgCode()){
- case kProtonBar: if(pPart->P()>cutPantiproton) {pProHitP->Fill(-pPart->P()); pFlux->Fill(-4.5,pHit->Ch());}break;
- case kKMinus : if(pPart->P()>cutPkaonminus) {pKaoHitP->Fill(-pPart->P()); pFlux->Fill(-3.5,pHit->Ch());}break;
- case kPiMinus : if(pPart->P()>cutPpionminus) {pPioHitP->Fill(-pPart->P()); pFlux->Fill(-2.5,pHit->Ch());}break;
- case kMuonMinus: if(pPart->P()>cutPmuonminus) {pMuoHitP->Fill(-pPart->P()); pFlux->Fill(-1.5,pHit->Ch());}break;
- case kPositron : if(pPart->P()>cutPpositron) {pEleHitP->Fill(-pPart->P()); pFlux->Fill(-0.5,pHit->Ch());
- pEleHitRP->Fill(-pPart->P(),pPart->R()); pEleHitRZ->Fill(pPart->Vz(),pPart->R()); }break;
-
- case kElectron : if(pPart->P()>cutPelectron) {pEleHitP->Fill( pPart->P()); pFlux->Fill( 0.5,pHit->Ch());
- pEleHitRP->Fill( pPart->P(),pPart->R()); pEleHitRZ->Fill(pPart->Vz(),pPart->R()); }break;
- case kMuonPlus : if(pPart->P()>cutPmuonplus) {pMuoHitP->Fill( pPart->P()); pFlux->Fill( 1.5,pHit->Ch());}break;
- case kPiPlus : if(pPart->P()>cutPpionplus) {pPioHitP->Fill( pPart->P()); pFlux->Fill( 2.5,pHit->Ch());}break;
- case kKPlus : if(pPart->P()>cutPkaonplus) {pKaoHitP->Fill( pPart->P()); pFlux->Fill( 3.5,pHit->Ch());}break;
- case kProton : if(pPart->P()>cutPproton) {pProHitP->Fill( pPart->P()); pFlux->Fill( 4.5,pHit->Ch());}break;
- }
- }//hits loop
- }//TreeH loop
- iCntPrimParts +=pStack->GetNprimary();
- iCntTotParts +=pStack->GetNtrack();
- }//events loop
-//unload all loaded staff
- fLoader->UnloadHits();
- fLoader->GetRunLoader()->UnloadHeader();
- fLoader->GetRunLoader()->UnloadKinematics();
-//Calculater some sums
- Stat_t sum=0;
-//sum row, sum over rows
- for(Int_t i=1;i<=pFlux->GetNbinsX();i++){
- sum=0; for(Int_t j=2;j<=8;j++) sum+=pFlux->GetBinContent(i,j);
- pFlux->SetBinContent(i,1,sum);
- }
-
-//display everything
- new TCanvas("canvas1",Form("Events %i Nprims=%i Nparticles=%i",iNevents,iCntPrimParts,iCntTotParts),1000,900); pFlux->Draw("text"); gPad->SetGrid();
-//total prims and particles
- TLatex latex; latex.SetTextSize(0.02);
- sum=0; for(Int_t i=1;i<=pFlux->GetNbinsX();i++) sum+=pFlux->GetBinContent(i,10); latex.DrawLatex(5.1,9.5,Form("%.0f",sum));
- sum=0; for(Int_t i=1;i<=pFlux->GetNbinsX();i++) sum+=pFlux->GetBinContent(i,9); latex.DrawLatex(5.1,8.5,Form("%.0f",sum));
- for(Int_t iCh=0;iCh<7;iCh++) {
- sum=0; for(Int_t i=1;i<=pFlux->GetNbinsX();i++) sum+=pFlux->GetBinContent(i,iCh+2);latex.DrawLatex(5.1,iCh+1.5,Form("%.0f",sum));
- }
- sum=0; for(Int_t i=1;i<=pFlux->GetNbinsX();i++) sum+=pFlux->GetBinContent(i,1); latex.DrawLatex(5.1,0.5,Form("%.0f",sum));
-
- new TCanvas("cEleAllP" ,"e" ,200,100); pEleAllP->Draw();
- new TCanvas("cEleHitRP" ,"e" ,200,100); pEleHitRP->Draw();
- new TCanvas("cEleHitRZ" ,"e" ,200,100); pEleHitRZ->Draw();
- new TCanvas("cEleHitP" ,"e" ,200,100); pEleHitP->Draw();
- new TCanvas("cMuoHitP" ,"mu",200,100); pMuoHitP->Draw();
- new TCanvas("cPioHitP" ,"pi",200,100); pPioHitP->Draw();
- new TCanvas("cKaoHitP" ,"K" ,200,100); pKaoHitP->Draw();
- new TCanvas("cProHitP" ,"p" ,200,100); pProHitP->Draw();
- new TCanvas("cVertex" ,"2d vertex" ,200,100); pVertex->Draw();
- new TCanvas("cRho" ,"Rho of sec" ,200,100); pRho->Draw();
-
- gBenchmark->Show("HitsPlots");
-}//HitQA()
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void CluQA(AliRunLoader *pAL)
-{
-// Quality assesment plots for clusters.
-// This methode takes list of digits and form list of clusters again in order to
-// calculate cluster shape and cluster particle mixture
- AliLoader *pRL=pAL->GetDetectorLoader("HMPID"); AliHMPID *pRich=(AliHMPID*)pAL->GetAliRun()->GetDetector("HMPID");//get pointers for HMPID and HMPID loader
- Int_t iNevt=pAL->GetNumberOfEvents(); if(iNevt==0) {AliInfoClass("No events");return;}
- if(pRL->LoadDigits()) {AliInfoClass("No digits file");return;}
- pAL->LoadHeader();
- pAL->LoadKinematics();
-// AliStack *pStack=pAL->Stack();
- TH1::AddDirectory(kFALSE);
-
-
- TH1F* pQ=new TH1F("HmpAllQ" ,"Charge All" ,4000 ,0 ,4000);// Q hists
- TH1F* pCerQ=new TH1F("HmpCerQ" ,"Charge Ckov" ,4000 ,0 ,4000);
- TH1F* pMipQ=new TH1F("HmpMipQ" ,"Charge MIP" ,4000 ,0 ,4000);
-
- TH1F* pS=new TH1F("HmpCluSize" ,"Cluster size;size" ,100 ,0 ,100 );// size hists
- TH1F* pCerS=new TH1F("HmpCluCerSize" ,"Ckov size;size" ,100 ,0 ,100 );
- TH1F* pMipS=new TH1F("HmpCluMipSize" ,"MIP size;size" ,100 ,0 ,100 );
-
- TH2F* pM=new TH2F("HmpCluMap" ,"Cluster map;x [cm];y [cm]" ,1000 ,0 ,AliHMPIDDigit::SizeAllX(),1000,0,AliHMPIDDigit::SizeAllY()); // maps
- TH2F* pMipM=new TH2F("HmpCluMipMap" ,"MIP map;x [cm];y [cm]" ,1000 ,0 ,AliHMPIDDigit::SizeAllX(),1000,0,AliHMPIDDigit::SizeAllY());
- TH2F* pCerM=new TH2F("HmpCluCerMap" ,"Ckov map;x [cm];y [cm]" ,1000 ,0 ,AliHMPIDDigit::SizeAllX(),1000,0,AliHMPIDDigit::SizeAllY());
-
-
-
- for(Int_t iEvt=0;iEvt<iNevt;iEvt++){
- pAL->GetEvent(iEvt);
- pRL->TreeD()->GetEntry(0);
- TClonesArray *pCluLst=new TClonesArray("AliHMPIDCluster");//tmp list of clusters for this event
-
- for(Int_t iCh=0;iCh<7;iCh++) AliHMPIDReconstructor::Dig2Clu(pRich->DigLst(iCh),pCluLst,kFALSE);//cluster finder for all chamber if any digits present
-
- for(Int_t iClu=0;iClu<pCluLst->GetEntriesFast();iClu++){
- AliHMPIDCluster *pClu = (AliHMPIDCluster*)pCluLst->At(iClu);
- Int_t cfm=0; for(Int_t iDig=0;iDig<pClu->Size();iDig++) cfm+=pClu->Dig(iDig)->Ch(); //collect ckov-fee-mip structure of current cluster ?????
- Int_t iNckov=cfm/1000000; Int_t iNfee =cfm%1000000/1000; Int_t iNmip =cfm%1000000%1000;
-
- pQ ->Fill(pClu->Q()) ; pS ->Fill(pClu->Size()) ; pM ->Fill(pClu->X(),pClu->Y()); //all clusters
- if(iNckov!=0 && iNfee==0 && iNmip==0) {pCerQ->Fill(pClu->Q()) ; pCerS->Fill(pClu->Size()) ; pCerM ->Fill(pClu->X(),pClu->Y());}//ckov only cluster
- if(iNckov==0 && iNfee==0 && iNmip!=0) {pMipQ->Fill(pClu->Q()) ; pMipS->Fill(pClu->Size()) ; pMipM ->Fill(pClu->X(),pClu->Y());}//mip only cluster
-
- }//clusters loop
- pCluLst->Clear();delete pCluLst;
- }//events loop
-
- pRL->UnloadDigits(); pAL->UnloadKinematics(); pAL->UnloadHeader();
- TCanvas *pC=new TCanvas("RichCluQA",Form("QA for cluster from %i events",iNevt),1000,900); pC->Divide(3,3);
- pC->cd(1); pM->Draw(); pC->cd(2); pQ->Draw(); pC->cd(3); pS->Draw();
- pC->cd(4); pMipM->Draw(); pC->cd(5); pMipQ->Draw(); pC->cd(6); pMipS->Draw();
- pC->cd(7); pCerM->Draw(); pC->cd(8); pCerQ->Draw(); pC->cd(9); pCerS->Draw();
-}//CluQA()
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void OccupancyPrint(Int_t iEvtNreq)
-{
-//prints occupancy for each chamber in a given event
- Int_t iEvtNmin,iEvtNmax;
- if(iEvtNreq==-1){
- iEvtNmin=0;
- iEvtNmax=gAlice->GetEventsPerRun();
- } else {
- iEvtNmin=iEvtNreq;iEvtNmax=iEvtNreq+1;
- }
-
- if(GetLoader()->GetRunLoader()->LoadHeader()) return;
- if(GetLoader()->GetRunLoader()->LoadKinematics()) return;
-
-// Info("Occupancy","for event %i",iEvtN);
- if(GetLoader()->LoadHits()) return;
- if(GetLoader()->LoadDigits()) return;
-
-
- for(Int_t iEvtN=iEvtNmin;iEvtN<iEvtNmax;iEvtN++){
- Int_t nDigCh[7]={0,0,0,0,0,0,0};
- Int_t iChHits[7]={0,0,0,0,0,0,0};
- Int_t nPrim[7]={0,0,0,0,0,0,0};
- Int_t nSec[7]={0,0,0,0,0,0,0};
- AliInfo(Form("events processed %i",iEvtN));
- if(GetLoader()->GetRunLoader()->GetEvent(iEvtN)) return;
- AliStack *pStack = GetLoader()->GetRunLoader()->Stack();
- for(Int_t iPrimN=0;iPrimN<GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop
- GetLoader()->TreeH()->GetEntry(iPrimN);
- for(Int_t iHitN=0;iHitN<Hits()->GetEntries();iHitN++){
- AliHMPIDHit *pHit = (AliHMPIDHit*)Hits()->At(iHitN);
- if(pHit->E()>0){
- iChHits[pHit->Ch()]++;
- if(pStack->Particle(pHit->GetTrack())->Rho()<0.01) nPrim[pHit->Ch()]++;else nSec[pHit->Ch()]++;
- }
- }
- }
-
- GetLoader()->TreeD()->GetEntry(0);
- for(Int_t iCh=0;iCh<7;iCh++){
- for(Int_t iDig=0;iDig<DigLst(iCh)->GetEntries();iDig++){
- AliHMPIDDigit *pDig=(AliHMPIDDigit*)DigLst(iCh)->At(iDig);
- nDigCh[pDig->Ch()]++;
- }
- Printf("Occupancy for chamber %i = %4.2f %% and charged prim tracks %i and sec. tracks %i with total %i",
- iCh,Float_t(nDigCh[iCh])*100/AliHMPIDDigit::kPadAll,nPrim[iCh],nSec[iCh],iChHits[iCh]);
- }
-
-
- }//events loop
- GetLoader()->UnloadHits();
- GetLoader()->UnloadDigits();
- GetLoader()->GetRunLoader()->UnloadHeader();
- GetLoader()->GetRunLoader()->UnloadKinematics();
-}//OccupancyPrint()
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void AliHMPID::SummaryOfEvent(Int_t iEvtN) const
-{
-//prints a summary for a given event
- AliInfo(Form("Summary of event %i",iEvtN));
- GetLoader()->GetRunLoader()->GetEvent(iEvtN);
- if(GetLoader()->GetRunLoader()->LoadHeader()) return;
- if(GetLoader()->GetRunLoader()->LoadKinematics()) return;
- AliStack *pStack=GetLoader()->GetRunLoader()->Stack();
-
- AliGenEventHeader* pGenHeader = gAlice->GetHeader()->GenEventHeader();
- if(pGenHeader->InheritsFrom("AliGenHijingEventHeader")) {
- AliInfo(Form(" Hijing event with impact parameter b = %.2f (fm)",((AliGenHijingEventHeader*) pGenHeader)->ImpactParameter()));
- }
- Int_t nChargedPrimaries=0;
- for(Int_t i=0;i<pStack->GetNtrack();i++) {
- TParticle *pParticle = pStack->Particle(i);
- if(pParticle->IsPrimary()&&pParticle->GetPDG()->Charge()!=0) nChargedPrimaries++;
- }
- AliInfo(Form("Total number of primaries %i",pStack->GetNprimary()));
- AliInfo(Form("Total number of charged primaries %i",nChargedPrimaries));
- AliInfo(Form("Total n. of tracks in stack(+sec) %i",pStack->GetNtrack()));
- GetLoader()->GetRunLoader()->UnloadHeader();
- GetLoader()->GetRunLoader()->UnloadKinematics();
-}//SummaryOfEvent()
+ TCanvas *pC=new TCanvas("seg","Segmentation as seen from electronics side");
+ DrawPc(1);
+ new TGedEditor(pC);
+ pC->ToggleEventStatus();
+ pC->SetEditable(0);
+ pC->AddExec("status","DrawStatus()");
+}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void DrawZoom()
-{
-// Show info about current cursur position in status bar of the canvas
-// Arguments: none
-// Returns: none
- TCanvas *pC=(TCanvas*)gPad;
+void DrawStatus()
+{// Show info about current cursur position in status bar of the canvas
+// Printf("event %i",gPad->GetEvent()); return;
+ TPad *pad=gPad;
+ TCanvas *pC=(TCanvas*)pad;
TRootCanvas *pRC= (TRootCanvas*)pC->GetCanvasImp();
TGStatusBar *pBar=pRC->GetStatusBar();
pBar->SetParts(5);
- Float_t x=gPad->AbsPixeltoX(gPad->GetEventX());
- Float_t y=gPad->AbsPixeltoY(gPad->GetEventY());
- AliHMPIDDigit dig;dig.Manual1(1,x,y); UInt_t w32=0;
- if(IsInDead(x,y))
+ Float_t x=pad->AbsPixeltoX(pad->GetEventX()); Float_t y=pad->AbsPixeltoY(pad->GetEventY());
+ if(AliHMPIDDigit::IsInDead(x,y))
pBar->SetText("Out of sensitive area",4);
else{
- Int_t ddl=dig.Raw(w32);
- pBar->SetText(Form("(p%i,x%i,y%i) ddl=%i 0x%x (r%i,d%i,a%i) (%.2f,%.2f)",
- dig.Pc(),dig.PadPcX(),dig.PadPcY(),
- ddl,w32,
- dig.Row(),dig.Dilogic(),dig.Addr(),
- dig.LorsX(),dig.LorsY() ),4);
+ Int_t pc,px,py,w32,ddl,r,d,a; AliHMPIDDigit::Lors2Pad(x,y,pc,px,py); AliHMPIDDigit dig; dig.Set(1,pc,px,py); dig.Raw(w32,ddl,r,d,a);
+ pBar->SetText(Form("(pc%i,px%i,py%i) (r%i,d%i,a%i) (%.2f,%.2f)",
+ pc ,px ,py, r ,d ,a ,dig.LorsX(),dig.LorsY()),4);
}
- if(gPad->GetEvent()==1){
- new TCanvas("zoom",Form("Row %i DILOGIC %i",dig.Row(),dig.Dilogic()));
- }
-}//DrawZoom()
+// if(pad->GetEvent()==1){
+// new TCanvas("zoom",Form("Row %i DILOGIC %i",dig.Row(),dig.Dilogic()));
+// }
+}//DrawStatus()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void DrawPc(Bool_t isFill)
{
-// Utility methode draws HMPID chamber PCs on event display.
-// Arguments: none
-// Returns: none
-// y6 ---------- ----------
-// | | | |
-// | 4 | | 5 |
-// y5 ---------- ----------
-//
-// y4 ---------- ----------
-// | | | |
-// | 2 | | 3 | view from electronics side
-// y3 ---------- ----------
-//
-// y2 ---------- ----------
-// | | | |
-// | 0 | | 1 |
-// y1 ---------- ----------
-// x1 x2 x3 x4
gPad->Range(-10,-10,AliHMPIDDigit::SizeAllX()+5,AliHMPIDDigit::SizeAllY()+5);
- Float_t x1=0,
- x2=AliHMPIDDigit::SizePcX(),
- x3=AliHMPIDDigit::SizePcX()+AliHMPIDDigit::SizeDead(),
- x4=AliHMPIDDigit::SizeAllX();
- Float_t y1=0,
- y2= AliHMPIDDigit::SizePcY(),
- y3= AliHMPIDDigit::SizePcY()+AliHMPIDDigit::SizeDead(),
- y4=2*AliHMPIDDigit::SizePcY()+AliHMPIDDigit::SizeDead(),
- y5= AliHMPIDDigit::SizeAllY()-AliHMPIDDigit::SizePcY(),
- y6= AliHMPIDDigit::SizeAllY();
-
- Float_t xL[5]={x1,x1,x2,x2,x1}; //clockwise
- Float_t xR[5]={x3,x3,x4,x4,x3};
- Float_t yD[5]={y1,y2,y2,y1,y1};
- Float_t yC[5]={y3,y4,y4,y3,y3};
- Float_t yU[5]={y5,y6,y6,y5,y5};
+
Float_t dX2=0.5*AliHMPIDDigit::SizePadX(),
dY2=0.5*AliHMPIDDigit::SizePadY() ;
TLatex txt; txt.SetTextSize(0.01);
- Int_t iColLeft=29,iColRight=41;
- TPolyLine *pc=0; TLine *pL;
- AliHMPIDDigit dig;
- for(Int_t iPc=0;iPc<AliHMPIDDigit::kPcAll;iPc++){
- if(iPc==4) pc=new TPolyLine(5,xL,yU); if(iPc==5) pc=new TPolyLine(5,xR,yU); //draw PCs
- if(iPc==2) pc=new TPolyLine(5,xL,yC); if(iPc==3) pc=new TPolyLine(5,xR,yC);
- if(iPc==0) pc=new TPolyLine(5,xL,yD); if(iPc==1) pc=new TPolyLine(5,xR,yD);
- (iPc%2)? pc->SetFillColor(iColLeft): pc->SetFillColor(iColRight);
- if(!isFill){ pc->Draw();continue;}
+ TLine *pL;
+
+ AliHMPIDDigit dig; UInt_t w32; Int_t ddl,r,d,a;
+
+ for(Int_t iPc=AliHMPIDDigit::kMinPc;iPc<=AliHMPIDDigit::kMaxPc;iPc++){
+ TBox *pBox=new TBox(AliHMPIDDigit::fMinPcX[iPc],AliHMPIDDigit::fMinPcY[iPc],
+ AliHMPIDDigit::fMaxPcX[iPc],AliHMPIDDigit::fMaxPcY[iPc]);
+
+ if(iPc==0||iPc==2||iPc==4) pBox->SetFillColor(29);
+ else pBox->SetFillColor(41);
+ pBox->Draw();
- pc->Draw("f");
- if(iPc%2) {dig.Manual2(0,iPc,79,25); txt.DrawText(dig.LorsX()+2,dig.LorsY(),Form("PC%i",dig.Pc()));}//print PC#
+ if(!isFill) continue;
+// if(iPc%2) {dig.Set(0,iPc,79,25); txt.DrawText(dig.LorsX()+2,dig.LorsY(),Form("PC%i",dig.Pc()));}//print PC#
+
txt.SetTextAlign(32);
for(Int_t iRow=0;iRow<8 ;iRow++){//draw row lines (horizontal)
- dig.Manual2(0,iPc,0,iRow*6); //set digit to the left-down pad of this row
+ dig.Set(0,iPc,0,iRow*6); dig.Raw(w32,ddl,r,d,a); //set digit to the left-down pad of this row
+
if(iPc%2) txt.DrawText(dig.LorsX()-1 ,dig.LorsY(),Form("%i",dig.PadPcY())); //print PadY#
- txt.DrawText(dig.LorsX()-1+(iPc%2)*67,dig.LorsY()+2,Form("r%i",dig.Row())); //print Row#
+ txt.DrawText(dig.LorsX()-1+(iPc%2)*67,dig.LorsY()+2,Form("r%i",r)); //print Row#
pL=new TLine(dig.LorsX()-dX2,dig.LorsY()-dY2,dig.LorsX()+AliHMPIDDigit::SizePcX()-dX2,dig.LorsY()-dY2);//draw horizontal line
if(iRow!=0) pL->Draw();
}//row loop
txt.SetTextAlign(13);
for(Int_t iDil=0;iDil<10;iDil++){//draw dilogic lines (vertical)
- dig.Manual2(0,iPc,iDil*8,0); //set this digit to the left-down pad of this dilogic
+ dig.Set(0,iPc,iDil*8,0); dig.Raw(w32,ddl,r,d,a); //set this digit to the left-down pad of this dilogic
+
txt.DrawText(dig.LorsX() ,dig.LorsY()-1,Form("%i",dig.PadPcX())); //print PadX#
- if(iPc==4 || iPc==5) txt.DrawText(dig.LorsX()+2,dig.LorsY()+42,Form("d%i",dig.Dilogic())); //print Dilogic#
+ if(iPc==4 || iPc==5) txt.DrawText(dig.LorsX()+2,dig.LorsY()+42,Form("d%i",d)); //print Dilogic#
pL=new TLine(dig.LorsX()-dX2,dig.LorsY()-dY2,dig.LorsX()-dX2,dig.LorsY()+AliHMPIDDigit::SizePcY()-dY2); //draw vertical line
if(iDil!=0)pL->Draw();
}//dilogic loop
static TFile *pEsdFl=0;
static TTree *pEsdTr=0;
static AliESD *pEsd=0;
+
+
+
+
if(!pC&&iEvt<iEvtTot){
iEvt=0;
iEvtTot=999;
pEsdFl=TFile::Open("AliESDs.root"); if(!pEsdFl || !pEsdFl->IsOpen()) return;//open AliESDs.root
pEsdTr=(TTree*) pEsdFl->Get("esdTree"); if(!pEsdTr) return;//get ESD tree
pEsdTr->SetBranchAddress("ESD", &pEsd);
- hl->LoadDigits(); hl->LoadRecPoints();
+ hl->LoadHits(); hl->LoadDigits(); hl->LoadRecPoints();
iEvtTot=pEsdTr->GetEntries();
pC=new TCanvas("hed","View from electronics side, IP is behind the picture.",1000,900); pC->ToggleEventStatus(); pC->Divide(3,3);
pC->cd(7); TButton *pBtn=new TButton("Next","hed()",0,0,0.2,0.1); pBtn->Draw();
pC->cd(7); TButton *pBtn=new TButton("Quit","Close_hed()",0.2,0,0.4,0.1); pBtn->Draw();
+ new TGedEditor(pC);
}
+ TLatex txt; txt.SetTextSize(0.1);
+ TClonesArray hits("AliHMPIDHit");
+
if(iEvt<iEvtTot){
- pEsdTr->GetEntry(iEvt); al->GetEvent(iEvt); hl->TreeD()->GetEntry(0); hl->TreeR()->GetEntry(0);
- TLatex txt; pC->cd(3); gPad->Clear();
- txt.SetTextSize(0.1);txt.DrawLatex(0.2,0.2,Form("Event %i (total %i)",iEvt,iEvtTot));
- DrawEvt(pC,h->DigLst(),h->CluLst(),pEsd);
+ pEsdTr->GetEntry(iEvt); al->GetEvent(iEvt);
+ hl->TreeD()->GetEntry(0); hl->TreeR()->GetEntry(0);
+ ReadHits(&hits);
+
+ pC->cd(3); gPad->Clear(); txt.DrawLatex(0.2,0.2,Form("Event %i (total %i)",iEvt,iEvtTot));
+ DrawEvt(pC,&hits,h->DigLst(),h->CluLst(),pEsd);
+
iEvt++;
}else{
Printf("--- No more events available...Bye.");
pC=0x0;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void ReadHits(TClonesArray *pHitLst)
+{
+ pHitLst->Delete(); Int_t cnt=0;
+ for(Int_t iEnt=0;iEnt<hl->TreeH()->GetEntries();iEnt++){ //TreeH loop
+ hl->TreeH()->GetEntry(iEnt); //get current entry (prim)
+ for(Int_t iHit=0;iHit<h->Hits()->GetEntries();iHit++){ //hits loop
+ AliHMPIDHit *pHit = (AliHMPIDHit*)h->Hits()->At(iHit); //get current hit
+ new((*pHitLst)[cnt++]) AliHMPIDHit(*pHit);
+ }//hits loop for this entry
+ }//tree entries loop
+
+}//ReadHits()
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void sed()
{
AliESD esd;
- EsdFromStack(&esd);
+ SimEsd(&esd);
HitsFromEsd(&esd,&lh);
AliHMPIDReconstructor::Dig2Clu(&ld,&lc);
// AliHMPIDTracker::Recon(&esd,&cl);
- DrawEvt(pC1,&ld,&lc,&esd);
+ DrawEvt(pC1,&lh,&ld,&lc,&esd);
}//SimEvt()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-void DrawEvt(TCanvas *pC,TObjArray *pDigLst,TObjArray *pCluLst,AliESD *pEsd)
+void DrawEvt(TCanvas *pC,TClonesArray *pHitLst,TObjArray *pDigLst,TObjArray *pCluLst,AliESD *pEsd)
{//draws all the objects of current event
AliHMPIDRecon rec;
}
gPad->SetEditable(kTRUE); gPad->Clear();
DrawPc(0);
+ for(Int_t iHit=0;iHit<pHitLst->GetEntries();iHit++){
+ AliHMPIDHit *pHit=(AliHMPIDHit*)pHitLst->At(iHit);
+ if(pHit->Ch()==iCh) pHit->Draw(); //draw hits
+ }
((TClonesArray*)pDigLst->At(iCh))->Draw(); //draw digits
((TClonesArray*)pCluLst->At(iCh))->Draw(); //draw clusters
pTxC[iCh]->Draw(); //draw intersections
pRin[iCh]->Draw(); //draw rings
gPad->SetEditable(kFALSE);
+// gPad->AddExec("zoom","DrawZoom()");
}//chambers loop
+
+
// TLatex txt; txt.SetTextSize(0.02);
// txt.DrawLatex(20,-5,Form("#theta=%.4f#pm%.5f with %2i #check{C}" ,simCkov,simErr,simN));
// txt.DrawLatex(25,-5,Form("#theta=%.4f#pm%.5f with %2i #check{C}" ,recCkov,recErr,recN));
// Printf("DIG------DIG---------DIG--------DIG------DIG------DIG");pDigLst->Print();Printf("");
}//DrawEvt()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void DrawZoom()
+{
+ TPad *pad=gPad; Float_t x=gPad->AbsPixeltoX(pad->GetEventX()); Float_t y=gPad->AbsPixeltoY(pad->GetEventY());
+ TCanvas *zoom = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("zoom");
+ if(!zoom) zoom=new TCanvas("zoom","");
+ zoom->SetTitle(Form("Zoom view around %.2f %.2f",x,y));
+ gPad->Range(x-20,y-20,x+20,y+20);
+ zoom->Update();
+}
void t1(Int_t case=1)
{
AliHMPIDDigit *d[10]; for(Int_t i=0;i<10;i++) d[i]=new AliHMPIDDigit;
cl->Print();
}//t1()
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void SimEsd(AliESD *pEsd)
+{
+ TParticle part; TLorentzVector mom;
+ for(Int_t iTrk=0;iTrk<25;iTrk++){//stack loop
+ part.SetPdgCode(kProton);
+ part.SetProductionVertex(0,0,0,0);
+ Double_t eta= -0.4+gRandom->Rndm()*0.8; //rapidity is random [-0.4,+0.4]
+ Double_t phi= gRandom->Rndm()*1.4; //phi is random [ 0 , 80 ] degrees
+ mom.SetPtEtaPhiM(2,eta,phi,part.GetMass());
+ part.SetMomentum(mom);
+ AliESDtrack trk(&part);
+ pEsd->AddTrack(&trk);
+ }//stack loop
+}//EsdFromStack()
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void EsdFromStack(AliESD *pEsd)
{
al->LoadHeader();al->LoadKinematics();
if(ch<0) continue; //this track does not hit HMPID
Float_t ckov=0.63;
- Float_t th,ph,xPc,yPc,; pTrk->GetHMPIDtrk(xPc,yPc,th,ph);
- new((*pHitLst)[hc++]) AliHMPIDHit(ch,200e-9,kProton ,iTrk,xPc ,yPc); //mip hit
- for(int i=0;i<4;i++) new((*pHitLst)[hc++]) AliHMPIDHit(ch,7.5e-9,kFeedback,iTrk,gRandom->Rndm()*130,gRandom->Rndm()*126); //bkg hits 4 per track
- for(int i=0;i<16;i++){
- rec.TracePhot(ckov,gRandom->Rndm()*TMath::TwoPi(),pos);
- new((*pHitLst)[hc++]) AliHMPIDHit(ch,7.5e-9,kCerenkov,iTrk,pos.X(),pos.Y());} //photon hits
+ Float_t th,ph,xPc,yPc,; pTrk->GetHMPIDtrk(xPc,yPc,th,ph); rec.SetTrack(xRa,yRa,th,ph);
+
+ if(!AliHMPIDDigit::IsInDead(xPc,yPc)) new((*pHitLst)[hc++]) AliHMPIDHit(ch,200e-9,kProton ,iTrk,xPc,yPc); //mip hit
+// for(int i=0;i<4;i++) new((*pHitLst)[hc++]) AliHMPIDHit(ch,7.5e-9,kFeedback,iTrk,gRandom->Rndm()*130,gRandom->Rndm()*126); //bkg hits 4 per track
+// for(int i=0;i<16;i++){
+// rec.TracePhot(ckov,gRandom->Rndm()*TMath::TwoPi(),pos);
+// new((*pHitLst)[hc++]) AliHMPIDHit(ch,7.5e-9,kCerenkov,iTrk,pos.X(),pos.Y());} //photon hits
}//tracks loop
+}//HitsFromEsd()
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void tst(Int_t nEvts=10000,Int_t type=999)
+{
+
+ TLegend *lQ=new TLegend(0.5,0.5,0.9,0.9);
+
+ TH1F *hQ7 =new TH1F("hQ7" ,"" ,300,-50,2000); hQ7 ->SetLineColor(kRed); lQ->AddEntry(hQ7 ,"Ckov 7 eV"); hQ7->SetStats(0);
+ TH1F *hQ200=new TH1F("hQ200","" ,300,-50,2000); hQ200->SetLineColor(kBlack); lQ->AddEntry(hQ200,"mip 200 eV");
+ TH1F *hQ500=new TH1F("hQ500","" ,300,-50,2000); hQ500->SetLineColor(kCyan); lQ->AddEntry(hQ500,"mip 500 eV");
+ TH1F *hQ900=new TH1F("hQ900","" ,300,-50,2000); hQ900->SetLineColor(kGreen); lQ->AddEntry(hQ900,"mip 900 eV");
+
+ TH1F *hCluPerEvt=new TH1F("hCluPerEvt","# clusters per event",11,-0.5,10.5);
+ TH1F *hCluChi2 =new TH1F("hChi2","Chi2 ",1000,0,100);
+ TH1F *hCluFlg =new TH1F("hCluFlg","Cluster flag",14,-1.5,12.5); hCluFlg->SetFillColor(5);
+ TH1F *hCluRawSize= new TH1F("hCluRawSize","Raw cluster size ",100,0,100);
+
+ TH2F *pCluMapSi1=new TH2F("cluMSi1","Size 1 map",1700,-10,160,1700,-10,160);
+ TH2F *pCluMapLo1=new TH2F("cluMLo1","Loc Max 1 map",1700,-10,160,1700,-10,160);
+ TH2F *pCluMapEmp=new TH2F("cluMEmp","Should be empty",1700,-10,160,1700,-10,160);
+ TH2F *pCluMapUnf=new TH2F("cluMUnf","Unfolded map",1700,-10,160,1700,-10,160);
+ TH2F *pCluMapMax=new TH2F("cluMMax","Max Loc Max map",1700,-10,160,1700,-10,160);
+ TH2F *pCluMapAbn=new TH2F("cluMAbn","Fit failed map",1700,-10,160,1700,-10,160);
+ TH2F *pCluMapEdg=new TH2F("cluMEdg","On edge map",1700,-10,160,1700,-10,160);
+ TH2F *pCluMapNoLoc=new TH2F("cluMNoLoc","No Loc map",1700,-10,160,1700,-10,160);
+ TH2F *pCluMapCoG=new TH2F("cluMCoG","CoG map",1700,-10,160,1700,-10,160);
+
+ TH1F *hHitCluDifX = new TH1F("hHitCluDifX" ,";entries;x_{Hit}-x_{Clu} [cm]" ,2000,-2,2); hHitCluDifX->Sumw2(); hHitCluDifX->SetFillColor(kYellow);
+ TH1F *hHitCluDifY = new TH1F("hHitCluDifY" ,";entries;y_{Hit}-y_{Clu} [cm]" ,2000,-2,2); hHitCluDifY->Sumw2(); hHitCluDifY->SetFillColor(kYellow);
+ TH2F *hHitCluDifXY= new TH2F("hHitCluDifXY",";x_{Hit}-x_{Clu};y_{Hit}-y_{Clu}",2000,-2,2,2000,-2,2);hHitCluDifXY->Sumw2();
+ TH1F *hHitCluDifQ = new TH1F("hHitCluDifQ" ,";entries;Q_{Hit}-Q_{Clu}" ,200 ,-100,100); hHitCluDifQ->Sumw2(); hHitCluDifQ->SetFillColor(kYellow);
+
+ Float_t e200=200e-9,e500=500e-9,e900=900e-9,e7=7e-9;//predefined Eloss
+
+ AliHMPIDHit hit(0,0,kProton,0,0,0);
+ TClonesArray hits("AliHMPIDHit"); TClonesArray sdigs("AliHMPIDDigit");
+ TObjArray digs(7); for(Int_t i=0;i<7;i++) digs.AddAt(new TClonesArray("AliHMPIDDigit"),i);
+ TObjArray clus(7); for(Int_t i=0;i<7;i++) clus.AddAt(new TClonesArray("AliHMPIDCluster"),i);
+
+
+ for(Int_t iEvt=0;iEvt<nEvts;iEvt++){//events loop
+ Printf("============> iEvt = %d ",iEvt);
+ Int_t nHits=(type==999)?1:40;
+ Int_t ch,pid; Float_t e,hitx,hity,hitq;
+ for(Int_t iHit=0;iHit<nHits;iHit++){//hits loop for the current event
+ Printf("In loop");
+ switch(iHit){
+ case 0: ch=0;pid=kProton;e=e200;hitx=16.0+gRandom->Rndm()*0.8;hity= 16.8+gRandom->Rndm()*0.84;break; //mip ramdomly distributed in one pad in the middle
+ case 1: ch=0;pid=kProton;e=e200;hitx=0.4;hity=0.42;break; //mip in left-hand bottom coner of chamber 0
+ case 2: ch=0;pid=kProton;e=e200;hitx=0.4;hity=30 ;break; //mip on left edge of chamber 0
+ case 3: ch=0;pid=kProton;e=e200;hitx=40; hity=0.42;break; //mip on bottom edge of chamber 0
+ default: ch=gRandom->Rndm()*6; pid=(gRandom->Rndm()>0.9)? kProton:kCerenkov;
+ if(pid==kProton)
+ e=gRandom->Rndm()*900e-9;
+ else
+ e=5.5e-9+3e-9*gRandom->Rndm();
+ x=gRandom->Rndm()*AliHMPIDDigit::SizeAllX(); y=gRandom->Rndm()*AliHMPIDDigit::SizeAllY();break; //random hit
+ }
+ new(hits[iHit]) AliHMPIDHit(ch,e,pid,iHit,hitx,hity);
+ hitq=e;
+ hQ200->Fill(hit.QdcTot(e200));
+ hQ500->Fill(hit.QdcTot(e500));
+ hQ900->Fill(hit.QdcTot(e900));
+ hQ7 ->Fill(hit.QdcTot(e7));
+ }//hits loop
+
+ hits.Print();
+
+ AliHMPIDv1::Hit2Sdi(&hits,&sdigs);
+ AliHMPIDDigitizer::Sdi2Dig(&sdigs,&digs);
+ AliHMPIDReconstructor::Dig2Clu(&digs,&clus);
+
+ for(Int_t iCh=AliHMPIDDigit::kMinCh;iCh<=AliHMPIDDigit::kMaxCh;iCh++){//chambers loop
+ TClonesArray *pDigs=(TClonesArray *)digs.UncheckedAt(iCh);
+ TClonesArray *pClus=(TClonesArray *)clus.UncheckedAt(iCh);
+
+ hCluPerEvt->Fill(pClus->GetEntriesFast());
+ for(Int_t iClu=0;iClu<pClus->GetEntriesFast();iClu++){//clusters loop
+ AliHMPIDCluster *pClu=(AliHMPIDCluster*)pClus->UncheckedAt(iClu);
+ Float_t clux=pClu->X(); Float_t cluy=pClu->Y(); Float_t cluq=pClu->Q();
+ hCluFlg->Fill(pClu->Status());
+ hCluChi2->Fill(pClu->Chi2());
+ hCluRawSize->Fill(pClu->Size());
+
+ switch(pClu->Status()){
+ case AliHMPIDCluster::kSi1: pCluMapSi1->Fill(clux,cluy); break;
+ case AliHMPIDCluster::kLo1: pCluMapLo1->Fill(clux,cluy); break;
+ case AliHMPIDCluster::kUnf: pCluMapUnf->Fill(clux,cluy); break;
+ case AliHMPIDCluster::kAbn: pCluMapAbn->Fill(clux,cluy); break;
+ case AliHMPIDCluster::kMax: pCluMapMax->Fill(clux,cluy); break;
+ case AliHMPIDCluster::kEdg: pCluMapEdg->Fill(clux,cluy); break;
+ case AliHMPIDCluster::kCoG: pCluMapCoG->Fill(clux,cluy); break;
+ case AliHMPIDCluster::kNoLoc: pCluMapNoLoc->Fill(clux,cluy); break;
+ default: pCluMapEmp->Fill(clux,cluy); break;
+ }
+
+ hHitCluDifX->Fill(hitx-clux); hHitCluDifY->Fill(hity-cluy); hHitCluDifXY->Fill(hitx-clux,hity-cluy); hHitCluDifQ->Fill(hitq-cluq);
+
+ }//clusters loop
+
+ }//chambers loop
+
+ hits.Delete(); sdigs.Delete(); for(int i = 0;i<7;i++){((TClonesArray*)digs.At(i))->Delete();((TClonesArray*)clus.At(i))->Delete();}
+ }//events loop
+
+ gStyle->SetPalette(1);
+ TCanvas *pC2=new TCanvas("Digit canvas","Digit canvas",1280,800); pC2->Divide(3,3);
+ pC2->cd(1);gPad->SetLogy(1);hHitCluDifX->Draw("hist");
+ pC2->cd(2);gPad->SetLogy(1);hHitCluDifY->Draw("hist");
+ pC2->cd(3);gPad->SetLogz(1);hHitCluDifXY->Draw("colz");
+ pC2->cd(4);gPad->SetLogy(1);hHitCluDifQ->Draw("hist");
+ pC2->cd(5);gPad->SetLogy(1);hCluFlg->Draw();
+ pC2->cd(6);gPad->SetLogy(1);hCluChi2->Draw();
+ pC2->cd(7); hCluRawSize->Draw();
+ pC2->cd(8); hCluPerEvt->Draw("colz");
+ pC2->cd(9); hQckov->Draw(); hQm200->Draw("same"); hQm500->Draw("same"); hQm900->Draw("same"); lQ->Draw();
+ TCanvas *pC1=new TCanvas("ClusterMaps","Cluster maps",1280,800); pC1->Divide(3,3);
+ pC1->cd(1); pCluMapSi1->Draw(); DrawPc(kFALSE);
+ pC1->cd(2); pCluMapLo1->Draw(); DrawPc(kFALSE);
+ pC1->cd(3); pCluMapUnf->Draw(); DrawPc(kFALSE);
+ pC1->cd(4); pCluMapAbn->Draw(); DrawPc(kFALSE);
+ pC1->cd(5); pCluMapMax->Draw(); DrawPc(kFALSE);
+ pC1->cd(6); pCluMapEdg->Draw(); DrawPc(kFALSE);
+ pC1->cd(7); pCluMapCoG->Draw(); DrawPc(kFALSE);
+ pC1->cd(8); pCluMapNoLoc->Draw();DrawPc(kFALSE);
+ pC1->cd(9); pCluMapEmp->Draw(); DrawPc(kFALSE);
+
+ pC1->SaveAs("$HOME/HitMaps.png"); //?????
+ pC2->SaveAs("$HOME/HitCluDif.gif");
+
+ Printf("Digits - raw -digits conversion...");
+
+ AliHMPIDDigit d1,d2; Int_t ddl,r,d,a;UInt_t w32;
+ for(Int_t iCh=AliHMPIDDigit::kMinCh;iCh<=AliHMPIDDigit::kMaxCh;iCh++)
+ for(Int_t iPc=AliHMPIDDigit::kMinPc;iPc<=AliHMPIDDigit::kMaxPc;iPc++)
+ for(Int_t iPx=AliHMPIDDigit::kMinPx;iPx<=AliHMPIDDigit::kMaxPx;iPx++)
+ for(Int_t iPy=AliHMPIDDigit::kMinPy;iPy<=AliHMPIDDigit::kMaxPy;iPy++){
+ d1.Set(iCh,iPc,iPx,iPy,3040); //set digit
+ d1.Raw(w32,ddl,r,d,a); //get raw word for this digit
+ d2.Raw(w32,ddl); //set another digit from that raw word
+ if(d1.Compare(&d2)) {d1.Print(); d2.Print(); Printf("");}//compare them
+ }
+ Printf("OK");
+}//tst()
+
+
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void hp(Int_t iEvt=0)
+{
+//Prints a list of HMPID hits for a given event. Default is event number 0.
+ Printf("List of HMPID hits for event %i",iEvt);
+ if(al->GetEvent(iEvt)) return;
+ if(hl->LoadHits()) return;
+
+ Int_t iTotHits=0;
+ for(Int_t iPrim=0;iPrim<hl->TreeH()->GetEntries();iPrim++){//prims loop
+ hl->TreeH()->GetEntry(iPrim);
+ h->Hits()->Print();
+ iTotHits+=h->Hits()->GetEntries();
+ }
+ hl->UnloadHits();
+ Printf("totally %i hits for event %i",iTotHits,iEvt);
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void sp(Int_t iEvt=0)
+{
+//prints a list of HMPID sdigits for a given event
+ Printf("List of HMPID sdigits for event %i",iEvt);
+ if(al->GetEvent(iEvt)) return;
+ if(hl->LoadSDigits()) return;
+
+ hl->TreeS()->GetEntry(0);
+ h->SdiLst()->Print();
+ hl->UnloadSDigits();
+ Printf("totally %i sdigits for event %i",h->SdiLst()->GetEntries(),iEvt);
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void dp(Int_t iEvt=0)
+{
+//prints a list of HMPID digits for a given event
+ Printf("List of HMPID digits for event %i",iEvt);
+ if(al->GetEvent(iEvt)) return;
+ if(hl->LoadDigits()) return;
+
+ hl->TreeD()->GetEntry(0);
+ h->DigLst()->Print();
+ Int_t totDigs=0;
+ for(Int_t i=0;i<7;i++) {totDigs+=h->DigLst(i)->GetEntries();}
+ hl->UnloadDigits();
+ Printf("totally %i digits for event %i",totDigs,iEvt);
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void cp(Int_t iEvt=0)
+{//prints a list of HMPID clusters for a given event
+ Printf("List of HMPID clusters for event %i",iEvt);
+ if(al->GetEvent(iEvt)) return;
+ if(hl->LoadRecPoints()) return;
+
+ hl->TreeR()->GetEntry(0);
+ h->CluLst()->Print();
+
+ Int_t iCluCnt=0; for(Int_t iCh=0;iCh<7;iCh++) iCluCnt+=h->CluLst(iCh)->GetEntries();
+
+ hl->UnloadRecPoints();
+ Printf("totally %i clusters for event %i",iCluCnt,iEvt);
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+void Gui()
+{
+ TGMainFrame *pMF =new TGMainFrame(gClient->GetRoot(),300,400);//main frame
+//1 level widgets: button and 2 horizontal frames
+ TGedFrame *pGedF;
+ pMF->AddFrame(pGedF=new TGedFrame(pMF),new TGLayoutHints(kLHintsExpandY));
+ pMF->AddFrame(pDis1=new TGEmbeddedCanvas(pMF,kSunkenFrame));
+
+ pMF->Layout();
+ pMF->MapSubwindows();
+ pMF->Resize(pMF->GetDefaultSize());
+ pMF->MapWindow();
+}
-SRCS:= AliHMPIDHelix.cxx AliHMPIDHit.cxx AliHMPIDDigit.cxx AliHMPIDCluster.cxx AliHMPIDParam.cxx AliHMPID.cxx
+SRCS:= AliHMPIDHit.cxx AliHMPIDDigit.cxx AliHMPIDCluster.cxx AliHMPIDParam.cxx AliHMPID.cxx
HDRS:= $(SRCS:.cxx=.h)
DHDR:= HMPIDbaseLinkDef.h
git://git.uio.no / u/mrichter/AliRoot.git / commitdiff