From 1d4857c5309d183422e07523c36b66bcdda36c92 Mon Sep 17 00:00:00 2001 From: kir Date: Wed, 28 Feb 2007 15:57:19 +0000 Subject: [PATCH] area of Mathieson parametrized; hit ctor change eloss to qdc; Helix removed --- HMPID/AliHMPID.cxx | 62 --- HMPID/AliHMPID.h | 10 +- HMPID/AliHMPIDCluster.cxx | 7 +- HMPID/AliHMPIDDigit.cxx | 118 ++--- HMPID/AliHMPIDDigit.h | 158 +++---- HMPID/AliHMPIDHit.cxx | 34 +- HMPID/AliHMPIDHit.h | 51 ++- HMPID/AliHMPIDParam.cxx | 32 +- HMPID/AliHMPIDParam.h | 8 +- HMPID/AliHMPIDReconstructor.cxx | 4 +- HMPID/AliHMPIDSelector.C | 4 +- HMPID/AliHMPIDv1.cxx | 163 ++----- HMPID/AliHMPIDv1.h | 2 +- HMPID/HMPIDbaseLinkDef.h | 4 +- HMPID/Hmenu.C | 782 +++++++++++++++----------------- HMPID/libHMPIDbase.pkg | 2 +- 16 files changed, 641 insertions(+), 800 deletions(-) diff --git a/HMPID/AliHMPID.cxx b/HMPID/AliHMPID.cxx index 83d5a64bc12..b8ac5c7b342 100644 --- a/HMPID/AliHMPID.cxx +++ b/HMPID/AliHMPID.cxx @@ -95,65 +95,3 @@ void AliHMPID::SetTreeAddress() // 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;iPrimNTreeH()->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); -} -//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ diff --git a/HMPID/AliHMPID.h b/HMPID/AliHMPID.h index 83ed19b4d1f..b1514d5f361 100644 --- a/HMPID/AliHMPID.h +++ b/HMPID/AliHMPID.h @@ -28,13 +28,11 @@ public: 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 @@ -42,7 +40,6 @@ public: 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 @@ -50,13 +47,12 @@ public: 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&); diff --git a/HMPID/AliHMPIDCluster.cxx b/HMPID/AliHMPIDCluster.cxx index 9fbb3cde96f..009b0f86b1a 100644 --- a/HMPID/AliHMPIDCluster.cxx +++ b/HMPID/AliHMPIDCluster.cxx @@ -63,8 +63,9 @@ void AliHMPIDCluster::CorrSin() // 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; } //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -117,7 +118,7 @@ void AliHMPIDCluster::Print(Option_t* opt)const } 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() diff --git a/HMPID/AliHMPIDDigit.cxx b/HMPID/AliHMPIDDigit.cxx index 3000c92d1c7..d9a8ab52130 100644 --- a/HMPID/AliHMPIDDigit.cxx +++ b/HMPID/AliHMPIDDigit.cxx @@ -13,19 +13,40 @@ // * provided "as is" without express or implied warranty. * // ************************************************************************** -#include "AliHMPIDDigit.h" //class header -#include //Hit2Sdi() -#include //Draw() +#include "AliHMPIDDigit.h" //class header +#include //WriteRaw() +#include //Draw() +#include //WriteRaw() +#include //WriteRaw() +#include //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. @@ -81,59 +102,46 @@ void AliHMPIDDigit::Draw(Option_t*) 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;iDigGetEntriesFast();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() //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - diff --git a/HMPID/AliHMPIDDigit.h b/HMPID/AliHMPIDDigit.h index 52563a37567..cddb49100a9 100644 --- a/HMPID/AliHMPIDDigit.h +++ b/HMPID/AliHMPIDDigit.h @@ -7,47 +7,40 @@ #include //Mathieson() #include //IsOverTh() #include //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 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(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 @@ -93,39 +82,20 @@ protected: //AliDigit has fTracks[3] 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 { @@ -166,30 +136,52 @@ Float_t AliHMPIDDigit::Mathieson(Float_t x,Float_t y)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(pxkMaxPx) return kTRUE; + if(pykMaxPy) return kTRUE; + + fPad=Abs(ch,pc,px,py);fQ=qdc;fTracks[0]=tid; + return kFALSE; +} #endif diff --git a/HMPID/AliHMPIDHit.cxx b/HMPID/AliHMPIDHit.cxx index 6b289fad1b2..591ee991f09 100644 --- a/HMPID/AliHMPIDHit.cxx +++ b/HMPID/AliHMPIDHit.cxx @@ -14,9 +14,9 @@ // ************************************************************************** #include "AliHMPIDHit.h" //class header -#include "AliHMPIDDigit.h" -#include -#include +#include //Draw() Print() +#include //Draw() +#include //Hit2Sdi() ClassImp(AliHMPIDHit) //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -28,10 +28,30 @@ void AliHMPIDHit::Draw(Option_t*) 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()); @@ -50,8 +70,8 @@ void AliHMPIDHit::Print(Option_t*)const 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":""); } //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ diff --git a/HMPID/AliHMPIDHit.h b/HMPID/AliHMPIDHit.h index e84d8da21f9..d623c444c93 100644 --- a/HMPID/AliHMPIDHit.h +++ b/HMPID/AliHMPIDHit.h @@ -6,29 +6,30 @@ #include //base class #include //ctor #include //QdcTot() -#include +#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] @@ -37,20 +38,30 @@ protected: / 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 diff --git a/HMPID/AliHMPIDParam.cxx b/HMPID/AliHMPIDParam.cxx index 1fa574cfc26..f5bc54b6fb0 100644 --- a/HMPID/AliHMPIDParam.cxx +++ b/HMPID/AliHMPIDParam.cxx @@ -35,7 +35,13 @@ AliHMPIDParam::AliHMPIDParam():TNamed("RichParam","default version") // 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 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -46,6 +52,30 @@ void AliHMPIDParam::Print(Option_t* opt) const 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 diff --git a/HMPID/AliHMPIDParam.h b/HMPID/AliHMPIDParam.h index 65c5c1bac1c..f85eee41fec 100644 --- a/HMPID/AliHMPIDParam.h +++ b/HMPID/AliHMPIDParam.h @@ -29,7 +29,8 @@ public: 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 @@ -55,11 +56,10 @@ AliHMPIDParam* AliHMPIDParam::Instance() // 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 diff --git a/HMPID/AliHMPIDReconstructor.cxx b/HMPID/AliHMPIDReconstructor.cxx index 385ecc8da9c..51eefc49d86 100644 --- a/HMPID/AliHMPIDReconstructor.cxx +++ b/HMPID/AliHMPIDReconstructor.cxx @@ -32,7 +32,7 @@ void AliHMPIDReconstructor::Dig2Clu(TObjArray *pDigAll,TObjArray *pCluAll,Bool_t // 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 @@ -129,7 +129,7 @@ void AliHMPIDReconstructor::Reconstruct(AliRunLoader *pAL,AliRawReader* pRR)cons 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 diff --git a/HMPID/AliHMPIDSelector.C b/HMPID/AliHMPIDSelector.C index 08fed6c3299..bf12db1fd0a 100644 --- a/HMPID/AliHMPIDSelector.C +++ b/HMPID/AliHMPIDSelector.C @@ -166,7 +166,7 @@ void caf() 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"); @@ -323,7 +323,7 @@ void cosmic() 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;iCluGetEntriesFast();iClu++){//clusters loop AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluCh->UncheckedAt(iClu); pCluQ->Fill(pClu->Q()); diff --git a/HMPID/AliHMPIDv1.cxx b/HMPID/AliHMPIDv1.cxx index 1df76788d1b..988d6d50a43 100644 --- a/HMPID/AliHMPIDv1.cxx +++ b/HMPID/AliHMPIDv1.cxx @@ -15,12 +15,10 @@ #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 //Hits2SDigits() -#include #include //StepManager() for gMC #include //StepHistory() #include //StepManager(),Hits2SDigits() @@ -29,16 +27,12 @@ #include #include #include //StepManager() -#include //Digits2Raw() -#include //Digits2Raw() #include //CreateMaterials() #include //CreateMaterials() #include //CreateGeometry() -#include //Optics() -#include //Optics() -#include //Optics() -#include //Optics() -#include //CreateMaterials() +#include //DefineOpticalProperties() +#include //DefineOpticalProperties() +#include //IsLostByFresnel() #include //CreateMaterials() #include //CreateMaterials() @@ -46,10 +40,10 @@ ClassImp(AliHMPIDv1) //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 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)); } //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -91,60 +85,7 @@ void AliHMPIDv1::CreateMaterials() AliMaterial(++matId,"W" ,aW ,zW ,dW ,radW ,absW ); AliMedium(kW ,"W" , matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin); AliMaterial(++matId,"Al" ,aAl ,zAl ,dAl ,radAl ,absAl ); AliMedium(kAl ,"Al" , matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin); - 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;iSetMarkerStyle(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() @@ -158,24 +99,9 @@ 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); } @@ -329,13 +255,13 @@ Bool_t AliHMPIDv1::IsLostByFresnel() 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]; @@ -403,7 +329,7 @@ void AliHMPIDv1::GenFee(Int_t iCh,Float_t eloss) //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 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."); @@ -413,8 +339,8 @@ void AliHMPIDv1::Hits2SDigits() if(!GetLoader()->TreeH()) {GetLoader()->LoadHits(); } if(!GetLoader()->TreeS()) {GetLoader()->MakeTree("S"); MakeBranch("S");}//to - for(Int_t iPrimN=0;iPrimNTreeH()->GetEntries();iPrimN++){//prims loop - GetLoader()->TreeH()->GetEntry(iPrimN); + for(Int_t iEnt=0;iEntTreeH()->GetEntries();iEnt++){//prims loop + GetLoader()->TreeH()->GetEntry(iEnt); Hit2Sdi(Hits(),SdiLst()); }//prims loop GetLoader()->TreeS()->Fill(); @@ -428,23 +354,19 @@ void AliHMPIDv1::Hits2SDigits() //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 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;iHitGetEntries();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."); @@ -456,31 +378,8 @@ void AliHMPIDv1::Digits2Raw() } 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;iGetEntriesFast();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;iUnloadDigits(); AliDebug(1,"Stop."); }//Digits2Raw() @@ -638,31 +537,31 @@ void AliHMPIDv1::StepManager() Int_t pid= gMC->TrackPid(); //take PID Float_t etot= gMC->Etot(); //total hpoton energy, [GeV] Double_t x[3]; gMC->TrackPosition(x[0],x[1],x[2]); //take MARS position at entrance to PC - Float_t 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() //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ diff --git a/HMPID/AliHMPIDv1.h b/HMPID/AliHMPIDv1.h index e9b2f51a71a..f56422db8c3 100644 --- a/HMPID/AliHMPIDv1.h +++ b/HMPID/AliHMPIDv1.h @@ -27,7 +27,7 @@ public: 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 diff --git a/HMPID/HMPIDbaseLinkDef.h b/HMPID/HMPIDbaseLinkDef.h index 14c0d9b3d29..c4c2cd501d2 100644 --- a/HMPID/HMPIDbaseLinkDef.h +++ b/HMPID/HMPIDbaseLinkDef.h @@ -3,14 +3,14 @@ #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+; diff --git a/HMPID/Hmenu.C b/HMPID/Hmenu.C index 09bca024277..df24aea3f3a 100644 --- a/HMPID/Hmenu.C +++ b/HMPID/Hmenu.C @@ -2,8 +2,9 @@ AliRun *a; AliRunLoader *al; TGeoManager *g; //globals for easy manual man 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() { @@ -18,18 +19,21 @@ void Hmenu() { 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(" ","" ,""); @@ -47,8 +51,8 @@ void Hmenu() 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(); @@ -56,15 +60,15 @@ void General() //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 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() @@ -82,14 +86,14 @@ 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() //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -101,25 +105,11 @@ void geo ( ) {gGeoManager->GetTopVolume()->Draw("ogl");} 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() { @@ -130,7 +120,8 @@ 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); @@ -175,399 +166,80 @@ void PrintMap() 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;iParticleNGetNtrack();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;iEvtGetEvent(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;iCluGetEntriesFast();iClu++){ - AliHMPIDCluster *pClu = (AliHMPIDCluster*)pCluLst->At(iClu); - Int_t cfm=0; for(Int_t iDig=0;iDigSize();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;iEvtNGetRunLoader()->GetEvent(iEvtN)) return; - AliStack *pStack = GetLoader()->GetRunLoader()->Stack(); - for(Int_t iPrimN=0;iPrimNTreeH()->GetEntries();iPrimN++){//prims loop - GetLoader()->TreeH()->GetEntry(iPrimN); - for(Int_t iHitN=0;iHitNGetEntries();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;iDigGetEntries();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;iGetNtrack();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;iPcSetFillColor(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 @@ -582,6 +254,10 @@ void hed() static TFile *pEsdFl=0; static TTree *pEsdTr=0; static AliESD *pEsd=0; + + + + if(!pC&&iEvtIsOpen()) 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(iEvtGetEntry(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."); @@ -619,6 +302,19 @@ void Close_hed() pC=0x0; } //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +void ReadHits(TClonesArray *pHitLst) +{ + pHitLst->Delete(); Int_t cnt=0; + for(Int_t iEnt=0;iEntTreeH()->GetEntries();iEnt++){ //TreeH loop + hl->TreeH()->GetEntry(iEnt); //get current entry (prim) + for(Int_t iHit=0;iHitHits()->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() { @@ -640,7 +336,7 @@ void sed() AliESD esd; - EsdFromStack(&esd); + SimEsd(&esd); HitsFromEsd(&esd,&lh); @@ -651,10 +347,10 @@ void sed() 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; @@ -689,12 +385,19 @@ void DrawEvt(TCanvas *pC,TObjArray *pDigLst,TObjArray *pCluLst,AliESD *pEsd) } gPad->SetEditable(kTRUE); gPad->Clear(); DrawPc(0); + for(Int_t iHit=0;iHitGetEntries();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)); @@ -702,6 +405,15 @@ void DrawEvt(TCanvas *pC,TObjArray *pDigLst,TObjArray *pCluLst,AliESD *pEsd) // 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; @@ -734,6 +446,21 @@ void t1(Int_t case=1) 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(); @@ -760,12 +487,231 @@ void HitsFromEsd(AliESD *pEsd, TClonesArray *pHitLst) 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 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;iHitRndm()*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;iCluGetEntriesFast();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;iPrimTreeH()->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(); +} diff --git a/HMPID/libHMPIDbase.pkg b/HMPID/libHMPIDbase.pkg index 65c4057cda7..13847cf1b8e 100644 --- a/HMPID/libHMPIDbase.pkg +++ b/HMPID/libHMPIDbase.pkg @@ -1,4 +1,4 @@ -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 -- 2.39.3