From 998b831fa16a477b96e57c7f856058baa61ed918 Mon Sep 17 00:00:00 2001 From: kir Date: Sat, 25 Sep 2004 09:48:53 +0000 Subject: [PATCH] New Geom, support for ESD- major change --- RICH/AliGenRadioactive.cxx | 45 ++ RICH/AliGenRadioactive.h | 28 ++ RICH/AliRICH.cxx | 530 ++++++++------------- RICH/AliRICH.h | 32 +- RICH/AliRICHChamber.cxx | 46 +- RICH/AliRICHChamber.h | 58 ++- RICH/AliRICHClusterFinder.cxx | 53 +-- RICH/AliRICHDisplFast.cxx | 8 +- RICH/AliRICHHelix.cxx | 12 + RICH/AliRICHHelix.h | 99 ++++ RICH/AliRICHParam.cxx | 17 +- RICH/AliRICHParam.h | 114 ++--- RICH/AliRICHRecon.cxx | 850 +++------------------------------- RICH/AliRICHRecon.h | 65 +-- RICH/AliRICHReconstructor.cxx | 89 ++-- RICH/AliRICHReconstructor.h | 29 +- RICH/AliRICHTracker.cxx | 72 +++ RICH/AliRICHTracker.h | 28 ++ RICH/AliRICHv0.cxx | 120 ++--- RICH/AliRICHv0.h | 1 - RICH/AliRICHv1.cxx | 45 +- RICH/CreateConfig.C | 32 +- RICH/Geom.C | 98 ++-- RICH/Opticals.h | 9 +- RICH/RICHLinkDef.h | 4 + RICH/RichBatch.C | 26 +- RICH/api.txt | 28 +- RICH/libRICH.pkg | 2 +- RICH/menu.C | 107 ++++- 29 files changed, 1061 insertions(+), 1586 deletions(-) create mode 100644 RICH/AliGenRadioactive.cxx create mode 100644 RICH/AliGenRadioactive.h create mode 100644 RICH/AliRICHHelix.cxx create mode 100644 RICH/AliRICHHelix.h create mode 100644 RICH/AliRICHTracker.cxx create mode 100644 RICH/AliRICHTracker.h diff --git a/RICH/AliGenRadioactive.cxx b/RICH/AliGenRadioactive.cxx new file mode 100644 index 00000000000..1244b351dc5 --- /dev/null +++ b/RICH/AliGenRadioactive.cxx @@ -0,0 +1,45 @@ +#include "AliGenRadioactive.h" +#include +#include +#include +#include + +//__________________________________________________________________________________________________ +AliGenRadioactive::AliGenRadioactive(Int_t iSrcNucleus,Int_t iNsecondaries):AliGenerator() +{ +//Main ctor. Used to define radioactive source. + Double_t e[100],a[100];//arrays to store experimental points + Int_t nPoints; + switch(iSrcNucleus){ + case kSr90: fPartId=kElectron; nPoints=46; //experimental part + a[ 0]=0.08605; a[ 1]=0.0878; a[ 2]=0.08705; a[ 3]=0.07855; a[ 4]=0.0709; a[ 5]=0.0647; a[ 6]=0.05015; a[ 7]=0.0372; a[ 8]=0.0268; a[ 9]=0.0215; + a[10]=0.0157; a[11]=0.01685; a[12]=0.01745; a[13]=0.01645; a[14]=0.0175; a[15]=0.01635; a[16]=0.01825; a[17]=0.0177; a[18]=0.01735; a[19]=0.0161; + a[20]=0.0159; a[21]=0.0176; a[22]=0.01605; a[23]=0.0161; a[24]=0.01495;a[25]=0.01595; a[26]=0.01525; a[27]=0.0138; a[28]=0.0121; a[29]=0.0101; + a[30]=0.01175; a[31]=0.01095; a[32]=0.0089; a[33]=0.0091; a[34]=0.0625; a[35]=0.0505; a[36]=0.0475; a[37]=0.0039; a[38]=0.0031; a[39]=0.0028; + a[40]=0.0025; a[41]=0.0017; a[42]=4.5e-4; a[43]=4.5e-4; a[44]=1.5e-4; a[45]=0; + break; + default: AliError("Wrong source nucleus specified"); return; + } + for(Int_t i=0;iFill(e[i],a[i]); + fNpart=iNsecondaries; +} +//__________________________________________________________________________________________________ +void AliGenRadioactive::Generate() +{ +// Generate one trigger + Int_t nt=0; + Double_t ekin=0,p=0,theta=0,phi=0,x=0,y=0,z=0,px=0,py=0,pz=0,polx=0,poly=0,polz=0; + Double_t m=gAlice->PDGDB()->GetParticle(fPartId)->Mass(); + for(Int_t i=0;iGetRandom(); p=TMath::Sqrt(ekin*(2*m+ekin)); + theta=Rndm()*fThetaMax; phi=Rndm()*fPhiMax; + px=p*TMath::Cos(theta)*TMath::Cos(phi); py=p*TMath::Cos(theta)*TMath::Sin(phi); pz=p*TMath::Sin(theta); +// AliDebug(1,Form("Origin=(%5.2f,%5.2f,%5.2f) Ekin=%5.3fMeV,P=%5.3fMeV (%5.2f,%5.2f,%5.2f)",x,y,z,1000*ekin,1000*p,1000*px,1000*py,1000*pz)); + PushTrack(fTrackIt,-1,fPartId,px,py,pz,ekin+m, + x, y, z,0, + polx=0,poly=0,polz=0,kPPrimary,nt);//cm, GeV + } +} diff --git a/RICH/AliGenRadioactive.h b/RICH/AliGenRadioactive.h new file mode 100644 index 00000000000..0a257f494b5 --- /dev/null +++ b/RICH/AliGenRadioactive.h @@ -0,0 +1,28 @@ +#ifndef AliGenRadioactive_h +#define AliGenRadioactive_h +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +// Simple radioactive source generator. +// Currently Sr90 only + +#include +#include + +typedef enum {kSr90} RadioNuclei_t; + +class AliGenRadioactive : public AliGenerator +{ + public: + AliGenRadioactive() : AliGenerator(),fPartId(0),fGenH1(0) {;}//default ctor + AliGenRadioactive(Int_t iSource,Int_t iNparts); //main ctor + virtual ~AliGenRadioactive() {if(fGenH1) delete fGenH1;}//dtor + virtual void Generate(); //interface from AliGenerator, generates current event + + +protected: + Int_t fPartId; //ID of secondary from radioactive nucleus + TH1F *fGenH1; //Histogram containg exp spectrum to be used to sample the energy + ClassDef(AliGenRadioactive,1) // Radioactive source generator +}; +#endif diff --git a/RICH/AliRICH.cxx b/RICH/AliRICH.cxx index 72e34ed217f..842bd21840a 100644 --- a/RICH/AliRICH.cxx +++ b/RICH/AliRICH.cxx @@ -35,21 +35,22 @@ #include #include #include +#include ClassImp(AliRICHhit) //__________________________________________________________________________________________________ void AliRICHhit::Print(Option_t*)const { - ::Info("hit","Ch=%1i, TID=%6i, eloss=%9.3f eV, in-out dist=%9.4f, OUT(%7.2f,%7.2f,%7.2f)" - ,fChamber,fTrack,fEloss*1e9,Length(),fOutX3.X(),fOutX3.Y(),fOutX3.Z()); + AliInfo(Form("Ch=%1i, TID=%6i, eloss=%9.3f eV, in-out dist=%9.4f, OUT(%7.2f,%7.2f,%7.2f)" + ,fChamber,fTrack,fEloss*1e9,Length(),fOutX3.X(),fOutX3.Y(),fOutX3.Z())); } //__________________________________________________________________________________________________ ClassImp(AliRICHdigit) //__________________________________________________________________________________________________ void AliRICHdigit::Print(Option_t*)const { - ::Info("digit","cfm=%9i, cs=%2i, x=%3i, y=%3i, q=%8.3f, TID1=%5i, TID2=%5i, TID3=%5i", - fCFM,fChamber,fPadX,fPadY,fQdc,fTracks[0],fTracks[1],fTracks[2]); + AliInfo(Form("cfm=%9i, cs=%2i, x=%3i, y=%3i, q=%8.3f, TID1=%5i, TID2=%5i, TID3=%5i", + fCFM,fChamber,fPadX,fPadY,fQdc,fTracks[0],fTracks[1],fTracks[2])); } //__________________________________________________________________________________________________ ClassImp(AliRICHcluster) @@ -66,8 +67,8 @@ void AliRICHcluster::Print(Option_t*)const ::Info("cluster","cfm=%10i, cs=%2i, SiMa=%6i, Shape=%5i, x=%7.3f, y=%7.3f, Q=%6i, %s with %i digits", fCFM,fChamber,fSize,fShape,fX,fY,fQdc,status,fDigits->GetEntriesFast()); else - ::Info("cluster","cfm=%10i, cs=%2i, SiMa=%6i, Shape=%5i, x=%7.3f, y=%7.3f, Q=%6i, %s with %i digits", - fCFM,fChamber,fSize,fShape,fX,fY,fQdc,status,0); + AliInfo(Form("cfm=%10i, cs=%2i, SiMa=%6i, Shape=%5i, x=%7.3f, y=%7.3f, Q=%6i, %s with %i digits", + fCFM,fChamber,fSize,fShape,fX,fY,fQdc,status,0)); } //__________________________________________________________________________________________________ @@ -93,17 +94,17 @@ AliRICH::AliRICH(const char *name, const char *title) :AliDetector(name,title),fpParam(new AliRICHParam),fSdigits(0),fNsdigits(0),fDigitsNew(0),fClusters(0) { //Named ctor - if(GetDebug())Info("named ctor","Start."); + AliDebug(1,"Start."); //AliDetector ctor deals with Hits and Digits (reset them to 0, does not create them) CreateHits(); gAlice->GetMCApp()->AddHitList(fHits); fCounters.ResizeTo(20); fCounters.Zero(); - if(GetDebug())Info("named ctor","Stop."); + AliDebug(1,"Stop."); }//AliRICH::AliRICH(const char *name, const char *title) //__________________________________________________________________________________________________ AliRICH::~AliRICH() { //dtor - if(GetDebug()) Info("dtor","Start."); + AliDebug(1,"Start."); if(fpParam) delete fpParam; @@ -112,13 +113,13 @@ AliRICH::~AliRICH() if(fDigits) delete fDigits; if(fDigitsNew) {fDigitsNew->Delete(); delete fDigitsNew;} if(fClusters) {fClusters->Delete(); delete fClusters;} - if(GetDebug()) Info("dtor","Stop."); + AliDebug(1,"Stop."); }//AliRICH::~AliRICH() //__________________________________________________________________________________________________ void AliRICH::Hits2SDigits() { // Create a list of sdigits corresponding to list of hits. Every hit generates one or more sdigits. - if(GetDebug()) Info("Hit2SDigits","Start."); + AliDebug(1,"Start."); for(Int_t iEventN=0;iEventNGetRunLoader()->GetAliRun()->GetEventsPerRun();iEventN++){//events loop GetLoader()->GetRunLoader()->GetEvent(iEventN);//get next event @@ -130,16 +131,16 @@ void AliRICH::Hits2SDigits() GetLoader()->TreeH()->GetEntry(iPrimN); for(Int_t iHitN=0;iHitNGetEntries();iHitN++){//hits loop AliRICHhit *pHit=(AliRICHhit*)Hits()->At(iHitN);//get current hit - TVector2 x2 = C(pHit->C())->Glob2Loc(pHit->OutX3());//hit position in the chamber local system + TVector2 x2 = C(pHit->C())->Mrs2Pc(pHit->OutX3());//hit position in photocathode plane Int_t iTotQdc=P()->TotQdc(x2,pHit->Eloss());//total charge produced by hit, 0 if hit in dead zone if(iTotQdc==0) continue; TVector area=P()->Loc2Area(x2);//determine affected pads, dead zones analysed inside - if(GetDebug()) Info("Hits2SDigits","hit(%6.2f,%6.2f)->area(%3.0f,%3.0f)-(%3.0f,%3.0f) QDC=%4i",x2.X(),x2.Y(),area[0],area[1],area[2],area[3],iTotQdc); + AliDebug(1,Form("hit(%6.2f,%6.2f)->area(%3.0f,%3.0f)-(%3.0f,%3.0f) QDC=%4i",x2.X(),x2.Y(),area[0],area[1],area[2],area[3],iTotQdc)); TVector pad(2); for(pad[1]=area[1];pad[1]<=area[3];pad[1]++)//affected pads loop for(pad[0]=area[0];pad[0]<=area[2];pad[0]++){ Double_t padQdc=iTotQdc*P()->FracQdc(x2,pad); - if(GetDebug()) Info("Hits2SDigits","current pad(%3.0f,%3.0f) with QDC =%6.2f",pad[0],pad[1],padQdc); + AliDebug(1,Form("current pad(%3.0f,%3.0f) with QDC =%6.2f",pad[0],pad[1],padQdc)); if(padQdc>0.1) AddSDigit(pHit->C(),pad,padQdc,GetLoader()->GetRunLoader()->Stack()->Particle(pHit->GetTrack())->GetPdgCode(),pHit->GetTrack()); }//affected pads loop }//hits loop @@ -150,127 +151,135 @@ void AliRICH::Hits2SDigits() }//events loop GetLoader()->UnloadHits(); GetLoader()->GetRunLoader()->UnloadHeader(); GetLoader()->GetRunLoader()->UnloadKinematics(); GetLoader()->UnloadSDigits(); - if(GetDebug()) Info("Hit2SDigits","Stop."); + AliDebug(1,"Stop."); }//Hits2SDigits() //__________________________________________________________________________________________________ void AliRICH::BuildGeometry() { //Builds a TNode geometry for event display - if(GetDebug())Info("BuildGeometry","Start."); + AliInfo("Start."); TNode *node, *subnode, *top; top=gAlice->GetGeometry()->GetNode("alice"); - - new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15); - Float_t wid=P()->SectorSizeX(); - Float_t len=P()->SectorSizeY(); - new TBRIK("PHOTO","PHOTO","void",wid/2,0.1,len/2); - + Float_t widx=P()->SectorSizeX(); + Float_t leny=P()->SectorSizeY(); + Float_t dz = P()->Zfreon()+P()->Zwin()+P()->Pc2Win(); + Float_t dead = P()->DeadZone(); + + new TBRIK("RICH","RICH","void",widx+dead/2,leny+leny/2+dead,dz+0.1); //RICH chamber + new TBRIK("RPC" ,"RPC" ,"void",widx/2,leny/2,0.01); //RICH sector + for(int i=1;i<=kNchambers;i++){ top->cd(); - node = new TNode(Form("RICH%i",i),Form("RICH%i",i),"S_RICH",C(i)->X(),C(i)->Y(),C(i)->Z(),C(i)->RotMatrixName()); + node = new TNode(Form("RICH%i",i),Form("RICH%i",i),"RICH",C(i)->Center().X(),C(i)->Center().Y(),C(i)->Center().Z(),C(i)->RotMatrixName()); node->SetLineColor(kRed); node->cd(); - subnode = new TNode("PHOTO1","PHOTO1","PHOTO",wid+P()->DeadZone(),5,len/2+P()->DeadZone()/2,""); + subnode = new TNode("PHOTO1","PHOTO1","RPC",-widx/2-dead/2,-leny-dead/2,dz,""); subnode->SetLineColor(kGreen); fNodes->Add(subnode); - subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,len/2+P()->DeadZone()/2,""); + subnode = new TNode("PHOTO1","PHOTO1","RPC", widx/2+dead/2,-leny-dead/2,dz,""); subnode->SetLineColor(kGreen); fNodes->Add(subnode); - subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-wid-P()->DeadZone(),5,len/2+P()->DeadZone()/2,""); + subnode = new TNode("PHOTO1","PHOTO1","RPC",-widx/2-dead/2, 0,dz,""); subnode->SetLineColor(kGreen); fNodes->Add(subnode); - subnode = new TNode("PHOTO1","PHOTO1","PHOTO",wid+P()->DeadZone(),5,-len/2-P()->DeadZone()/2,""); + subnode = new TNode("PHOTO1","PHOTO1","RPC", widx/2+dead/2, 0,dz,""); subnode->SetLineColor(kGreen); fNodes->Add(subnode); - subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-len/2 -P()->DeadZone()/2,""); + subnode = new TNode("PHOTO1","PHOTO1","RPC",-widx/2-dead/2, leny+dead/2,dz,""); subnode->SetLineColor(kGreen); fNodes->Add(subnode); - subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-wid-P()->DeadZone(),5,-len/2 - P()->DeadZone()/2,""); + subnode = new TNode("PHOTO1","PHOTO1","RPC", widx/2+dead/2, leny+dead/2,dz,""); subnode->SetLineColor(kGreen); fNodes->Add(subnode); fNodes->Add(node); - } - if(GetDebug())Info("BuildGeometry","Stop."); + } + + AliDebug(1,"Stop."); }//void AliRICH::BuildGeometry() //______________________________________________________________________________ void AliRICH::CreateMaterials() { // Definition of available RICH materials -#include "Opticals.h" - Float_t a=0,z=0,den=0,radl=0,absl=0; + Int_t material=0; //tmp material id number + Float_t a=0,z=0,den=0,radl=0,absl=0; //tmp material parameters Float_t tmaxfd=-10.0, deemax=-0.2, stemax=-0.1,epsil=0.001, stmin=-0.001; Int_t isxfld = gAlice->Field()->Integ(); Float_t sxmgmx = gAlice->Field()->Max(); - Int_t material; - AliMaterial(material=1, "Air $",a=14.61,z=7.3, den=0.001205,radl=30420.0,absl=67500);//(Air) - AliMedium(1, "DEFAULT MEDIUM AIR$",material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); + Float_t aAir[4]={12.,14.,16.,36.}; Float_t zAir[4]={6.,7.,8.,18.}; Float_t wAir[4]={0.000124,0.755267,0.231781,0.012827};//total 0.9999999 + AliMixture(++material, "RichAir",aAir,zAir,den=0.00120479,4,wAir); //1 (Air) 0.01% C 75% N 23% O 1% Ar + AliMedium(kAir, "RichAir",material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); - AliMaterial( 6, "HON", a=12.01,z=6.0, den=0.1, radl=18.8, absl=0); //(C)-equivalent radl - AliMedium(2, "HONEYCOMB$", 6, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); + AliMaterial(++material, "RichRohacell", a=12.01,z=6.0, den=0.1, radl=18.8, absl=0); //2 Rohacell 51 C-equiv radl rad cover + AliMedium(kRoha, "RichRohacell", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); - AliMaterial(16, "CSI", a=12.01,z=6.0, den=0.1, radl=18.8, absl=0); //CsI-radl equivalent - AliMedium(kCSI, "CSI$", 16, 1, isxfld, sxmgmx,tmaxfd, stemax, deemax, epsil, stmin); + Float_t aQuartz[2]={28.09,16.0}; Float_t zQuartz[2]={14.00, 8.0}; Float_t wQuartz[2]={1,2}; + AliMixture(++material, "RichSiO2",aQuartz,zQuartz,den=2.64,-2, wQuartz); //3 Quarz (SiO2) -trasparent rad window + AliMedium(kSiO2, "RichSiO2",material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); - AliMaterial(11, "GRI", a=63.54,z=29.0,den=8.96, radl=1.43, absl=0); //anode grid (Cu) - AliMedium(7, "GRID$", 11, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); + Float_t aFreon[2]={12,19}; Float_t zFreon[2]={6,9}; Float_t wmatFreon[2]={6,14}; + AliMixture(++material, "RichC6F14",aFreon,zFreon,den=1.68,-2,wmatFreon); //4 Freon (C6F14) + AliMedium(kC6F14, "RichC6F14",material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); + + Float_t aMethane[2]={12.01,1}; Float_t zMethane[2]={6,1}; Float_t wMethane[2]={1,4}; + AliMixture (++material, "RichCH4", aMethane, zMethane, den=7.17e-4,-2, wMethane); //5,9 methane (CH4) normal and for Gap + AliMedium(kCH4, "RichCH4" , material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); + AliMixture (++material, "RichCH4", aMethane, zMethane, den=7.17e-4,-2, wMethane); //5,9 methane (CH4) normal and for Gap + AliMedium(kGap, "RichCH4gap$", material, 1, isxfld, sxmgmx, tmaxfd, 0.1 , -deemax, epsil, -stmin); + + AliMaterial(++material, "RichCsI", a=12.01,z=6.0, den=0.1, radl=18.8, absl=0); //6 CsI-radl equivalent + AliMedium(kCsI, "RichCsI$", material, 1, isxfld, sxmgmx,tmaxfd, stemax, deemax, epsil, stmin); - AliMaterial(50, "ALUM", a=26.98,z=13.0,den=2.699, radl=8.9, absl=0); //aluminium sheet (Al) - AliMedium(10, "ALUMINUM$", 50, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); + AliMaterial(++material, "GridCu", a=63.54,z=29.0,den=8.96, radl=1.43, absl=0); //7 anode grid (Cu) + AliMedium(kGridCu, "GRID$", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); - AliMaterial(material=31, "COPPER$", a=63.54,z=29.0,den=8.96, radl=1.4, absl=0); //(Cu) - AliMedium(12, "PCB_COPPER", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); + AliMixture (++material, "OpSiO2",aQuartz, zQuartz, den=2.64, -2, wQuartz); //8 Quarz (SiO2) - opaque + AliMedium(kOpSiO2, "QUARTZO$",material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); - Float_t aQuartz[2]={28.09,16.0}; Float_t zQuartz[2]={14.00, 8.0}; Float_t wmatQuartz[2]={1,2}; - AliMixture (20, "QUA",aQuartz,zQuartz,den=2.64,-2, wmatQuartz);//Quarz (SiO2) - trasnparent - AliMedium(3, "QUARTZ$", 20, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); + AliMaterial(++material, "ALUM", a=26.98,z=13.0,den=2.699, radl=8.9, absl=0); //10 aluminium sheet (Al) + AliMedium(kAl, "ALUMINUM$", material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); - AliMixture (21, "QUAO",aQuartz, zQuartz, den=2.64, -2, wmatQuartz);//Quarz (SiO2) - opaque - AliMedium(8, "QUARTZO$", 21, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); + Float_t aGlass[5]={12.01,28.09,16,10.8,23}; Float_t zGlass[5]={6,14,8,5,11}; Float_t wGlass[5]={0.5,0.105,0.355,0.03,0.01}; + AliMixture(++material,"GLASS",aGlass, zGlass, den=1.74, 5, wGlass); //11 Glass 50%-C 10.5%-Si 35.5%-O 3%-B 1%-Na + AliMedium(kGlass, "GLASS", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); - Float_t aFreon[2]={12,19}; Float_t zFreon[2]={6,9}; Float_t wmatFreon[2]={6,14}; - AliMixture (material=30, "C6F14",aFreon,zFreon,den=1.68,-2,wmatFreon);//Freon (C6F14) - AliMedium(4, "C6F14$",material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); + AliMaterial(++material, "COPPER$", a=63.54,z=29.0,den=8.96, radl=1.4, absl=0); //12 Cu + AliMedium(kCu, "PCB_COPPER", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); + + + AliMaterial(++material, "W$", a=183.84,z=74.0,den=19.3, radl=0.35, absl=185.0/den); //13 W - anod wires + AliMedium(kW, "W", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); - Float_t aMethane[2]={12.01,1}; Float_t zMethane[2]={6,1}; Float_t wmatMethane[2]={1,4}; - AliMixture (material=40, "CH4", aMethane, zMethane, den=7.17e-4,-2, wmatMethane);//methane (CH4) - AliMedium(kCH4, "CH4$" , material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); - AliMedium(kGAP, "CH4GAP$", material, 1, isxfld, sxmgmx,tmaxfd, 0.1, -deemax, epsil, -stmin); if(P()->IsRadioSrc()){ - AliMaterial(material=45, "STEEL$", a=63.54,z=29.0,den=8.96, radl=1.4, absl=0); //Steel + AliInfo("Special radioactive source materials"); + AliMaterial(++material, "STEEL$", a=63.54,z=29.0,den=8.96, radl=1.4, absl=0); //14 Steel AliMedium(kSteel, "STEEL", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); - AliMaterial(material=46, "PERPEX$", a=63.54,z=29.0,den=8.96, radl=1.4, absl=0); //Perpex + AliMaterial(++material, "PERPEX$", a=63.54,z=29.0,den=8.96, radl=1.4, absl=0); //15 Perpex AliMedium(kPerpex, "PERPEX", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); - AliMaterial(material=47, "STRONZIUM$", a=87.62,z=38.0,den=2.54, radl=4.24, absl=0); //Sr90 + AliMaterial(++material, "STRONZIUM$", a=87.62,z=38.0,den=2.54, radl=4.24, absl=0); //16 Sr90 AliMedium(kSr90, "STRONZIUM", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); } - Float_t aGlass[5]={12.01, 28.09, 16., 10.8, 23.}; - Float_t zGlass[5]={ 6., 14., 8., 5., 11.}; - Float_t wGlass[5]={ 0.5, 0.105, 0.355, 0.03, 0.01}; - AliMixture(material=32, "GLASS",aGlass, zGlass, den=1.74, 5, wGlass);//Glass 50%C+10.5%Si+35.5%O+3% + 1% - AliMedium(11, "GLASS", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); - - Int_t *idtmed = fIdtmed->GetArray()-999; - gMC->SetCerenkov(idtmed[1000], kNbins, aPckov, aAbsCH4, aQeAll, aIdxCH4); - gMC->SetCerenkov(idtmed[1001], kNbins, aPckov, aAbsCH4, aQeAll, aIdxCH4); - gMC->SetCerenkov(idtmed[1002], kNbins, aPckov, aAbsSiO2, aQeAll, aIdxSiO2); - gMC->SetCerenkov(idtmed[1003], kNbins, aPckov, aAbsC6F14, aQeAll, aIdxC6F14); - gMC->SetCerenkov(idtmed[1004], kNbins, aPckov, aAbsCH4, aQeAll, aIdxCH4); - gMC->SetCerenkov(idtmed[1005], kNbins, aPckov, aAbsCsI, aQeCsI, aIdxCH4); - gMC->SetCerenkov(idtmed[1006], kNbins, aPckov, aAbsGrid, aQeAll, aIdxGrid); - gMC->SetCerenkov(idtmed[1007], kNbins, aPckov, aAbsOpSiO2, aQeAll, aIdxOpSiO2); - gMC->SetCerenkov(idtmed[1008], kNbins, aPckov, aAbsCH4, aQeAll, aIdxCH4); - gMC->SetCerenkov(idtmed[1009], kNbins, aPckov, aAbsGrid, aQeAll, aIdxGrid); - gMC->SetCerenkov(idtmed[1010], kNbins, aPckov, aAbsOpSiO2, aQeAll, aIdxOpSiO2); - +//Optical properties: +#include "Opticals.h" + gMC->SetCerenkov((*fIdtmed)[kAir] , kNbins, aPckov, aAbsCH4, aQeAll, aIdxCH4); //1 Air + gMC->SetCerenkov((*fIdtmed)[kRoha] , kNbins, aPckov, aAbsCH4, aQeAll, aIdxCH4); //2 Honeycomb + gMC->SetCerenkov((*fIdtmed)[kSiO2] , kNbins, aPckov, aAbsSiO2, aQeAll, aIdxSiO2); //3 Quartz SiO2 + gMC->SetCerenkov((*fIdtmed)[kC6F14] , kNbins, aPckov, aAbsC6F14, aQeAll, aIdxC6F14); //4 Freon C6F14 + gMC->SetCerenkov((*fIdtmed)[kCH4] , kNbins, aPckov, aAbsCH4, aQeAll, aIdxCH4); //5 Methane CH4 + gMC->SetCerenkov((*fIdtmed)[kCsI] , kNbins, aPckov, aAbsCsI, aQeCsI, aIdxCH4); //6 CsI + gMC->SetCerenkov((*fIdtmed)[kGridCu] , kNbins, aPckov, aAbsGrid, aQeAll, aIdxGrid); //7 grid Cu + gMC->SetCerenkov((*fIdtmed)[kOpSiO2] , kNbins, aPckov, aAbsOpSiO2, aQeAll, aIdxOpSiO2); //8 Opaque quartz SiO2 + gMC->SetCerenkov((*fIdtmed)[kGap] , kNbins, aPckov, aAbsCH4, aQeAll, aIdxCH4); //9 Special methane gap + gMC->SetCerenkov((*fIdtmed)[kAl] , kNbins, aPckov, aAbsGrid, aQeAll, aIdxGrid); //10 Aluminium + gMC->SetCerenkov((*fIdtmed)[kGlass] , kNbins, aPckov, aAbsOpSiO2, aQeAll, aIdxOpSiO2); //11 Glass }//void AliRICH::CreateMaterials() //__________________________________________________________________________________________________ Float_t AliRICH::Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)const @@ -365,7 +374,7 @@ Float_t AliRICH::AbsoCH4(Float_t x)const void AliRICH::MakeBranch(Option_t* option) { //Create Tree branches for the RICH. - if(GetDebug())Info("MakeBranch","Start with option= %s.",option); + AliDebug(1,Form("Start with option= %s.",option)); const Int_t kBufferSize = 4000; @@ -395,29 +404,29 @@ void AliRICH::MakeBranch(Option_t* option) for(Int_t i=0;iTreeR(),Form("%sClusters%d",GetName(),i+1), &((*fClusters)[i]), kBufferSize, 0); }//R - if(GetDebug())Info("MakeBranch","Stop."); + AliDebug(1,"Stop."); }//void AliRICH::MakeBranch(Option_t* option) //__________________________________________________________________________________________________ void AliRICH::SetTreeAddress() { //Set branch address for the Hits and Digits Tree. - if(GetDebug())Info("SetTreeAddress","Start."); + AliDebug(1,"Start."); TBranch *branch; if(fLoader->TreeH()){//H - if(GetDebug())Info("SetTreeAddress","tree H is requested."); + AliDebug(1,"tree H is requested."); CreateHits();//branch map will be in AliDetector::SetTreeAddress }//H AliDetector::SetTreeAddress();//this is after TreeH because we need to guarantee that fHits array is created if(fLoader->TreeS()){//S - if(GetDebug())Info("SetTreeAddress","tree S is requested."); + AliDebug(1,"tree S is requested."); branch=fLoader->TreeS()->GetBranch(GetName()); if(branch){CreateSDigits(); branch->SetAddress(&fSdigits);} }//S if(fLoader->TreeD()){//D - if(GetDebug())Info("SetTreeAddress","tree D is requested."); + AliDebug(1,"tree D is requested."); for(int i=0;iTreeD()->GetBranch(Form("%s%d",GetName(),i+1)); if(branch){CreateDigits(); branch->SetAddress(&((*fDigitsNew)[i]));} @@ -425,13 +434,13 @@ void AliRICH::SetTreeAddress() }//D if(fLoader->TreeR()){//R - if(GetDebug())Info("SetTreeAddress","tree R is requested."); + AliDebug(1,"tree R is requested."); for(int i=0;iTreeR()->GetBranch(Form("%sClusters%d" ,GetName(),i+1)); if(branch){CreateClusters(); branch->SetAddress(&((*fClusters)[i]));} } }//R - if(GetDebug())Info("SetTreeAddress","Stop."); + AliDebug(1,"Stop."); }//void AliRICH::SetTreeAddress() //__________________________________________________________________________________________________ void AliRICH::Print(Option_t *option)const @@ -445,216 +454,100 @@ void AliRICH::Print(Option_t *option)const void AliRICH::CreateGeometry() { //Creates detailed geometry simulation (currently GEANT volumes tree) - if(GetDebug())Info("CreateGeometry","Start main."); -//Opaque quartz thickness - Float_t oquaThickness = .5; -//CsI dimensions - Float_t pcX=P()->PcSizeX(); - Float_t pcY=P()->PcSizeY(); - - Int_t *idtmed = fIdtmed->GetArray()-999; - - Int_t i; - Float_t zs; - Int_t idrotm[1099]; + AliDebug(1,"Start main."); + Double_t cm=1,mm=0.1*cm,mkm=0.001*cm; Float_t par[3]; - par[0]=68.8;par[1]=13 ;par[2]=70.86; gMC->Gsvolu("RICH","BOX ",(*fIdtmed)[kAl], par,3);//External aluminium box - par[0]=66.3;par[1]=13 ;par[2]=68.35; gMC->Gsvolu("SRIC","BOX ",(*fIdtmed)[kAir],par,3);//Air - par[0]=66.3;par[1]=0.025;par[2]=68.35; gMC->Gsvolu("ALUM","BOX ",(*fIdtmed)[kAl], par,3);//Aluminium sheet -//Air 2 (cutting the lower part of the box) -// par[0]=1.25; par[1] = 3; par[2] = 70.86; gMC->Gsvolu("AIR2", "BOX ", idtmed[1000], par, 3); -//Air 3 (cutting the lower part of the box) -// par[0]=66.3; par[1] = 3; par[2] = 1.2505; gMC->Gsvolu("AIR3", "BOX ", idtmed[1000], par, 3); -//Honeycomb - par[0]=66.3;par[1]=0.188; par[2] = 68.35; gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3); - - //par[0] = 66.5; par[1] = .025; par[2] = 63.1; -//Quartz - par[0]=P()->QuartzWidth()/2;par[1]=P()->QuartzThickness()/2;par[2]=P()->QuartzLength()/2; - gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3); -//Spacers (cylinders) - par[0]=0.;par[1]=.5;par[2]=P()->FreonThickness()/2; gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3); -//Feet (freon slabs supports) - par[0] = .7; par[1] = .3; par[2] = 1.9; gMC->Gsvolu("FOOT", "BOX", idtmed[1009], par, 3); -//Opaque quartz - par[0]=P()->QuartzWidth()/2;par[1]= .2;par[2]=P()->QuartzLength()/2; - gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3); -//Frame of opaque quartz - par[0]=P()->OuterFreonWidth()/2;par[1]=P()->FreonThickness()/2;par[2]=P()->OuterFreonLength()/2; - gMC->Gsvolu("OQF1", "BOX ", idtmed[1007], par, 3); - par[0]=P()->InnerFreonWidth()/2;par[1]=P()->FreonThickness()/2;par[2]=P()->InnerFreonLength()/2; - gMC->Gsvolu("OQF2", "BOX ", idtmed[1007], par, 3); -//Freon - par[0]=P()->OuterFreonWidth()/2 - oquaThickness; - par[1]=P()->FreonThickness()/2; - par[2]=P()->OuterFreonLength()/2 - 2*oquaThickness; - gMC->Gsvolu("FRE1", "BOX ", idtmed[1003], par, 3); - - par[0]=P()->InnerFreonWidth()/2 - oquaThickness; - par[1]=P()->FreonThickness()/2; - par[2]=P()->InnerFreonLength()/2 - 2*oquaThickness; - gMC->Gsvolu("FRE2", "BOX ", idtmed[1003], par, 3); -//Methane - par[0]=pcX/2;par[1]=P()->GapThickness()/2;par[2]=pcY/2; gMC->Gsvolu("META","BOX ",idtmed[1004], par, 3); -//Methane gap - par[0]=pcX/2;par[1]=P()->GapAmp()/2;par[2]=pcY/2;gMC->Gsvolu("GAP ","BOX ",(*fIdtmed)[kGAP],par,3); -//CsI PC - par[0]=pcX/2;par[1]=.25;par[2]=pcY/2; gMC->Gsvolu("CSI ", "BOX ", (*fIdtmed)[kCSI], par, 3); -//Anode grid - par[0] = 0.;par[1] = .001;par[2] = 20.; gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3); - -//Wire supports -//Bar of metal - par[0]=pcX/2;par[1]=1.05;par[2]=1.05; gMC->Gsvolu("WSMe", "BOX ", idtmed[1009], par, 3); -//Ceramic pick up (base) - par[0]=pcX/2;par[1]= .25;par[2]=1.05; gMC->Gsvolu("WSG1", "BOX ", idtmed[1010], par, 3); -//Ceramic pick up (head) - par[0] = pcX/2;par[1] = .1;par[2] = .1; gMC->Gsvolu("WSG2", "BOX ", idtmed[1010], par, 3); - -//Aluminium supports for methane and CsI -//Short bar - par[0]=pcX/2;par[1]=P()->GapThickness()/2 + .25; par[2] = (68.35 - pcY/2)/2; - gMC->Gsvolu("SMSH", "BOX", idtmed[1009], par, 3); -//Long bar - par[0]=(66.3 - pcX/2)/2;par[1]=P()->GapThickness()/2+.25;par[2]=pcY/2+68.35-pcY/2; - gMC->Gsvolu("SMLG", "BOX", idtmed[1009], par, 3); - -//Aluminium supports for freon -//Short bar - par[0] = P()->QuartzWidth()/2; par[1] = .3; par[2] = (68.35 - P()->QuartzLength()/2)/2; - gMC->Gsvolu("SFSH", "BOX", idtmed[1009], par, 3); -//Long bar - par[0] = (66.3 - P()->QuartzWidth()/2)/2; par[1] = .3; - par[2] = P()->QuartzLength()/2 + 68.35 - P()->QuartzLength()/2; - gMC->Gsvolu("SFLG", "BOX", idtmed[1009], par, 3); -//PCB backplane - par[0] = pcX/2;par[1] = .25; par[2] = pcY/4 -.5025; gMC->Gsvolu("PCB ", "BOX", idtmed[1011], par, 3); - -//Backplane supports -//Aluminium slab - par[0] = 33.15;par[1] = 2;par[2] = 21.65; gMC->Gsvolu("BACK", "BOX", idtmed[1009], par, 3); -//Big hole - par[0] = 9.05; par[1] = 2; par[2] = 4.4625; gMC->Gsvolu("BKHL", "BOX", idtmed[1000], par, 3); -//Small hole - par[0] = 5.7;par[1] = 2;par[2] = 4.4625; gMC->Gsvolu("BKHS", "BOX", idtmed[1000], par, 3); -//Place holes inside backplane support - gMC->Gspos("BKHS", 1, "BACK", .8 + 5.7,0., .6 + 4.4625, 0, "ONLY"); - gMC->Gspos("BKHS", 2, "BACK", -.8 - 5.7,0., .6 + 4.4625, 0, "ONLY"); - gMC->Gspos("BKHS", 3, "BACK", .8 + 5.7,0., -.6 - 4.4625, 0, "ONLY"); - gMC->Gspos("BKHS", 4, "BACK", -.8 - 5.7,0., -.6 - 4.4625, 0, "ONLY"); - gMC->Gspos("BKHS", 5, "BACK", .8 + 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY"); - gMC->Gspos("BKHS", 6, "BACK", -.8 - 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY"); - gMC->Gspos("BKHS", 7, "BACK", .8 + 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY"); - gMC->Gspos("BKHS", 8, "BACK", -.8 - 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY"); - gMC->Gspos("BKHL", 1, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., .6 + 4.4625, 0, "ONLY"); - gMC->Gspos("BKHL", 2, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 4.4625, 0, "ONLY"); - gMC->Gspos("BKHL", 3, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 4.4625, 0, "ONLY"); - gMC->Gspos("BKHL", 4, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 4.4625, 0, "ONLY"); - gMC->Gspos("BKHL", 5, "BACK", .8 + 11.4+ 1.6 + 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY"); - gMC->Gspos("BKHL", 6, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY"); - gMC->Gspos("BKHL", 7, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY"); - gMC->Gspos("BKHL", 8, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY"); -//Place material inside RICH - gMC->Gspos("SRIC", 1, "RICH", 0.,0., 0., 0, "ONLY"); -// gMC->Gspos("AIR2", 1, "RICH", 66.3 + 1.2505, 1.276-P()->GapThickness()/2-P()->QuartzThickness()-P()->FreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 0., 0, "ONLY"); -// gMC->Gspos("AIR2", 2, "RICH", -66.3 - 1.2505,1.276-P()->GapThickness()/2-P()->QuartzThickness()-P()->FreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 0., 0, "ONLY"); -// gMC->Gspos("AIR3", 1, "RICH", 0., 1.276-P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, -68.35 - 1.25, 0, "ONLY"); -// gMC->Gspos("AIR3", 2, "RICH", 0., 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 68.35 + 1.25, 0, "ONLY"); - gMC->Gspos("ALUM", 1, "SRIC", 0., 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .6 - .05 - .376 -.025, 0., 0, "ONLY"); - gMC->Gspos("HONE", 1, "SRIC", 0., 1.276- P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .6 - .05 - .188, 0., 0, "ONLY"); - gMC->Gspos("ALUM", 2, "SRIC", 0., 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .6 - .025, 0., 0, "ONLY"); - gMC->Gspos("FOOT", 1, "SRIC", 64.95, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, 36.9, 0, "ONLY"); - gMC->Gspos("FOOT", 2, "SRIC", 21.65, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3 , 36.9, 0, "ONLY"); - gMC->Gspos("FOOT", 3, "SRIC", -21.65, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, 36.9, 0, "ONLY"); - gMC->Gspos("FOOT", 4, "SRIC", -64.95, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, 36.9, 0, "ONLY"); - gMC->Gspos("FOOT", 5, "SRIC", 64.95, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, -36.9, 0, "ONLY"); - gMC->Gspos("FOOT", 6, "SRIC", 21.65, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, -36.9, 0, "ONLY"); - gMC->Gspos("FOOT", 7, "SRIC", -21.65, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, -36.9, 0, "ONLY"); - gMC->Gspos("FOOT", 8, "SRIC", -64.95, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, -36.9, 0, "ONLY"); - gMC->Gspos("OQUA", 1, "SRIC", 0., 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .2, 0., 0, "ONLY"); -// Methane supports - gMC->Gspos("SMLG", 1, "SRIC", pcX/2 + (66.3 - pcX/2)/2, 1.276 + .25, 0., 0, "ONLY"); - gMC->Gspos("SMLG", 2, "SRIC", - pcX/2 - (66.3 - pcX/2)/2, 1.276 + .25, 0., 0, "ONLY"); - gMC->Gspos("SMSH", 1, "SRIC", 0., 1.276 + .25, pcY/2 + (68.35 - pcY/2)/2, 0, "ONLY"); - gMC->Gspos("SMSH", 2, "SRIC", 0., 1.276 + .25, - pcY/2 - (68.35 - pcY/2)/2, 0, "ONLY"); -//Freon supports - Float_t suppY = 1.276 - P()->GapThickness()/2- P()->QuartzThickness() -P()->FreonThickness() - .2 + .3; //y position of freon supports - gMC->Gspos("SFLG", 1, "SRIC", P()->QuartzWidth()/2 + (66.3 - P()->QuartzWidth()/2)/2, suppY, 0., 0, "ONLY"); - gMC->Gspos("SFLG", 2, "SRIC", - P()->QuartzWidth()/2 - (66.3 - P()->QuartzWidth()/2)/2, suppY, 0., 0, "ONLY"); - gMC->Gspos("SFSH", 1, "SRIC", 0., suppY, P()->QuartzLength()/2 + (68.35 - P()->QuartzLength()/2)/2, 0, "ONLY"); - gMC->Gspos("SFSH", 2, "SRIC", 0., suppY, - P()->QuartzLength()/2 - (68.35 - P()->QuartzLength()/2)/2, 0, "ONLY"); - AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.); -//Place spacers - Int_t nspacers = 30; - for (i = 0; i < nspacers/3; i++) { - zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2; - gMC->Gspos("SPAC", i, "FRE1", 10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings - } - for (i = nspacers/3; i < (nspacers*2)/3; i++) { - zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2; - gMC->Gspos("SPAC", i, "FRE1", 0, 0., zs, idrotm[1019], "ONLY"); //Original settings - } - for (i = (nspacers*2)/3; i < nspacers; ++i) { - zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2; - gMC->Gspos("SPAC", i, "FRE1", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings - } - for (i = 0; i < nspacers/3; i++) { - zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2; - gMC->Gspos("SPAC", i, "FRE2", 10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings - } - for (i = nspacers/3; i < (nspacers*2)/3; i++) { - zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2; - gMC->Gspos("SPAC", i, "FRE2", 0, 0., zs, idrotm[1019], "ONLY"); //Original settings - } - for (i = (nspacers*2)/3; i < nspacers; ++i) { - zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2; - gMC->Gspos("SPAC", i, "FRE2", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings - } - gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY"); - gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY"); - gMC->Gspos("OQF1", 1, "SRIC", P()->OuterFreonWidth()/2 + P()->InnerFreonWidth()/2 + 2, 1.276 - P()->GapThickness()/2- P()->QuartzThickness() -P()->FreonThickness()/2, 0., 0, "ONLY"); //Original settings (31.3) - gMC->Gspos("OQF2", 2, "SRIC", 0., 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()/2, 0., 0, "ONLY"); //Original settings - gMC->Gspos("OQF1", 3, "SRIC", - (P()->OuterFreonWidth()/2 + P()->InnerFreonWidth()/2) - 2, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()/2, 0., 0, "ONLY"); //Original settings (-31.3) - gMC->Gspos("QUAR", 1, "SRIC", 0., 1.276 - P()->GapThickness()/2 - P()->QuartzThickness()/2, 0., 0, "ONLY"); - gMC->Gspos("GAP ", 1, "META", 0., P()->GapThickness()/2 - P()->GapAmp()/2 - 0.0001, 0., 0, "ONLY"); - gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY"); - gMC->Gspos("CSI ", 1, "SRIC", 0., 1.276 + P()->GapThickness()/2 + .25, 0., 0, "ONLY"); -//Wire support placing - gMC->Gspos("WSG2", 1, "GAP ", 0., P()->GapAmp()/2 - .1, 0., 0, "ONLY"); - gMC->Gspos("WSG1", 1, "CSI ", 0., 0., 0., 0, "ONLY"); - gMC->Gspos("WSMe", 1, "SRIC ", 0., 1.276 + P()->GapThickness()/2 + .5 + 1.05, 0., 0, "ONLY"); -//Backplane placing - gMC->Gspos("BACK", 1, "SRIC ", -33.15, 1.276 + P()->GapThickness()/2 + .5 + 2.1 + 2, 43.3, 0, "ONLY"); - gMC->Gspos("BACK", 2, "SRIC ", 33.15, 1.276 + P()->GapThickness()/2 + .5 + 2.1 + 2 , 43.3, 0, "ONLY"); - gMC->Gspos("BACK", 3, "SRIC ", -33.15, 1.276 + P()->GapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY"); - gMC->Gspos("BACK", 4, "SRIC ", 33.15, 1.276 + P()->GapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY"); - gMC->Gspos("BACK", 5, "SRIC ", 33.15, 1.276 + P()->GapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY"); - gMC->Gspos("BACK", 6, "SRIC ", -33.15, 1.276 + P()->GapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY"); -//PCB placing - gMC->Gspos("PCB ", 1, "SRIC ", 0., 1.276 + P()->GapThickness()/2 + .5 + 1.05, pcX/4 + .5025 + 2.5, 0, "ONLY"); - gMC->Gspos("PCB ", 2, "SRIC ", 0., 1.276 + P()->GapThickness()/2 + .5 + 1.05, -pcX/4 - .5025 - 2.5, 0, "ONLY"); - + Int_t id=0; + //place chambers into mother volume ALIC - for(int i=1;i<=kNchambers;i++){ - AliMatrix(idrotm[1000+i],C(i)->ThetaXd(),C(i)->PhiXd(), - C(i)->ThetaYd(),C(i)->PhiYd(), - C(i)->ThetaZd(),C(i)->PhiZd()); - gMC->Gspos("RICH",i,"ALIC",C(i)->X(),C(i)->Y(),C(i)->Z(),idrotm[1000+i], "ONLY"); - } + par[0]=(6*mm+1681*mm+6*mm)/2;par[1]=(6*mm+1466*mm+6*mm)/2;par[2]=(80*mm+40*mm)*2/2;gMC->Gsvolu("RICH","BOX ",(*fIdtmed)[kCH4],par,3);//2033P1 z p84 TDR + if(P()->IsRadioSrc()){ + AliInfo("Special test beam geometry"); + gMC->Gspos("RICH",1,"ALIC",0,0,0,0, "ONLY"); //single RICH chamber without rotation test beam geometry + }else + for(int i=1;i<=kNchambers;i++){ //full RICH with 7 chambers + AliMatrix(id,C(i)->ThetaXd(),C(i)->PhiXd(), C(i)->ThetaYd(),C(i)->PhiYd(), C(i)->ThetaZd(),C(i)->PhiZd()); + gMC->Gspos("RICH",i,"ALIC",C(i)->Center().X(),C(i)->Center().Y(),C(i)->Center().Z(),id, "ONLY"); + } +//Pad Panel frame 6 sectors + par[0]=648*mm/2;par[1]= 411*mm/2;par[2]=40 *mm/2;gMC->Gsvolu("RPPF","BOX ",(*fIdtmed)[kAl] ,par,3);//PPF 2001P2 inner size of the slab by 1mm more + par[0]=181*mm/2;par[1]=89.25*mm/2;par[2]=38.3*mm/2;gMC->Gsvolu("PPFL","BOX ",(*fIdtmed)[kAir] ,par,3);//large whole + par[0]=114*mm/2;par[1]=89.25*mm/2;par[2]=38.3*mm/2;gMC->Gsvolu("PPFS","BOX ",(*fIdtmed)[kAir] ,par,3);//small whole + par[0]=644*mm/2;par[1]= 407*mm/2;par[2]= 1.7*mm/2;gMC->Gsvolu("RPC ","BOX ",(*fIdtmed)[kCsI] ,par,3);//by 0.2 mm more then actual size (PCB 2006P1) + + gMC->Gspos("RPPF",1,"RICH", -335*mm, -433*mm, 8*cm+20*mm, 0,"ONLY");//F1 2040P1 z p.84 TDR + gMC->Gspos("RPPF",2,"RICH", +335*mm, -433*mm, 8*cm+20*mm, 0,"ONLY"); + gMC->Gspos("RPPF",3,"RICH", -335*mm, 0*mm, 8*cm+20*mm, 0,"ONLY"); + gMC->Gspos("RPPF",4,"RICH", +335*mm, 0*mm, 8*cm+20*mm, 0,"ONLY"); + gMC->Gspos("RPPF",5,"RICH", -335*mm, +433*mm, 8*cm+20*mm, 0,"ONLY"); + gMC->Gspos("RPPF",6,"RICH", +335*mm, +433*mm, 8*cm+20*mm, 0,"ONLY"); + gMC->Gspos("RPC ",1,"RPPF", 0*mm, 0*mm, -19.15*mm, 0,"ONLY");//PPF 2001P2 + gMC->Gspos("PPFL",1,"RPPF", -224.5*mm, -151.875*mm, 0.85*mm, 0,"ONLY"); + gMC->Gspos("PPFL",2,"RPPF", -224.5*mm, - 50.625*mm, 0.85*mm, 0,"ONLY"); + gMC->Gspos("PPFL",3,"RPPF", -224.5*mm, + 50.625*mm, 0.85*mm, 0,"ONLY"); + gMC->Gspos("PPFL",4,"RPPF", -224.5*mm, +151.875*mm, 0.85*mm, 0,"ONLY"); + gMC->Gspos("PPFS",1,"RPPF", - 65.0*mm, -151.875*mm, 0.85*mm, 0,"ONLY"); + gMC->Gspos("PPFS",2,"RPPF", - 65.0*mm, - 50.625*mm, 0.85*mm, 0,"ONLY"); + gMC->Gspos("PPFS",3,"RPPF", - 65.0*mm, + 50.625*mm, 0.85*mm, 0,"ONLY"); + gMC->Gspos("PPFS",4,"RPPF", - 65.0*mm, +151.875*mm, 0.85*mm, 0,"ONLY"); + gMC->Gspos("PPFS",5,"RPPF", + 65.0*mm, -151.875*mm, 0.85*mm, 0,"ONLY"); + gMC->Gspos("PPFS",6,"RPPF", + 65.0*mm, - 50.625*mm, 0.85*mm, 0,"ONLY"); + gMC->Gspos("PPFS",7,"RPPF", + 65.0*mm, + 50.625*mm, 0.85*mm, 0,"ONLY"); + gMC->Gspos("PPFS",8,"RPPF", + 65.0*mm, +151.875*mm, 0.85*mm, 0,"ONLY"); + gMC->Gspos("PPFL",5,"RPPF", +224.5*mm, -151.875*mm, 0.85*mm, 0,"ONLY"); + gMC->Gspos("PPFL",6,"RPPF", +224.5*mm, - 50.625*mm, 0.85*mm, 0,"ONLY"); + gMC->Gspos("PPFL",7,"RPPF", +224.5*mm, + 50.625*mm, 0.85*mm, 0,"ONLY"); + gMC->Gspos("PPFL",8,"RPPF", +224.5*mm, +151.875*mm, 0.85*mm, 0,"ONLY"); +//Gap - anod wires 6 copies to RICH + par[0]=648*mm/2;par[1]= 411*mm/2 ;par[2]=4.45*mm/2;gMC->Gsvolu("RGAP ","BOX ",(*fIdtmed)[kCH4] ,par,3);//xy as PPF 2001P2 z WP 2099P1 + par[0]= 0*mm ;par[1]= 20*mkm/2 ;par[2]= 648*mm/2;gMC->Gsvolu("RANO","TUBE",(*fIdtmed)[kW] ,par,3);//WP 2099P1 z = gap x PPF 2001P2 + AliMatrix(id,180,0, 90,90, 90,0); //wires along x + + gMC->Gspos("RGAP",1,"RICH", -335*mm, -433*mm,8*cm-2.225*mm, 0,"ONLY"); //F1 2040P1 z WP 2099P1 + gMC->Gspos("RGAP",2,"RICH", +335*mm, -433*mm,8*cm-2.225*mm, 0,"ONLY"); + gMC->Gspos("RGAP",3,"RICH", -335*mm, 0*mm,8*cm-2.225*mm, 0,"ONLY"); + gMC->Gspos("RGAP",4,"RICH", +335*mm, 0*mm,8*cm-2.225*mm, 0,"ONLY"); + gMC->Gspos("RGAP",5,"RICH", -335*mm, +433*mm,8*cm-2.225*mm, 0,"ONLY"); + gMC->Gspos("RGAP",6,"RICH", +335*mm, +433*mm,8*cm-2.225*mm, 0,"ONLY"); + for(int i=1;i<=96;i++) + gMC->Gspos("RANO",i,"RGAP", 0*mm, -411/2*mm+i*4*mm, 0.185*mm, id,"ONLY"); //WP 2099P1 +//Radiator 3 modules + par[0]=1330*mm/2 ;par[1]= 413*mm/2 ;par[2]= 24*mm/2; gMC->Gsvolu("RRAD","BOX ",(*fIdtmed)[kC6F14] ,par,3); // Rad 2011P1 + par[0]=1330*mm/2 ;par[1]= 413*mm/2 ;par[2]= 4*mm/2; gMC->Gsvolu("RRFR","BOX ",(*fIdtmed)[kRoha] ,par,3); //front + par[0]=1330*mm/2 ;par[1]= 413*mm/2 ;par[2]= 5*mm/2; gMC->Gsvolu("RRWI","BOX ",(*fIdtmed)[kSiO2] ,par,3); //window + par[0]=1330*mm/2 ;par[1]= 5*mm/2 ;par[2]= 15*mm/2; gMC->Gsvolu("RRLO","BOX ",(*fIdtmed)[kRoha] ,par,3); //long side + par[0]= 10*mm/2 ;par[1]= 403*mm/2 ;par[2]= 15*mm/2; gMC->Gsvolu("RRSH","BOX ",(*fIdtmed)[kRoha] ,par,3); //short side + par[0]= 0 ;par[1]= 10*mm/2 ;par[2]= 15*mm/2; gMC->Gsvolu("RRSP","TUBE",(*fIdtmed)[kSiO2] ,par,3); //spacer + + gMC->Gspos("RRAD",1,"RICH", 0*mm,-434*mm, -12*mm, 0,"ONLY"); //radiator to RICH + gMC->Gspos("RRAD",2,"RICH", 0*mm, 0*mm, -12*mm, 0,"ONLY"); + gMC->Gspos("RRAD",3,"RICH", 0*mm,+434*mm, -12*mm, 0,"ONLY"); + gMC->Gspos("RRFR",1,"RRAD", 0*mm, 0*mm, -10.0*mm, 0,"ONLY"); //front cover + gMC->Gspos("RRWI",1,"RRAD", 0*mm, 0*mm, 9.5*mm, 0,"ONLY"); //quartz window (back cover) + gMC->Gspos("RRLO",1,"RRAD", 0*mm,-204*mm, -0.5*mm, 0,"ONLY"); //long side + gMC->Gspos("RRLO",2,"RRAD", 0*mm,+204*mm, -0.5*mm, 0,"ONLY"); //long side + gMC->Gspos("RRSH",1,"RRAD",-660*mm, 0*mm, -0.5*mm, 0,"ONLY"); //short side + gMC->Gspos("RRSH",2,"RRAD",+660*mm, 0*mm, -0.5*mm, 0,"ONLY"); //short side + for(int i=0;i<3;i++) + for(int j=0;j<10;j++) + gMC->Gspos("RRSP",10*i+j,"RRAD",-1330*mm/2+116*mm+j*122*mm,(i-1)*105*mm,-0.5*mm,0,"ONLY"); - - if(GetDebug())Info("CreateGeometry","Stop main."); +//Sandbox + par[0]=1419*mm/2 ;par[1]=1378*mm/2;par[2]=50.5*mm/2; gMC->Gsvolu("RSNB","BOX ",(*fIdtmed)[kAir] ,par,3); //2072P1 + par[0]=1419*mm/2 ;par[1]=1378*mm/2;par[2]= 0.5*mm/2; gMC->Gsvolu("RSCO","BOX ",(*fIdtmed)[kAl] ,par,3); //cover + par[0]=1359*mm/2 ;par[1]=1318*mm/2;par[2]=49.5*mm/2; gMC->Gsvolu("RSHO","BOX ",(*fIdtmed)[kRoha] ,par,3); //honeycomb structure + + gMC->Gspos("RSNB",1,"RICH", 0*mm, 0*mm, -73.75*mm, 0,"ONLY"); //p.84 TDR sandbox to rich + gMC->Gspos("RSHO",1,"RSNB", 0*mm, 0*mm, 0*mm, 0,"ONLY"); //2072P1 honeycomv to sandbox + gMC->Gspos("RSCO",1,"RSNB", 0*mm, 0*mm, +25*mm, 0,"ONLY"); //cover to sandbox + gMC->Gspos("RSCO",2,"RSNB", 0*mm, 0*mm, -25*mm, 0,"ONLY"); //cover to sandbox + + AliDebug(1,"Stop main."); }//void AliRICH::CreateGeometry() //__________________________________________________________________________________________________ -void AliRICH::Reconstruct()const -{ - AliRICHClusterFinder finder(const_cast(this)); - finder.Exec(); -} -//__________________________________________________________________________________________________ void AliRICH::ControlPlots() { -//Creates a set of hists to control the results of simulation +// Creates a set of hists to control the results of simulation. Hists are in file $HOME/RCP.root TH1F *pHxD=0,*pHyD=0,*pNumClusH1=0, *pQdcH1=0, *pSizeH1=0, @@ -667,14 +560,14 @@ void AliRICH::ControlPlots() Bool_t isDig =!GetLoader()->LoadDigits(); Bool_t isClus=!GetLoader()->LoadRecPoints(); - if(!isDig && !isClus){Error("ControlPlots","No digits and clusters! Nothing to do.");return;} + if(!isDig && !isClus){AliError("No digits and clusters! Nothing to do.");return;} TStopwatch sw;TDatime time; TFile *pFile = new TFile("$(HOME)/RCP.root","RECREATE"); if(isDig){ - cout<<"Digits available\n"; + AliInfo("Digits available"); pHxD=new TH1F("HitDigitDiffX","Hit-Digits diff X all chambers;diff [cm]",100,-10,10); pHyD=new TH1F("HitDigitDiffY","Hit-Digits diff Y all chambers;diff [cm]",100,-10,10); }//isDig @@ -742,8 +635,8 @@ void AliRICH::ControlPlots() if(isDig){ for(Int_t iDigN=0;iDigNDigits(iChamN)->GetEntries();iDigN++){//digits loop AliRICHdigit *pDig=(AliRICHdigit*)R()->Digits(iChamN)->At(iDigN); - AliRICHhit *pHit=Hit(pDig->GetTrack(0)); - TVector2 hitV2=R()->C(iChamN)->Glob2Loc(pHit->OutX3()); TVector2 digV2=R()->P()->Pad2Loc(pDig->Pad()); + AliRICHhit *pHit=Hit(pDig->GetTrack(0));//get first hit of this digit + TVector2 hitV2=R()->C(iChamN)->Mrs2Pc(pHit->OutX3()); TVector2 digV2=R()->P()->Pad2Loc(pDig->Pad());//center of pad for digit pHxD->Fill(hitV2.X()-digV2.X()); pHyD->Fill(hitV2.Y()-digV2.Y()); }//digits loop }//isDig @@ -776,7 +669,7 @@ AliRICHhit* AliRICH::Hit(Int_t tid) void AliRICH::PrintHits(Int_t iEvtN) { //Prints a list of RICH hits for a given event. Default is event number 0. - Info("PrintHits","List of RICH hits for event %i",iEvtN); + AliInfo(Form("List of RICH hits for event %i",iEvtN)); R()->GetLoader()->GetRunLoader()->GetEvent(iEvtN); if(R()->GetLoader()->LoadHits()) return; @@ -787,7 +680,7 @@ void AliRICH::PrintHits(Int_t iEvtN) iTotalHits+=R()->Hits()->GetEntries(); } R()->GetLoader()->UnloadHits(); - Info("PrintHits","totally %i hits",iTotalHits); + AliInfo(Form("totally %i hits",iTotalHits)); } //__________________________________________________________________________________________________ void AliRICH::PrintSDigits(Int_t iEvtN) @@ -837,40 +730,13 @@ void AliRICH::PrintClusters(Int_t iEvtN) Info("PrintClusters","totally %i clusters",iTotalClusters); } //__________________________________________________________________________________________________ -void AliRICH::FillESD(AliESD *pESD)const +void AliRICH::PrintTracks(Int_t iEvtN) { -//This methode fills AliESDtrack with information from RICH - Info("FillESD","Start with %i tracks",pESD->GetNumberOfTracks()); - const Double_t masses[5]={0.000511,0.105658,0.139567,0.493677,0.93828};//electron,muon,pion,kaon,proton - const Double_t refIndex = 1.29052; - - Double_t thetaExp = 0.7; - Double_t thetaTh[5]; - Double_t sinThetaThNorm; - Double_t sigmaThetaTh[5]; - Double_t height[5]; - Double_t totalHeight=0; - - for(Int_t iTrackN=0;iTrackNGetNumberOfTracks();iTrackN++){//ESD tracks loop - AliESDtrack *pTrack = pESD->GetTrack(iTrackN); - - pTrack->Print(""); -// TVector2 x2=P()->HelixCross(pTrack);//returns cross point of track with RICH PC in LRS - Double_t pmod = pTrack->GetP(); - - for(Int_t iPart=4;iPart>=0;iPart--){ - Double_t cosThetaTh = TMath::Sqrt(masses[iPart]*masses[iPart]+pmod*pmod)/(refIndex*pmod); - if(cosThetaTh>=1) {break;} - thetaTh[iPart] = TMath::ACos(cosThetaTh); - sinThetaThNorm = TMath::Sin(thetaTh[iPart])/TMath::Sqrt(1-1/(refIndex*refIndex)); - sigmaThetaTh[iPart] = (0.014*(1/sinThetaThNorm-1) + 0.0043)*1.25; - height[iPart] = TMath::Gaus(thetaExp,thetaTh[iPart],sigmaThetaTh[iPart]); - totalHeight +=height[iPart]; - } - - pTrack->SetRICHsignal(thetaExp); - Double_t richPID[5]; - for(Int_t iPart=0;iPart<5;iPart++) richPID[iPart] = height[iPart]/totalHeight; - pTrack->SetRICHpid(richPID); - }//ESD tracks loop +//prints a list of tracks (including secondaru) for a given event + Info("PrintTracks","List of all tracks for event %i",iEvtN); + R()->GetLoader()->GetRunLoader()->GetEvent(iEvtN); + if(R()->GetLoader()->GetRunLoader()->LoadHeader()) return; + if(R()->GetLoader()->GetRunLoader()->LoadKinematics()) return; + Int_t iTracksCounter=0; + Info("PrintTracks","totally %i tracks",iTracksCounter); } diff --git a/RICH/AliRICH.h b/RICH/AliRICH.h index bba4b603e05..6ab61c68ac2 100644 --- a/RICH/AliRICH.h +++ b/RICH/AliRICH.h @@ -77,9 +77,9 @@ class AliRICHcluster :public TObject public: enum ClusterStatus {kEdge,kShape,kSize,kRaw,kResolved,kEmpty=kBad}; AliRICHcluster():TObject(),fCFM(0),fSize(0),fShape(0),fQdc(0),fChamber(0),fX(0),fY(0),fStatus(kEmpty),fDigits(0) {} - virtual ~AliRICHcluster() {Reset();} + virtual ~AliRICHcluster() {AliDebug(1,"Start");/*Reset();*/} void Reset() {DeleteDigits();fCFM=fSize=fShape=fQdc=fChamber=0;fX=fY=0;fStatus=kEmpty;} //cleans the cluster - void DeleteDigits() {if(fDigits) delete fDigits; fDigits=0;} //deletes the list of digits + void DeleteDigits() {if(fDigits) {delete fDigits;} fDigits=0;} //deletes the list of digits AliRICHcluster& operator=(const AliRICHcluster&) {return *this;} Int_t Nlocals() const{return fSize-10000*(fSize/10000);} //number of local maximums Int_t Size() const{return fSize/10000;} //number of digits in cluster @@ -108,11 +108,12 @@ public: Int_t CombiPid() const{return fCFM;} // void CFM(Int_t c,Int_t f,Int_t m) {fCFM=1000000*c+1000*f+m;} //cluster contributors TObjArray* Digits() const{return fDigits;} // - virtual void Print(Option_t *option="")const; // - inline void AddDigit(AliRICHdigit *pDig); // - inline void CoG(Int_t nLocals); //calculates center of gravity - void Fill(AliRICHcluster *pRaw,Double_t x,Double_t y,Double_t q,Int_t cfm) //form new resolved cluster from raw one - {fCFM=cfm;fChamber=pRaw->Fchamber();fSize=pRaw->Fsize();fQdc=(Int_t)(q*pRaw->Q());fX=x;fY=y;fStatus=kResolved;} // + virtual void Print(Option_t *option="")const; // + inline void AddDigit(AliRICHdigit *pDig); // + inline void CoG(Int_t nLocals); //calculates center of gravity + void Fill(AliRICHcluster *pRaw,Double_t x,Double_t y,Double_t q,Int_t cfm) //form new resolved cluster from raw one + {fCFM=cfm;fChamber=pRaw->Fchamber();fSize=pRaw->Fsize();fQdc=(Int_t)(q*pRaw->Q());fX=x;fY=y;fStatus=kResolved;} + Double_t DistTo(TVector2 x) const{return TMath::Sqrt((x.X()-fX)*(x.X()-fX)+(x.Y()-fY)*(x.Y()-fY));} //distance to given point protected: Int_t fCFM; //1000000*Ncerenkovs+1000*Nfeedbacks+Nmips Int_t fSize; //10000*(how many digits belong to this cluster) + nLocalMaxima @@ -171,9 +172,8 @@ public: virtual void StepManager() =0; //interface from AliMC virtual void Hits2SDigits(); //interface from AliSimulation virtual AliDigitizer* CreateDigitizer(AliRunDigitizer* man) const {return new AliRICHDigitizer(man);} //interface from AliSimulation - virtual void Reconstruct() const; //interface from AliReconstruction - virtual void FillESD(AliESD *pESD) const; //interface from AliReconstruction - virtual void Print(Option_t *option) const; //prints current RICH status +// virtual void Reconstruct() const; //interface from AliReconstruction +// virtual void FillESD(AliESD *pESD) const; //interface from AliReconstruction virtual void SetTreeAddress(); //interface from AliLoader virtual void MakeBranch(Option_t *opt=" "); //interface from AliLoader virtual void CreateMaterials(); //interface from AliMC @@ -198,10 +198,12 @@ public: AliRICH* R() {return this;} //provides pointer to RICH main object TVector Counters() const{return fCounters;} //provides a set of counters void ControlPlots(); //utility + virtual void Print(Option_t *option="") const; //prints current RICH status void PrintHits (Int_t iEvent=0); //utility void PrintSDigits (Int_t iEvent=0); //utility void PrintDigits (Int_t iEvent=0); //utility void PrintClusters(Int_t iEvent=0); //utility + void PrintTracks (Int_t iEvent=0); //utility void AddHit(Int_t c,Int_t tid,TVector3 i3,TVector3 o3,Double_t eloss=0){TClonesArray &tmp=*fHits;new(tmp[fNhits++])AliRICHhit(c,tid,i3,o3,eloss);} inline void AddSDigit(Int_t c,TVector pad,Double_t q,Int_t pid,Int_t tid); @@ -213,7 +215,7 @@ public: {Int_t c=cl.C()-1;TClonesArray &tmp=*((TClonesArray*)fClusters->At(c));new(tmp[fNclusters[c]++])AliRICHcluster(cl);} AliRICHhit* Hit(Int_t tid); //returns pointer ot RICH hit for a given tid protected: - enum {kAir=1,kCSI=6,kGAP=9,kAl=10,kCH4=5,kSteel=15,kPerpex=16,kSr90=17}; + enum {kAir=1,kRoha,kSiO2,kC6F14,kCH4,kCsI,kGridCu,kOpSiO2,kGap,kAl,kGlass,kCu,kW,kSteel,kPerpex,kSr90}; AliRICHParam *fpParam; //main RICH parametrization //fHits and fDigits belong to AliDetector TClonesArray *fSdigits; //! list of sdigits @@ -233,21 +235,21 @@ protected: void AliRICH::CreateHits() { if(fHits) return; - if(GetDebug())Info("CreateHits","creating hits container."); + AliDebug(1,"creating hits container."); fHits=new TClonesArray("AliRICHhit",10000); fNhits=0; } //__________________________________________________________________________________________________ void AliRICH::CreateSDigits() { if(fSdigits) return; - if(GetDebug())Info("CreateSDigits","creating sdigits container."); + AliDebug(1,"creating sdigits container."); fSdigits=new TClonesArray("AliRICHdigit",10000); fNsdigits=0; } //__________________________________________________________________________________________________ void AliRICH::CreateDigits() { if(fDigitsNew) return; - if(GetDebug())Info("CreateDigits","creating digits containers."); + AliDebug(1,"creating digits containers."); fDigitsNew = new TObjArray(kNchambers); for(Int_t i=0;iAddAt(new TClonesArray("AliRICHdigit",10000), i); fNdigitsNew[i]=0;} } @@ -255,7 +257,7 @@ void AliRICH::CreateDigits() void AliRICH::CreateClusters() { if(fClusters) return; - if(GetDebug())Info("CreateClusters","creating clusters containers."); + AliDebug(1,"creating clusters containers."); fClusters = new TObjArray(kNchambers); for(Int_t i=0;iAddAt(new TClonesArray("AliRICHcluster",10000), i); fNclusters[i]=0;} } diff --git a/RICH/AliRICHChamber.cxx b/RICH/AliRICHChamber.cxx index 7f7f86f3205..c945d3afeb9 100644 --- a/RICH/AliRICHChamber.cxx +++ b/RICH/AliRICHChamber.cxx @@ -15,47 +15,50 @@ #include "AliRICHChamber.h" #include +#include ClassImp(AliRICHChamber) //______________________________________________________________________________ AliRICHChamber::AliRICHChamber(Int_t iChamber):TNamed() { //main ctor. Defines all geometry parameters for the given module. - SetToZenith();//put to up position +// 7 6 ^ +// 5 4 3 | +// 2 1 | +// <-----z y . x +// horizontal angle between chambers 19.5 grad +// vertical angle between chambers 20 grad + RotY(90);//rotate around y + fCenterX3.SetXYZ(490,0,0);fPcX3.SetXYZ(490+8,0,0);fRadX3.SetXYZ(490-2,0,0); //shift center along x by 490 cm + switch(iChamber){ - case 1: - RotateX(-AliRICHParam::AngleYZ()); - RotateZ(-AliRICHParam::AngleXY()); - fName="RICHc1";fTitle="RICH chamber 1"; + case 0: //special test beam configuration without rotation. + break; + case 1: + RotY(19.5); RotZ(-20); //right and down break; case 2: - RotateZ(-AliRICHParam::AngleXY()); - fName="RICHc2";fTitle="RICH chamber 2"; + RotZ(-20); //down break; case 3: - RotateX(-AliRICHParam::AngleYZ()); - fName="RICHc3";fTitle="RICH chamber 3"; + RotY(19.5); //right break; case 4: - fName="RICHc4";fTitle="RICH chamber 4"; //no turns - break; + break; //no rotation case 5: - RotateX(AliRICHParam::AngleYZ()); - fName="RICHc5";fTitle="RICH chamber 5"; + RotY(-19.5); //left break; case 6: - RotateZ(AliRICHParam::AngleXY()); - fName="RICHc6";fTitle="RICH chamber 6"; + RotZ(20); //up break; case 7: - RotateX(AliRICHParam::AngleYZ()); - RotateZ(AliRICHParam::AngleXY()); - fName="RICHc7";fTitle="RICH chamber 7"; + RotY(-19.5); RotZ(20); //left and up break; default: Fatal("named ctor","Wrong chamber number %i, check CreateChamber ctor",iChamber); }//switch(iModuleN) - RotateZ(AliRICHParam::AngleRot());//apply common rotation + fName=Form("RICHc%i",iChamber);fTitle=Form("RICH chamber %i",iChamber); + RotZ(30); //apply common rotation fpRotMatrix=new TRotMatrix("rot"+fName,"rot"+fName, Rot().ThetaX()*TMath::RadToDeg(), Rot().PhiX()*TMath::RadToDeg(), Rot().ThetaY()*TMath::RadToDeg(), Rot().PhiY()*TMath::RadToDeg(), Rot().ThetaZ()*TMath::RadToDeg(), Rot().PhiZ()*TMath::RadToDeg()); @@ -64,9 +67,6 @@ AliRICHChamber::AliRICHChamber(Int_t iChamber):TNamed() void AliRICHChamber::Print(Option_t *opt) const { // Debug printout -// printf("%s r=%8.3f theta=%5.1f phi=%5.1f x=%8.3f y=%8.3f z=%8.3f %6.2f,%6.2f %6.2f,%6.2f %6.2f,%6.2f\n",fName.Data(), -// Rho(), ThetaD(),PhiD(), X(), Y(), Z(), -// ThetaXd(),PhiXd(),ThetaYd(),PhiYd(),ThetaZd(),PhiZd()); - fCenterV3.Print(opt);fPcX3.Print(opt); + ToAliInfo(fCenterX3.Print(opt)); }//Print() //__________________________________________________________________________________________________ diff --git a/RICH/AliRICHChamber.h b/RICH/AliRICHChamber.h index 08ebdb9db52..2942481fca8 100644 --- a/RICH/AliRICHChamber.h +++ b/RICH/AliRICHChamber.h @@ -20,49 +20,43 @@ public: AliRICHChamber(const AliRICHChamber &chamber):TNamed(chamber) {;} virtual ~AliRICHChamber() {;} AliRICHChamber& operator=(const AliRICHChamber&) {return *this;} - + + static Double_t AlphaFeedback(Int_t ) {return 0.030;} //determines number of feedback photons + TRotMatrix* RotMatrix() const{return fpRotMatrix;} - TString RotMatrixName() const{return "rot"+fName;} + TString RotMatrixName() const{return "rot"+fName;} TRotation Rot() const{return fRot;} - Double_t Rho() const{return fCenterV3.Mag();} //gives distance to chamber center in MRS - Double_t ThetaD() const{return fCenterV3.Theta()*TMath::RadToDeg();} //gives polar angle of chamber center in MRS - Double_t PhiD() const{return fCenterV3.Phi() *TMath::RadToDeg();} //gives azimuthal angle of chamber center in MRS + Double_t Rho() const{return fCenterX3.Mag();} //gives distance to chamber center in MRS + Double_t ThetaD() const{return fCenterX3.Theta()*TMath::RadToDeg();} //gives polar angle of chamber center in MRS + Double_t PhiD() const{return fCenterX3.Phi() *TMath::RadToDeg();} //gives azimuthal angle of chamber center in MRS Double_t ThetaXd() const{return fRot.ThetaX() *TMath::RadToDeg();} Double_t PhiXd() const{return fRot.PhiX() *TMath::RadToDeg();} Double_t ThetaYd() const{return fRot.ThetaY() *TMath::RadToDeg();} Double_t PhiYd() const{return fRot.PhiY() *TMath::RadToDeg();} Double_t ThetaZd() const{return fRot.ThetaZ() *TMath::RadToDeg();} Double_t PhiZd() const{return fRot.PhiZ() *TMath::RadToDeg();} - void RotateX(Double_t a) {fRot.RotateX(a);fCenterV3.RotateX(a);fPcX3.RotateX(a);} //rotate chamber around X by "a" degrees - void RotateY(Double_t a) {fRot.RotateY(a);fCenterV3.RotateY(a);fPcX3.RotateY(a);} //rotate chamber around Y by "a" degrees - void RotateZ(Double_t a) {fRot.RotateZ(a);fCenterV3.RotateZ(a);fPcX3.RotateZ(a);} //rotate chamber around Z by "a" degrees - Double_t X() const{return fCenterV3.X();} - Double_t Y() const{return fCenterV3.Y();} - Double_t Z() const{return fCenterV3.Z();} - TVector2 Glob2Loc(TVector3 x3)const{x3-=fPcX3;x3.Transform(fRot.Inverse());return TVector2(x3.Z()+0.5*AliRICHParam::PcSizeX(),-x3.X()+0.5*AliRICHParam::PcSizeY());}//Y and Z are misplaced????? - TVector3 Loc2Glob(TVector2 x2)const{TVector3 x3(-x2.Y()+0.5*AliRICHParam::PcSizeY(),0,x2.X()-0.5*AliRICHParam::PcSizeX());x3.Transform(fRot); x3+=fPcX3;return x3;} - - TVector2 Glob2Loc(TLorentzVector x4) const{return Glob2Loc(x4.Vect());} - - void Print(Option_t *sOption)const;//virtual - - - inline void SetToZenith(); - TRotMatrix *GetRotMatrix() const{return fpRotMatrix;} + void RotX(Double_t a) {a*=TMath::DegToRad();fRot.RotateX(a);fCenterX3.RotateX(a);fRadX3.RotateX(a);fPcX3.RotateX(a);}//degrees around X + void RotY(Double_t a) {a*=TMath::DegToRad();fRot.RotateY(a);fCenterX3.RotateY(a);fRadX3.RotateY(a);fPcX3.RotateY(a);}//degrees around Y + void RotZ(Double_t a) {a*=TMath::DegToRad();fRot.RotateZ(a);fCenterX3.RotateZ(a);fRadX3.RotateZ(a);fPcX3.RotateZ(a);}//degrees around Z + TVector3 Rad() const{return fRadX3;} //provides center of radiator position in MRS, cm + TVector3 Pc() const{return fPcX3;} //provides center of photocathond position in MRS, cm + TVector3 Center() const{return fCenterX3;} //provides center of chamber position in MRS, cm + void Print(Option_t *sOption)const; //virtual interface from TObject +//Transformations for photcathode plane + TVector2 Mrs2Pc(TVector3 x3)const{x3-=fPcX3;x3.Transform(fRot.Inverse());return TVector2(-x3.X()+0.5*AliRICHParam::PcSizeX(),x3.Y()+0.5*AliRICHParam::PcSizeY());} + TVector3 Pc2Mrs(TVector2 x2)const{TVector3 x3(-x2.X()+0.5*AliRICHParam::PcSizeX(),x2.Y()-0.5*AliRICHParam::PcSizeY(),0);x3.Transform(fRot); x3+=fPcX3;return x3;} + TVector2 Mrs2Pc(TLorentzVector x4) const{return Mrs2Pc(x4.Vect());} +//Transformations for radiator plane + TVector2 Mrs2Rad(TVector3 x3)const{x3-=fRadX3;x3.Transform(fRot.Inverse());return TVector2(-x3.X()+0.5*AliRICHParam::PcSizeX(),x3.Y()+0.5*AliRICHParam::PcSizeY());} + TVector3 Rad2Mrs(TVector2 x2)const{TVector3 x3(-x2.X()+0.5*AliRICHParam::PcSizeX(),x2.Y()-0.5*AliRICHParam::PcSizeY(),0);x3.Transform(fRot); x3+=fRadX3;return x3;} + TVector3 PMrs2Loc(TVector3 p3)const{TVector3 ploc=Rot().Invert()*p3;ploc.SetXYZ(-ploc.Px(),ploc.Py(),ploc.Pz()); return ploc;} protected: - TVector3 fCenterV3; //chamber center position in MRS (cm) + TVector3 fCenterX3; //chamber center position in MRS (cm) + TVector3 fRadX3; //radiator entrance center position in MRS (cm) TVector3 fPcX3; //PC center position in MRS (cm) TRotation fRot; //chamber rotation in MRS TRotMatrix *fpRotMatrix; //rotation matrix of the chamber with respect to MRS - ClassDef(AliRICHChamber,6) //single RICH chamber description + ClassDef(AliRICHChamber,7) //single RICH chamber description };//class AliRICHChamber -//__________________________________________________________________________________________________ -void AliRICHChamber::SetToZenith() -{ -//Put the chamber to zenith. Position of PC is shifted in X-Z plane since the origin of chamber local system is in -//left hand down coner. - fCenterV3.SetXYZ(0,AliRICHParam::Offset()-AliRICHParam::GapThickness()/2 ,0); - fPcX3.SetXYZ(0,AliRICHParam::Offset()-AliRICHParam::GapThickness()/2+5.276+0.25,0); -} -//__________________________________________________________________________________________________ + #endif //AliRICHChamber_h diff --git a/RICH/AliRICHClusterFinder.cxx b/RICH/AliRICHClusterFinder.cxx index c12ff0e9ad0..e22e715c50a 100644 --- a/RICH/AliRICHClusterFinder.cxx +++ b/RICH/AliRICHClusterFinder.cxx @@ -29,33 +29,33 @@ ClassImp(AliRICHClusterFinder) //__________________________________________________________________________________________________ AliRICHClusterFinder::AliRICHClusterFinder(AliRICH *pRICH) {//main ctor - fRICH = pRICH;//this goes first as GetDebug depends on it. - if(GetDebug())Info("main ctor","Start."); + fRICH = pRICH; + AliDebug(1,"main ctor Start."); fDigitMap = 0; fRawCluster.Reset(); fResolvedCluster.Reset(); - if(GetDebug())Info("main ctor","Stop."); + AliDebug(1,"main ctor Stop."); }//main ctor //__________________________________________________________________________________________________ void AliRICHClusterFinder::Exec() { //Main method of cluster finder. Loops on events and chambers, everything else is done in FindClusters() - if(GetDebug()) Info("Exec","Start."); + AliDebug(1,"Exec Start."); R()->GetLoader() ->LoadDigits(); - R()->GetLoader()->GetRunLoader()->LoadHeader(); - R()->GetLoader()->GetRunLoader()->LoadKinematics(); +// R()->GetLoader()->GetRunLoader()->LoadHeader(); + R()->GetLoader()->GetRunLoader()->LoadKinematics(); //header is already loaded for(Int_t iEventN=0;iEventNGetEventsPerRun();iEventN++){//events loop - if(GetDebug()) Info("Exec","Processing event %i...",iEventN); + AliDebug(1,Form("Processing event %i...",iEventN)); R()->GetLoader()->GetRunLoader()->GetEvent(iEventN); R()->GetLoader()->MakeTree("R"); R()->MakeBranch("R"); R()->ResetDigits(); R()->ResetClusters(); R()->GetLoader()->TreeD()->GetEntry(0); - for(Int_t iChamber=1;iChamber<=kNCH;iChamber++){//chambers loop + for(Int_t iChamber=1;iChamber<=kNchambers;iChamber++){//chambers loop FindClusters(iChamber); }//chambers loop R()->GetLoader()->TreeR()->Fill(); R()->GetLoader()->WriteRecPoints("OVERWRITE");//write out clusters for current event @@ -65,17 +65,17 @@ void AliRICHClusterFinder::Exec() R()->ResetClusters(); R()->GetLoader() ->UnloadDigits(); R()->GetLoader() ->UnloadRecPoints(); - R()->GetLoader()->GetRunLoader()->UnloadHeader(); +// R()->GetLoader()->GetRunLoader()->UnloadHeader(); R()->GetLoader()->GetRunLoader()->UnloadKinematics(); - if(GetDebug()) Info("Exec","Stop."); + AliDebug(1,"Stop."); }//Exec() //__________________________________________________________________________________________________ void AliRICHClusterFinder::FindClusters(Int_t iChamber) { //Loops on digits for a given chamber, forms raw clusters, then tries to resolve them if requested Int_t iNdigits=R()->Digits(iChamber)->GetEntriesFast(); - if(GetDebug())Info("FindClusters","Start for chamber %i with %i digits.",iChamber,iNdigits); + AliDebug(1,Form("Start for chamber %i with %i digits.",iChamber,iNdigits)); if(iNdigits==0)return;//no digits for a given chamber, nothing to do @@ -87,9 +87,9 @@ void AliRICHClusterFinder::FindClusters(Int_t iChamber) if(fDigitMap->TestHit(i,j)==kUsed) continue;//this digit is already taken, go after next digit FormRawCluster(i,j);//form raw cluster starting from (i,j) pad - if(GetDebug()) {Info("FindClusters","After FormRawCluster:");fRawCluster.Print();} + AliDebug(1,"After FormRawCluster:");ToAliDebug(1,fRawCluster.Print()); FindLocalMaxima(); //find number of local maxima and initial center of gravity - if(GetDebug()) {Info("FindClusters","After FindLocalMaxima:");fRawCluster.Print();} + AliDebug(1,"After FindLocalMaxima:");ToAliDebug(1,fRawCluster.Print()); if(AliRICHParam::IsResolveClusters()&&fRawCluster.Size()>1&&fRawCluster.Size()<6){ FitCoG(); //serialization of resolved clusters will happen inside @@ -101,13 +101,13 @@ void AliRICHClusterFinder::FindClusters(Int_t iChamber) delete fDigitMap; - if(GetDebug())Info("FindClusters","Stop."); + AliDebug(1,"Stop."); }//FindClusters() //__________________________________________________________________________________________________ void AliRICHClusterFinder::FindClusterContribs(AliRICHcluster *pCluster) { //Finds cerenkov-feedback-mip mixture for a given cluster - if(GetDebug()) {Info("FindClusterContribs","Start.");pCluster->Print();} + AliDebug(1,"Start.");ToAliDebug(1,pCluster->Print()); TObjArray *pDigits = pCluster->Digits(); if(!pDigits) return; //?????????? @@ -124,13 +124,13 @@ void AliRICHClusterFinder::FindClusterContribs(AliRICHcluster *pCluster) }//loop on digits of a given cluster TMath::Sort(contribs.GetSize(),contribs.GetArray(),pindex); for(Int_t iDigN=0;iDigN<3*pCluster->Size()-1;iDigN++) {//loop on digits to sort tids - if(GetDebug()) Info("FindClusterContribs","%4i for digit n. %4i",contribs[pindex[iDigN]],iDigN); + AliDebug(1,Form("%4i for digit n. %4i",contribs[pindex[iDigN]],iDigN)); if(contribs[pindex[iDigN]]!=contribs[pindex[iDigN+1]]) { TParticle* particle = R()->GetLoader()->GetRunLoader()->Stack()->Particle(contribs[pindex[iDigN]]); Int_t code = particle->GetPdgCode(); Double_t charge = 0; if(particle->GetPDG()) charge=particle->GetPDG()->Charge(); - if(GetDebug()) Info("FindClusterContribs"," charge of particle %f",charge); + AliDebug(1,Form(" charge of particle %f",charge)); if(code==50000050) iNckovs++; if(code==50000051) iNfeeds++; @@ -144,7 +144,7 @@ void AliRICHClusterFinder::FindClusterContribs(AliRICHcluster *pCluster) Int_t code = particle->GetPdgCode(); Double_t charge = 0; if(particle->GetPDG()) charge=particle->GetPDG()->Charge(); - if(GetDebug()) Info("FindClusterContribs"," charge of particle %f",charge); + AliDebug(1,Form(" charge of particle %f",charge)); if(code==50000050) iNckovs++; if(code==50000051) iNfeeds++; if(charge!=0) iNmips++; @@ -153,13 +153,14 @@ void AliRICHClusterFinder::FindClusterContribs(AliRICHcluster *pCluster) pCluster->CFM(iNckovs,iNfeeds,iNmips); // delete [] pindex; - if(GetDebug()) {pCluster->Print();Info("FindClusterContribs","Stop.");} + ToAliDebug(1,pCluster->Print()); + AliDebug(1,"Stop."); }//FindClusterContribs() //__________________________________________________________________________________________________ void AliRICHClusterFinder::FormRawCluster(Int_t i, Int_t j) { //Builds the raw cluster (before deconvolution). Starts from the first pad (i,j) then calls itself recursevly for all neighbours. - if(GetDebug()) Info("FormRawCluster","Start with digit(%i,%i) Q=%f",i,j,((AliRICHdigit*)fDigitMap->GetHit(i,j))->Q()); + AliDebug(1,Form("Start with digit(%i,%i) Q=%f",i,j,((AliRICHdigit*)fDigitMap->GetHit(i,j))->Q())); fRawCluster.AddDigit((AliRICHdigit*) fDigitMap->GetHit(i,j));//take this pad in cluster fDigitMap->FlagHit(i,j);//flag this pad as taken @@ -204,30 +205,30 @@ void AliRICHClusterFinder::FindLocalMaxima() void AliRICHClusterFinder::WriteRawCluster() { //Add the current raw cluster to the list of clusters - if(GetDebug()) Info("WriteRawCluster","Start."); + AliDebug(1,"Start."); FindClusterContribs(&fRawCluster); R()->AddCluster(fRawCluster); - if(GetDebug()){fRawCluster.Print(); Info("WriteRawCluster","Stop.");} + ToAliDebug(1,fRawCluster.Print()); AliDebug(1,"Stop."); }//WriteRawCluster() //__________________________________________________________________________________________________ void AliRICHClusterFinder::WriteResolvedCluster() { //Add the current resolved cluster to the list of clusters - if(GetDebug()) Info("WriteResolvedCluster","Start."); + AliDebug(1,"Start."); FindClusterContribs(&fResolvedCluster); R()->AddCluster(fResolvedCluster); - if(GetDebug()){fResolvedCluster.Print(); Info("WriteResolvedCluster","Stop.");} + ToAliDebug(1,fResolvedCluster.Print()); AliDebug(1,"Stop."); }//WriteResolvedCluster() //__________________________________________________________________________________________________ void AliRICHClusterFinder::FitCoG() { //Fits cluster of size by the corresponding number of Mathieson shapes. //This methode is only invoked in case everything is ok to start deconvolution - if(GetDebug()) {Info("FitCoG","Start with:"); fRawCluster.Print();} + AliDebug(1,"Start with:"); ToAliDebug(1,fRawCluster.Print()); Double_t arglist; Int_t ierflag = 0; @@ -293,7 +294,7 @@ void AliRICHClusterFinder::FitCoG() fResolvedCluster.Fill(&fRawCluster,xCoG[i],yCoG[i],qfracCoG[i],fLocalC[i]); WriteResolvedCluster(); } - if(GetDebug()) Info("CoG","Stop."); + AliDebug(1,"Stop."); }//FitCoG() //__________________________________________________________________________________________________ void RICHMinMathieson(Int_t &npar, Double_t *, Double_t &chi2, Double_t *par, Int_t ) diff --git a/RICH/AliRICHDisplFast.cxx b/RICH/AliRICHDisplFast.cxx index 656872f5339..54a253f616c 100644 --- a/RICH/AliRICHDisplFast.cxx +++ b/RICH/AliRICHDisplFast.cxx @@ -80,7 +80,7 @@ void AliRICHDisplFast::Exec(Option_t *) AliRICHhit *pHit = (AliRICHhit*)Hits[i].At(j); if(pHit->C()==iChamber){ TVector3 hitGlobX3= pHit->OutX3(); - TVector2 hitLocX2 = pRich->C(iChamber)->Glob2Loc(hitGlobX3); + TVector2 hitLocX2 = pRich->C(iChamber)->Mrs2Pc(hitGlobX3); pHitsH2->Fill(hitLocX2.X(),hitLocX2.Y(),200); }//if }//hits loop @@ -95,7 +95,7 @@ void AliRICHDisplFast::Exec(Option_t *) if(!isClusters) {l.SetTextColor(kBlue) ;l.DrawLatex(0.8,0.01,"No CLUSTERS");} Display->Update(); Display->Modified(); - getchar(); + gPad->WaitPrimitive(); //deals with digits if(isDigits){ pRich->GetLoader()->TreeD()->GetEntry(0); @@ -108,7 +108,7 @@ void AliRICHDisplFast::Exec(Option_t *) pDigitsH2->Draw("same"); Display->Update(); Display->Modified(); - getchar(); + gPad->WaitPrimitive(); }//if(isDigits) //deals with clusters if(isClusters){ @@ -121,7 +121,7 @@ void AliRICHDisplFast::Exec(Option_t *) pClustersH2->Draw("same"); Display->Update(); Display->Modified(); - getchar(); + gPad->WaitPrimitive(); }//if(isClusters) }//chambers loop delete [] Hits; diff --git a/RICH/AliRICHHelix.cxx b/RICH/AliRICHHelix.cxx new file mode 100644 index 00000000000..c5ad6d6501a --- /dev/null +++ b/RICH/AliRICHHelix.cxx @@ -0,0 +1,12 @@ +#include "AliRICHHelix.h" +#include + +ClassImp(AliRICHHelix) + +void AliRICHHelix::Print(Option_t *opt) const +{ +// Debug printout + fX0.Print(opt);fP0.Print(opt); + AliInfo("Point of interest:"); + fX.Print(opt); fP.Print(opt); +}//Print() diff --git a/RICH/AliRICHHelix.h b/RICH/AliRICHHelix.h new file mode 100644 index 00000000000..70cee594806 --- /dev/null +++ b/RICH/AliRICHHelix.h @@ -0,0 +1,99 @@ +#ifndef AliRICHHelix_h +#define AliRICHHelix_h + +#include +#include +#include "AliRICHParam.h" +#include "AliRICHChamber.h" + +class AliRICHHelix: public TObject +{ +public: + AliRICHHelix() {;} + AliRICHHelix(TVector3 x0,TVector3 p0,Int_t q,Double_t b=0.4) {fX0=x0;fP0=p0;fQ=q;fBz=b;} //takes position and momentum at parametrised point + virtual ~AliRICHHelix() {;} + inline void Propagate(Double_t lenght); + Bool_t Intersection(TVector3 plane) {return Intersection(plane,plane.Unit());} // special plane given by point only + inline Bool_t Intersection(TVector3 planePoint,TVector3 planeNorm); //intersection with plane given by point and normal vector + inline Int_t RichIntersect(AliRICHParam *pParam); + TVector3 X()const {return fX;} + TVector3 P()const {return fP;} + TVector3 X0()const {return fX0;} + TVector3 P0()const {return fP0;} + TVector2 PosPc() const {return fPosPc;} //returns position of intersection with PC + TVector2 PosRad()const {return fPosRad;} //returns position of intersection with radiator + TVector3 Ploc()const {return fPloc;} //returns momentum at the position of intersection with radiator + + void Print(Option_t *sOption)const; //virtual interface from TObject +protected: + TVector3 fX0; //helix position in parametrised point, cm in MRS + TVector3 fP0; //helix momentum in parametrised point, GeV/c in MRS + TVector3 fX; //helix position in point of interest, cm in MRS + TVector3 fP; //helix momentum in point of interest, GeV/c in MRS + Double_t fLength; //helix length in point of interest + Int_t fQ; //sign of track charge (value not provided by current ESD) + Double_t fBz; //magnetic field along z value in Tesla under assumption of uniformity + TVector2 fPosPc; //position on PC in local system + TVector2 fPosRad; //position on radiator in local system + TVector3 fPloc; //momentum in local system + ClassDef(AliRICHHelix,0) //General helix +};//class AliRICHHelix +//__________________________________________________________________________________________________ +void AliRICHHelix::Propagate(Double_t length) +{ +// Propogates the helix to the position of interest defined by helix length s +// Assumes uniform magnetic field along z direction. + const Double_t c = 0.00299792458;//this value provides that coordinates are in cm momentum in GeV/c + Double_t a = -c*fBz*fQ; + + Double_t rho = a/fP0.Mag(); + fX.SetX( fX0.X()+fP0.X()*TMath::Sin(rho*length)/a-fP0.Y()*(1-TMath::Cos(rho*length))/a ); + fX.SetY( fX0.Y()+fP0.Y()*TMath::Sin(rho*length)/a+fP0.X()*(1-TMath::Cos(rho*length))/a ); + fX.SetZ( fX0.Z()+fP0.Z()*length/fP0.Mag() ); + fP.SetX( fP0.X()*TMath::Cos(rho*length)-fP0.Y()*TMath::Sin(rho*length) ); + fP.SetY( fP0.Y()*TMath::Cos(rho*length)+fP0.X()*TMath::Sin(rho*length) ); + fP.SetZ( fP0.Z() ); + fLength=length; +} +//__________________________________________________________________________________________________ +Bool_t AliRICHHelix::Intersection(TVector3 planePoint,TVector3 planeNorm) +{ +// Finds point of intersection (if exists) of the helix to the plane given by point and normal vector. +// Returns kTrue if helix intersects the plane, kFALSE otherwise. +// Stores result in current helix fields fX and fP. + + Double_t s=(planePoint-fX0)*planeNorm,dist=99999,distPrev=dist; + + while(TMath::Abs(dist)>0.0001){ + Propagate(s); //calculates helix at the distance s from x0 ALONG the helix + dist=(fX-planePoint)*planeNorm; //distance between current helix position and plane + if(TMath::Abs(dist) > TMath::Abs(distPrev)) { return kFALSE;} + distPrev=dist; + s-=dist; + } + return kTRUE; +} +//__________________________________________________________________________________________________ +Int_t AliRICHHelix::RichIntersect(AliRICHParam *pParam) +{ +// Searchs for intersection of this helix with all RICH chambers, returns chamber number or 0 if no intersection +// On exit fPosRad contain position of intersection in Local System with radiator +// fPosPc contains the same for photocathode + for(Int_t iChamberN=1;iChamberN<=kNchambers;iChamberN++){//chamber loop + if(Intersection(pParam->C(iChamberN)->Rad())){//there is intersection with radiator plane + fPosRad=pParam->C(iChamberN)->Mrs2Rad(fX);//position on radiator plane + if(pParam->IsAccepted(fPosRad)){//intersection within radiator (even if in dead zone) + if(Intersection(pParam->C(iChamberN)->Pc())){//there is intersection with photocathode + fPosPc=pParam->C(iChamberN)->Mrs2Pc(fX);//position on photcathode plane + if(pParam->IsAccepted(fPosPc)){//intersection within pc (even if in dead zone) + fPloc=pParam->C(iChamberN)->PMrs2Loc(fP); + return iChamberN; + }//if inside PC + }//if for PC + }//if inside radiator + }//if for radiator + }//chamber loop + return 0; +} +//__________________________________________________________________________________________________ +#endif diff --git a/RICH/AliRICHParam.cxx b/RICH/AliRICHParam.cxx index ad1679fb9c8..02712443f0b 100644 --- a/RICH/AliRICHParam.cxx +++ b/RICH/AliRICHParam.cxx @@ -28,15 +28,20 @@ Float_t AliRICHParam::fgSigmaThSpread =0.035; // //__________________________________________________________________________________________________ void AliRICHParam::Print(Option_t*) { - ::Info("","Pads in chamber (%3i,%3i) in sector (%2i,%2i)",NpadsX(),NpadsY(),NpadsXsec(),NpadsYsec()); - fpChambers->Print(); -} + AliInfo(Form("Pads in chamber (%3i,%3i) in sector (%2i,%2i)",NpadsX(),NpadsY(),NpadsXsec(),NpadsYsec())); + ToAliInfo(fpChambers->Print()); +}//Print() //__________________________________________________________________________________________________ void AliRICHParam::CreateChambers() { //Create all RICH Chambers on each call. Previous chambers deleted. if(fpChambers) delete fpChambers; - fpChambers=new TObjArray(kNchambers); + if(IsRadioSrc()){ + fpChambers=new TObjArray(1);//test beam configuration 1 chamber + fpChambers->AddAt(new AliRICHChamber(0),0); + }else{ + fpChambers=new TObjArray(kNchambers);//normal configuration 7 chambers + for(int iChamberN=0;iChamberNAddAt(new AliRICHChamber(iChamberN+1),iChamberN); + } fpChambers->SetOwner(); - for(int iChamberN=0;iChamberNAddAt(new AliRICHChamber(iChamberN+1),iChamberN); -}//void AliRICH::CreateChambers() +}//CreateChambers() diff --git a/RICH/AliRICHParam.h b/RICH/AliRICHParam.h index 4557e9e5c14..2a664c17a5b 100644 --- a/RICH/AliRICHParam.h +++ b/RICH/AliRICHParam.h @@ -5,19 +5,24 @@ #include #include #include +#include #include #include #include #include +#include +#include +#include +#include +#include + -static const int kNCH=7; //number of RICH chambers ??? static const int kNchambers=7; //number of RICH chambers static const int kNpadsX = 160; //number of pads along X in single chamber static const int kNpadsY = 144; //number of pads along Y in single chamber static const int kBad=-101; //useful static const to mark initial (uninitalised) values static const int kNsectors=6; //number of sectors per chamber -static const int kadc_satm = 4096; //dynamic range (10 bits) static const int kCerenkov=50000050; //??? go to something more general like TPDGCode static const int kFeedback=50000051; //??? go to something more general like TPDGCode @@ -29,7 +34,7 @@ public: AliRICHParam():TObject(),fpChambers(0) {CreateChambers();} virtual ~AliRICHParam() {delete fpChambers;} void CreateChambers(); - AliRICHChamber* C(Int_t i) {return (AliRICHChamber*)fpChambers->UncheckedAt(i-1);} //returns pointer to chamber i + AliRICHChamber* C(Int_t i) {return (AliRICHChamber*)fpChambers->UncheckedAt(i-1);} //returns pointer to chamber i static Int_t NpadsX() {return kNpadsX;} //pads along X in chamber static Int_t NpadsY() {return kNpadsY;} //pads along Y in chamber static Int_t NpadsXsec() {return NpadsX()/2;} //pads along X in sector @@ -37,41 +42,26 @@ public: static Double_t DeadZone() {return 2.6;} //dead zone size in cm static Double_t PadSizeX() {return 0.8;} //pad size x in cm static Double_t PadSizeY() {return 0.84;} //pad size y in cm + static Double_t SectorSizeX() {return NpadsX()*PadSizeX()/2;} //sector size x in cm static Double_t SectorSizeY() {return NpadsY()*PadSizeY()/3;} //sector size y in cm - static Double_t PcSizeX() {return NpadsX()*PadSizeX()+DeadZone();} //photocathode size x in cm - static Double_t PcSizeY() {return NpadsY()*PadSizeY()+2*DeadZone();} //photocathode size y in cm - static Double_t SizeX() {return 132.6;} - static Double_t SizeY() {return 26;} - static Double_t SizeZ() {return 136.7;} - static Double_t Offset() {return 490+1.267;} //distance from IP to center of chamber in cm - static Double_t AngleYZ() {return 19.5*TMath::DegToRad();} //angle between chambers in YZ plane, rad - static Double_t AngleXY() {return 20*TMath::DegToRad();} //angle between chambers in XY plane, rad - static Double_t AngleRot() {return fgAngleRot*TMath::DegToRad();} //RICH rotation around Z, rad - static Double_t FreonThickness() {return 1.5;} - static Double_t QuartzThickness() {return 0.5;} - - static Double_t GapProx() {return 8.0;} //cm between CsI PC and radiator quartz window - static Double_t GapColl() {return 7.0;} //cm between CsI PC and third wire grid (collection wires) - static Double_t GapAnod() {return 0.204;} //cm between CsI PC and first wire grid (anod wires) - static Double_t GapAmp() {return 0.445;} //cm between CsI PC and second wire grid (cathode wires) + static Double_t PcSizeX() {return NpadsX()*PadSizeX()+DeadZone();} //PC size x, cm + static Double_t PcSizeY() {return NpadsY()*PadSizeY()+2*DeadZone();} //PC size y, cm + + static Double_t Zfreon() {return 1.5;} //freon thinkness, cm + static Double_t Zwin() {return 0.5;} //radiator quartz window, cm + static Double_t Pc2Win() {return 8.0;} //cm between CsI PC and radiator quartz window + static Double_t Pc2Coll() {return 7.0;} //cm between CsI PC and third wire grid (collection wires) + static Double_t Pc2Anod() {return 0.204;} //cm between CsI PC and first wire grid (anod wires) + static Double_t Pc2Cath() {return 0.445;} //cm between CsI PC and second wire grid (cathode wires) + static Double_t Freon2Pc() {return Zfreon()+Zwin()+Pc2Win();} //cm between CsI PC and entrance to freon static Double_t PitchAnod() {return PadSizeY()/2;} //cm between anode wires static Double_t PitchCath() {return PadSizeY()/4;} //cm between cathode wires - static Double_t PitchColl() {return 0.5;} //cm between collect wires + static Double_t PitchColl() {return 0.5;} //cm between collection wires - static Double_t GapThickness() {return 8.0;} - static Double_t RadiatorToPads() {return FreonThickness()+QuartzThickness()+GapThickness();} - static Double_t AnodeCathodeGap() {return 0.2;} //between CsI PC and first wire grid - static Double_t QuartzLength() {return 133;} - static Double_t QuartzWidth() {return 127.9;} - static Double_t OuterFreonLength() {return 133;} - static Double_t OuterFreonWidth() {return 41.3;} - static Double_t InnerFreonLength() {return 133;} - static Double_t InnerFreonWidth() {return 41.3;} static Double_t IonisationPotential() {return 26.0e-9;} //for CH4 in GeV taken from ???? static TVector2 MathiesonDelta() {return TVector2(5*0.18,5*0.18);} //area of 5 sigmas of Mathieson distribution (cm) static Int_t MaxQdc() {return 4095;} //QDC number of channels - static Double_t AlphaFeedback(Int_t ) {return 0.030;} //determines number of feedback photons static Bool_t IsResolveClusters() {return fgIsResolveClusters;} //go after resolved clusters? static Bool_t IsWireSag() {return fgIsWireSag;} //take wire sagita in account? @@ -89,6 +79,10 @@ public: static void SetWireSag(Bool_t status) {fgIsWireSag=status;} static void SetHV(Int_t sector,Int_t hv){fgHV[sector-1]=hv;} static void SetAngleRot(Double_t rot) {fgAngleRot =rot;} + static Double_t IndOfRefC6F14(Double_t eV) {return eV*0.0172+1.177;} // eV = photon energy in eV + static Double_t IndOfRefSiO2(Double_t eV) {Double_t e1=10.666,e2=18.125,f1=46.411,f2= 228.71; + return TMath::Sqrt(1.+f1/(e1*e1-TMath::Power(eV,2))+f2/(e2*e2-TMath::Power(eV,2)));}//TDR p.35 + static Double_t IndOfRefCH4() {return 1.000444;} inline static TVector Loc2Area(TVector2 x2); //return area affected by hit x2 inline static TVector Loc2Pad(TVector2 x2); //return pad containing given position @@ -107,10 +101,11 @@ public: static Float_t SigmaThMean() {return fgSigmaThMean;} //QDC electronic noise mean static Float_t SigmaThSpread() {return fgSigmaThSpread;} //QDC electronic noise width void Print(const Option_t *opt=""); //virtual - inline static void PropogateHelix(TVector3 x0,TVector3 p0,Double_t s,TVector3 *x,TVector3 *p); inline static Int_t Loc2Sec(TVector2 &x2); //return sector, x2->Sector RS - inline static Int_t Pad2Sec(TVector pad); //return sector + inline static Int_t Pad2Sec(const TVector &pad); //return sector + static Bool_t IsAccepted(const TVector2 &x2) {return ( x2.X()>=0 && x2.X()<=PcSizeX() && x2.Y()>=0 && x2.Y()<=PcSizeY() ) +? kTRUE:kFALSE;} protected: TObjArray *fpChambers; //list of chambers static Bool_t fgIsWireSag; //wire sagitta ON/OFF flag @@ -156,12 +151,12 @@ Int_t AliRICHParam::Loc2Sec(TVector2 &v2) Double_t x=v2.X(),y=v2.Y(); if (v2.X() >= x0 && v2.X() <= x1 ) {sector=1;} else if(v2.X() >= x2 && v2.X() <= x3 ) {sector=2; x=v2.X()-x2;} - else {sector=kBad; ::Error("Loc2Sec","Position %6.2f,%6.2f is out of chamber in X",v2.X(),v2.Y());return kBad;} + else {return kBad;} if (v2.Y() >= y0 && v2.Y() <= y1 ) {} //sectors 1 or 2 else if(v2.Y() >= y2 && v2.Y() <= y3 ) {sector+=2; y=v2.Y()-y2;} //sectors 3 or 4 else if(v2.Y() >= y4 && v2.Y() <= y5 ) {sector+=4; y=v2.Y()-y4;} //sectors 5 or 6 - else {sector=kBad; ::Error("Loc2Sec","Position %6.2f,%6.2f is out of chamber in Y",v2.X(),v2.Y());return kBad;} + else {return kBad;} v2.Set(x,y); return sector; }//Loc2Sec(Double_t x, Double_t y) @@ -188,27 +183,27 @@ TVector AliRICHParam::Loc2Pad(TVector2 x2) return pad; } //__________________________________________________________________________________________________ -Int_t AliRICHParam::Pad2Sec(TVector pad) +Int_t AliRICHParam::Pad2Sec(const TVector &pad) { // Determines sector containing the given pad. Int_t sector=kBad; if (pad[0] >= 1 && pad[0] <= NpadsXsec() ) {sector=1;} else if(pad[0] > NpadsXsec() && pad[0] <= NpadsX() ) {sector=2;} - else ::Error("Pad2Sec","Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]); + else AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1])); if (pad[1] >= 1 && pad[1] <= NpadsYsec() ) {} else if(pad[1] > NpadsYsec() && pad[1] <= 2*NpadsYsec() ) {sector+=2;} else if(pad[1] > 2*NpadsYsec() && pad[1] <= NpadsY() ) {sector+=4;} - else ::Error("Pad2Sec","Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]); + else AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1])); return sector; }//Pad2Sec() //__________________________________________________________________________________________________ TVector2 AliRICHParam::Pad2Loc(TVector pad) { -// Returns position of the center of the given pad in local system of the chamber +// Returns position of the center of the given pad in local system of the chamber (cm) // y ^ 5 6 -// | 3 4 chamber structure +// | 3 4 sector numbers // | 1 2 // -------> x Double_t x=kBad,y=kBad; @@ -217,7 +212,7 @@ TVector2 AliRICHParam::Pad2Loc(TVector pad) else if(pad[0] > NpadsXsec() && pad[0] <= NpadsX())//it's 2 or 4 or 6 x=(pad[0]-0.5)*PadSizeX()+DeadZone(); else - ::Error("Pad2Loc","Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]); + AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1])); if(pad[1] > 0 && pad[1] <= NpadsYsec())//it's 1 or 2 y=(pad[1]-0.5)*PadSizeY(); @@ -226,7 +221,7 @@ TVector2 AliRICHParam::Pad2Loc(TVector pad) else if(pad[1] > 2*NpadsYsec() && pad[1]<= NpadsY())//it's 5 or 6 y=(pad[1]-0.5)*PadSizeY()+2*DeadZone(); else - ::Error("Pad2Loc","Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]); + AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1])); return TVector2(x,y); } @@ -234,7 +229,8 @@ TVector2 AliRICHParam::Pad2Loc(TVector pad) Double_t AliRICHParam::GainSag(Double_t x,Int_t sector) { // Returns % of gain variation due to wire sagita. -// All curves are parametrized as per sector basis, so y must be apriory transformed to the Sector RS. +// All curves are parametrized as per sector basis, so x must be apriory transformed to the Sector RS. +// Here x is a distance along wires. x-=SectorSizeX()/2; if(x>SectorSizeX()) x-=SectorSizeX(); switch(HV(sector)){ @@ -250,8 +246,7 @@ Int_t AliRICHParam::TotQdc(TVector2 x2,Double_t eloss) { // Calculates the total charge produced by the eloss in point x2 (Chamber RS). // Returns this change parametrised in QDC channels, or 0 if the hit in the dead zone. -// eloss=0 means photons which provided for only 1 electron -// eloss > 0 for Mip +// eloss=0 means photon which produces 1 electron only eloss > 0 for Mip if(Loc2Sec(x2)==kBad) return 0; //hit in the dead zone Int_t iNelectrons=Int_t(eloss/IonisationPotential()); if(iNelectrons==0) iNelectrons=1; Double_t qdc=0; @@ -264,13 +259,13 @@ Double_t AliRICHParam::FracQdc(TVector2 x2,TVector pad) // Calculates the charge fraction induced to given pad by the hit from the given point. // Integrated Mathieson distribution is used. TVector2 center2=Pad2Loc(pad);//gives center of requested pad - Double_t normXmin=(x2.X()-center2.X()-PadSizeX()/2) /AnodeCathodeGap(); - Double_t normXmax=(x2.X()-center2.X()+PadSizeX()/2) /AnodeCathodeGap(); - Double_t normYmin=(x2.Y()-center2.Y()-PadSizeY()/2) /AnodeCathodeGap(); - Double_t normYmax=(x2.Y()-center2.Y()+PadSizeY()/2) /AnodeCathodeGap(); - - if(Loc2Sec(x2)!=Pad2Sec(pad)) return 0;//requested pad does not belong to the sector of the given hit position - else return Mathieson(normXmin, normYmin, normXmax, normYmax); + Double_t normXmin=(x2.X()-center2.X()-PadSizeX()/2) /Pc2Cath();//parametrise for Mathienson + Double_t normXmax=(x2.X()-center2.X()+PadSizeX()/2) /Pc2Cath(); + Double_t normYmin=(x2.Y()-center2.Y()-PadSizeY()/2) /Pc2Cath(); + Double_t normYmax=(x2.Y()-center2.Y()+PadSizeY()/2) /Pc2Cath(); + +//requested pad might not belong to the sector of the given hit position, hence the check: + return (Loc2Sec(x2)!=Pad2Sec(pad)) ? 0:Mathieson(normXmin, normYmin, normXmax, normYmax); } //__________________________________________________________________________________________________ Double_t AliRICHParam::Mathieson(Double_t xMin,Double_t yMin,Double_t xMax,Double_t yMax) @@ -306,21 +301,4 @@ Bool_t AliRICHParam::IsOverTh(Int_t ,TVector ,Double_t q) // Checks if the current q is over threshold and FEE will save this value to data concentrator. return (q>NsigmaTh()*(SigmaThMean()+(1.-2*gRandom->Rndm())*SigmaThSpread())); } -//__________________________________________________________________________________________________ -void AliRICHParam::PropogateHelix(TVector3 x0,TVector3 p0,Double_t s,TVector3 *x,TVector3 *p) -{ -// Propogates the helix given by (x0,p0) in MRS to the position of interest defined by helix length s - const Double_t c = 0.00299792458; - const Double_t Bz = 0.5; //field in Tesla - const Double_t q = 1; //charge in electron units - Double_t a = -c*Bz*q; - - Double_t rho = a/p0.Mag(); - p->SetX(p0.X()*TMath::Cos(rho*s)-p0.Y()*TMath::Sin(rho*s)); - p->SetY(p0.Y()*TMath::Cos(rho*s)+p0.X()*TMath::Sin(rho*s)); - p->SetZ(p0.Z()); - x->SetX(x0.X()+p0.X()*TMath::Sin(rho*s)/a-p0.Y()*(1-TMath::Cos(rho*s))/a); - x->SetY(x0.Y()+p0.Y()*TMath::Sin(rho*s)/a+p0.X()*(1-TMath::Cos(rho*s))/a); - x->SetZ(x0.Z()+p0.Z()*s/p->Mag()); -} #endif //AliRICHParam_h diff --git a/RICH/AliRICHRecon.cxx b/RICH/AliRICHRecon.cxx index 16e4bc2756b..5ceada666dc 100644 --- a/RICH/AliRICHRecon.cxx +++ b/RICH/AliRICHRecon.cxx @@ -22,371 +22,82 @@ ////////////////////////////////////////////////////////////////////////// #include -#include -#include -#include #include -#include -#include -#include -#include -#include #include #include -#include "AliLoader.h" #include "AliRICH.h" -#include "AliRICHChamber.h" #include "AliRICHParam.h" #include "AliRICHRecon.h" -#include "AliRun.h" -#include "AliStack.h" +#include "AliRICHHelix.h" +#include #define NPointsOfRing 201 -TMinuit *gAliRICHminuit ; - -void fcnrecon(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag); //__________________________________________________________________________________________________ -AliRICHRecon::AliRICHRecon(const char*name, const char*title) - :TTask(name,title) +AliRICHRecon::AliRICHRecon(AliRICHHelix *pHelix,TClonesArray *pClusters,Int_t iMipId) + :TTask("RichRec","RichPat") { - // main ctor +// main ctor + SetFreonScaleFactor(1); + fIsWEIGHT = kFALSE; fThetaBin=750; fThetaMin = 0.0; fThetaMax = 0.75; - fDTheta = 0.001; fWindowWidth = 0.060; - fNpadX = AliRICHParam::NpadsY(); - fNpadY = AliRICHParam::NpadsX(); - fPadSizeX = AliRICHParam::PadSizeY(); - fPadSizeY = AliRICHParam::PadSizeX(); - fRadiatorWidth = AliRICHParam::FreonThickness(); - fQuartzWidth = AliRICHParam::QuartzThickness(); - fGapWidth = AliRICHParam::RadiatorToPads(); - fXmin = -AliRICHParam::PcSizeY()/2.; - fXmax = AliRICHParam::PcSizeY()/2.; - fYmin = -AliRICHParam::PcSizeX()/2.; - fYmax = AliRICHParam::PcSizeX()/2.; - fRich = (AliRICH*)gAlice->GetDetector("RICH"); - fOutFile=new TFile("Anal.root","RECREATE","My Analysis histos"); - if(fIsDISPLAY) fDisplay = new TCanvas("Display","RICH Display",0,0,1200,750); - fNtuple=new TNtuple("hn","ntuple", -"Run:Trig:VertZ:Pmod:Pt:Eta:TrackTheta:TrackPhi:TrackThetaFit:TrackPhiFit:Charge::NPhotons:NPhotonsFit:InRing:MassOfParticle:HoughArea:Multiplicity:TPCLastZ"); + fDTheta = 0.001; fWindowWidth = 0.060; + fRadiatorWidth = AliRICHParam::Zfreon(); + fQuartzWidth = AliRICHParam::Zwin(); + fGapWidth = AliRICHParam::Freon2Pc() - fRadiatorWidth - fQuartzWidth; + fXmin = -AliRICHParam::PcSizeX()/2.; + fXmax = AliRICHParam::PcSizeX()/2.; + fYmin = -AliRICHParam::PcSizeY()/2.; + fYmax = AliRICHParam::PcSizeY()/2.; + SetTrackTheta(pHelix->Ploc().Theta()); + SetTrackPhi(pHelix->Ploc().Phi()); + SetMipIndex(iMipId); + SetShiftX(pHelix->PosRad().X()); + SetShiftY(pHelix->PosRad().Y()); + fpClusters = pClusters; } //__________________________________________________________________________________________________ -void AliRICHRecon::StartProcessEvent() +Double_t AliRICHRecon::ThetaCerenkov() { - //start to process for pattern recognition - - Float_t trackThetaStored = 0; - Float_t trackPhiStored = 0; - Float_t thetaCerenkovStored = 0; - Int_t houghPhotonsStored = 0; - - SetFreonScaleFactor(0.994); - - if(fIsDISPLAY) - { - DrawEvent(0); -// Waiting(); - } - - AliLoader * richLoader = Rich()->GetLoader(); - AliRunLoader * runLoader = richLoader->GetRunLoader(); - - if (richLoader->TreeH() == 0x0) richLoader->LoadHits(); - if (richLoader->TreeR() == 0x0) richLoader->LoadRecPoints(); - if (richLoader->TreeD() == 0x0) richLoader->LoadDigits(); - if (runLoader->TreeE() == 0x0) runLoader->LoadHeader(); - if (runLoader->TreeK() == 0x0) runLoader->LoadKinematics(); - - richLoader->TreeR()->GetEntry(0); - - Float_t clusX[7][500],clusY[7][500]; - Int_t clusQ[7][500],clusMul[7][500]; - Int_t nClusters[7]; - - for (Int_t ich=0;ich<7;ich++) { - nClusters[ich] = Rich()->Clusters(ich+1)->GetEntries(); - for(Int_t k=0;kClusters(ich+1)->At(k); - clusX[ich][k] = pCluster->X(); - clusY[ich][k] = pCluster->Y(); - clusQ[ich][k] = pCluster->Q(); - clusMul[ich][k] = pCluster->Size(); - // pCluster->Print(); - } - } - - Int_t nPrimaries = (Int_t)richLoader->TreeH()->GetEntries(); - - cout << " N. primaries " << nPrimaries << endl; - - for(Int_t i=0;iTreeH()->GetEntry(i); - - // Rich()->Hits()->Print(); - Int_t iPrim = 0; - - AliRICHhit* pHit=0; - - for(Int_t j=0;jHits()->GetEntries();j++) { - - pHit = (AliRICHhit*)Rich()->Hits()->At(j); - if(pHit->GetTrack() < nPrimaries) break; - iPrim++; - } - - cout << " iPrim " << iPrim << " pHit " << pHit << endl; - - if (!pHit) return; - - // pHit->Print(); - - TParticle *pParticle = runLoader->Stack()->Particle(pHit->GetTrack()); - Float_t pmod = pParticle->P(); - Float_t pt = pParticle->Pt(); - Float_t trackEta = pParticle->Eta(); - Int_t q = (Int_t)TMath::Sign(1.,pParticle->GetPDG()->Charge()); - - // pParticle->Print(); - - cout << " pmod " << pmod << " pt " << pt << " Eta " << trackEta << " charge " << q << endl; - - SetTrackMomentum(pmod); - SetTrackPt(pt); - SetTrackEta(trackEta); - SetTrackCharge(q); - - TVector3 pLocal(0,0,0);//????? - - TVector2 primLocal =Rich()->C(pHit->C())->Glob2Loc(pHit->InX3()); - - // Float_t pmodFreo = pLocal.Mag(); - Float_t trackTheta = pLocal.Theta(); - Float_t trackPhi = pLocal.Phi(); - - // cout << " trackTheta " << trackTheta << " trackPhi " << trackPhi << endl; - - SetTrackTheta(trackTheta); - SetTrackPhi(trackPhi); - - Int_t maxInd = 0; - Float_t minDist = 999.; - - // cout << " n Clusters " << nClusters[pHit->Chamber()-1] << " for chamber n. " << pHit->Chamber() << endl; - - for(Int_t j=0;jChamber()-1];j++) - { - Float_t diffx = primLocal.X() - clusX[pHit->Chamber()-1][j]; - Float_t diffy = primLocal.Y() - clusY[pHit->Chamber()-1][j]; - - - Float_t diff = sqrt(diffx*diffx + diffy*diffy); - - if(diff < minDist) - { - minDist = diff; - maxInd = j; - } - - } - - Float_t diffx = primLocal.X() - clusX[pHit->Chamber()-1][maxInd]; - Float_t diffy = primLocal.Y() - clusY[pHit->Chamber()-1][maxInd]; - - cout << " diffx " << diffx << " diffy " << diffy << endl; - - - SetMipIndex(maxInd); - SetTrackIndex(i); - - Float_t shiftX = 0;//primLocal.X()/primLocal.Z()*(fRadiatorWidth+fQuartzWidth+fGapWidth) + primLocal.X(); ????? - Float_t shiftY = 0;//primLocal.Y()/primLocal.Z()*(fRadiatorWidth+fQuartzWidth+fGapWidth) + primLocal.Y(); ????? - - SetShiftX(shiftX); - SetShiftY(shiftY); - - Float_t *pclusX = &clusX[pHit->Chamber()-1][0]; - Float_t *pclusY = &clusY[pHit->Chamber()-1][0]; - - SetCandidatePhotonX(pclusX); - SetCandidatePhotonY(pclusY); - SetCandidatePhotonsNumber(nClusters[pHit->Chamber()-1]); - - Int_t qch = clusQ[pHit->Chamber()-1][maxInd]; - - - if(minDist < 3.0 && qch > 120 && maxInd !=0) - { - - if(fIsBACKGROUND) - { - - Float_t xrndm = fXmin + (fXmax-fXmin)*gRandom->Rndm(280964); - Float_t yrndm = fYmin + (fYmax-fYmin)*gRandom->Rndm(280964); - SetShiftX(xrndm); - SetShiftY(yrndm); - - } - - PatRec(); - - trackThetaStored = GetTrackTheta(); - trackPhiStored = GetTrackPhi(); - thetaCerenkovStored = GetThetaCerenkov(); - houghPhotonsStored = GetHoughPhotons(); - - Int_t diffNPhotons = 999; - Int_t nsteps = 0; - Float_t diffTrackTheta = 999.; - Float_t diffTrackPhi = 999.; - - while(fIsMINIMIZER && GetHoughPhotons() > 2 - && diffNPhotons !=0 - && diffTrackTheta > 0.0001 - && nsteps < 10) - { - - Int_t houghPhotonsBefore = GetHoughPhotons(); - - Float_t trackThetaBefore = GetTrackTheta(); - Float_t trackPhiBefore = GetTrackPhi(); - - Minimization(); - - PatRec(); - - diffNPhotons = TMath::Abs(houghPhotonsBefore - GetHoughPhotons()); - - Float_t trackThetaAfter = GetTrackTheta(); - Float_t trackPhiAfter = GetTrackPhi(); - - diffTrackTheta = TMath::Abs(trackThetaAfter - trackThetaBefore); - diffTrackPhi = TMath::Abs(trackPhiAfter - trackPhiBefore); - - if(fDebug) - cout << " houghPhotonsBefore " << houghPhotonsBefore - << " GetHoughPhotons() " << GetHoughPhotons(); - - nsteps++; - } - - SetFittedThetaCerenkov(GetThetaCerenkov()); - SetFittedHoughPhotons(GetHoughPhotons()); - - SetTrackTheta(trackThetaStored); - SetTrackPhi(trackPhiStored); - SetThetaCerenkov(thetaCerenkovStored); - SetHoughPhotons(houghPhotonsStored); - - SetMinDist(minDist); - - FillHistograms(); - - if(fIsDISPLAY) DrawEvent(1); - - Waiting(); - - } - } - if(fIsDISPLAY) fDisplay->Print("display.ps"); -}//StartProcessEvent() -//__________________________________________________________________________________________________ -void AliRICHRecon::EndProcessEvent() -{ - // function called at the end of the event loop - - fOutFile->Write(); - fOutFile->Close(); -} -//__________________________________________________________________________________________________ -void AliRICHRecon::PatRec() -{ - //pattern recognition method based on Hough transform - - - Float_t trackTheta = GetTrackTheta(); - Float_t trackPhi = GetTrackPhi(); - Float_t pmod = GetTrackMomentum(); - Int_t iMipIndex = GetMipIndex(); +// Pattern recognition method based on Hough transform +// Return theta Cerenkov for a given track and list of clusters which are set in ctor + if(fpClusters->GetEntries()==0) return kBad;//no clusters at all for a given track Bool_t kPatRec = kFALSE; Int_t candidatePhotons = 0; - Float_t shiftX = GetShiftX(); - Float_t shiftY = GetShiftY(); - - Float_t* candidatePhotonX = GetCandidatePhotonX(); - Float_t* candidatePhotonY = GetCandidatePhotonY(); - - Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber(); - - if(fDebug) cout << " n " << candidatePhotonsNumber << endl; - SetThetaCerenkov(999.); SetHoughPhotons(0); SetHoughPhotonsNorm(0); - SetHoughRMS(999.); - - for (Int_t j=0; j < candidatePhotonsNumber; j++) - { - - SetPhotonIndex(j); - - SetPhotonFlag(0); - SetPhotonEta(-999.); - SetPhotonWeight(0.); - - if (j == iMipIndex) continue; - - - if(candidatePhotonX[j] < -64.) continue; /* avoid artificial clusters from edge uesd by Yale.... */ - - Float_t xtoentr = candidatePhotonX[j] - shiftX; - Float_t ytoentr = candidatePhotonY[j] - shiftY; - - // Float_t chargehit = fHits_charge[j]; - // if(chargehit > 150) continue; - - SetEntranceX(xtoentr); - SetEntranceY(ytoentr); - - FindPhiPoint(); - - Int_t photonStatus = PhotonInBand(); - - if(fDebug) - { - cout << " Photon n. " << j << " Status " << photonStatus << " accepted " << endl; - cout << " CandidatePhotonX[j] " << candidatePhotonX[j] << " CandidatePhotonY[j] " << candidatePhotonY[j] << endl; - } - - if(photonStatus == 0) continue; - - SetPhotonFlag(1); - - FindThetaPhotonCerenkov(); - - Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov(); - - if(fDebug) cout << " theta photon " << thetaPhotonCerenkov << endl; - SetPhotonEta(thetaPhotonCerenkov); - - candidatePhotons++; - - - } + for (Int_t j=0; j < fpClusters->GetEntries(); j++){//clusters loop + SetPhotonIndex(j); + SetPhotonFlag(0); + SetPhotonEta(-999.); + SetPhotonWeight(0.); + if (j == GetMipIndex()) continue; // do not consider MIP cluster as a candidate photon + Float_t xtoentr = ((AliRICHcluster*)fpClusters->UncheckedAt(j))->X() - GetShiftX(); + Float_t ytoentr = ((AliRICHcluster*)fpClusters->UncheckedAt(j))->Y() - GetShiftY(); + SetEntranceX(xtoentr); + SetEntranceY(ytoentr); + FindPhiPoint(); +// Int_t photonStatus = PhotonInBand(); +// if(photonStatus == 0) continue; + SetPhotonFlag(1); + FindThetaPhotonCerenkov(); + Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov(); + AliDebug(1,Form("THETA CERENKOV ---> %i",thetaPhotonCerenkov)); + SetPhotonEta(thetaPhotonCerenkov); + candidatePhotons++; + }//clusters loop if(candidatePhotons >= 1) kPatRec = kTRUE; - if(!kPatRec) return; - { - SetThetaCerenkov(999.); - SetHoughPhotons(0); - } - SetPhotonsNumber(candidatePhotonsNumber); + if(!kPatRec) return kBad; + + SetPhotonsNumber(fpClusters->GetEntries()); HoughResponse(); @@ -399,21 +110,23 @@ void AliRICHRecon::PatRec() { SetThetaCerenkov(999.); SetHoughPhotonsNorm(0.); - return; + return kBad; } if(fIsWEIGHT) FindWeightThetaCerenkov(); Float_t thetaCerenkov = GetThetaCerenkov(); + AliDebug(1,Form("Number of clusters accepted ---> %i",nPhotonHough)); + SetThetaOfRing(thetaCerenkov); - FindAreaAndPortionOfRing(); +// FindAreaAndPortionOfRing(); - Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing(); - SetHoughPhotonsNorm(nPhotonHoughNorm); +// Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing(); +// SetHoughPhotonsNorm(nPhotonHoughNorm); // Calculate the area where the photon are accepted... - +/* Float_t thetaInternal = thetaCerenkov - 0.5*fWindowWidth; SetThetaOfRing(thetaInternal); FindAreaAndPortionOfRing(); @@ -427,63 +140,11 @@ void AliRICHRecon::PatRec() Float_t houghArea = externalArea - internalArea; SetHoughArea(houghArea); +*/ + return GetThetaCerenkov(); - if(fDebug) - { - cout << " ----- SUMMARY OF RECONSTRUCTION ----- " << endl; - cout << " Rings found " << fNrings << " with thetac " << thetaCerenkov << endl; - - - cout << " Nphotons " << GetPhotonsNumber() - << " Hough " << nPhotonHough - << " norm " << nPhotonHoughNorm << endl; - - cout << " In PatRec:p " << pmod << " theta " << trackTheta << " phi " << trackPhi << endl; - cout << " ------------------------------------- " << endl; - } - - Int_t nPhotons = GetPhotonsNumber(); - - Float_t xmean = 0.; - Float_t x2mean = 0.; - Int_t nev = 0; - - for (Int_t j=0; j < nPhotons;j++) - { - SetPhotonIndex(j); - - Float_t eta = GetPhotonEta(); - - if(eta != -999.) - { - if(GetPhotonFlag() == 2) - { - - - xmean += eta; - x2mean += eta*eta; - nev++; - } - } - } - - if(nev > 0) - { - xmean /=(Float_t)nev; - x2mean /=(Float_t)nev; - } else { - xmean = 0.; - x2mean = 0.; - } - - Float_t vRMS = sqrt(x2mean - xmean*xmean); - - SetHoughRMS(vRMS); - - if(fDebug) cout << " RMS " << vRMS << endl; - -} - +}//ThetaCerenkov() +//__________________________________________________________________________________________________ void AliRICHRecon::FindEmissionPoint() { //estimate the emission point in radiator @@ -501,6 +162,7 @@ void AliRICHRecon::FindEmissionPoint() Float_t emissionPoint = fRadiatorWidth + photonLenghtMin - photonLenghtMax; SetEmissionPoint(emissionPoint); + SetEmissionPoint(fRadiatorWidth/2); // tune the emission point } @@ -509,7 +171,6 @@ Int_t AliRICHRecon::PhotonInBand() //search band fro photon candidates // Float_t massOfParticle; - Float_t beta; Float_t nfreon; Float_t thetacer; @@ -522,26 +183,16 @@ Int_t AliRICHRecon::PhotonInBand() Float_t phpad = GetPhiPoint(); - // Float_t pmod = GetTrackMomentum(); - // Float_t trackTheta = GetTrackTheta(); - // Float_t trackPhi = GetTrackPhi(); // inner radius // SetPhotonEnergy(5.6); SetEmissionPoint(fRadiatorWidth -0.0001); - SetMassHypotesis(0.93828); - - SetBetaOfParticle(); SetFreonRefractiveIndex(); - beta = GetBetaOfParticle(); nfreon = GetFreonRefractiveIndex(); - - thetacer = Cerenkovangle(nfreon,beta); - thetacer = 0.; - if(fDebug) cout << " thetacer in photoninband min " << thetacer << endl; + AliDebug(1,Form("thetacer in photoninband min %f",thetacer)); FindThetaAtQuartz(thetacer); @@ -569,19 +220,15 @@ Int_t AliRICHRecon::PhotonInBand() SetPhotonEnergy(7.7); SetEmissionPoint(0.); // SetMassHypotesis(0.139567); - SetMassHypotesis(0.); - - SetBetaOfParticle(); SetFreonRefractiveIndex(); - beta = GetBetaOfParticle(); nfreon = GetFreonRefractiveIndex(); - thetacer = Cerenkovangle(nfreon,beta); + thetacer = Cerenkovangle(nfreon,1); // thetacer = 0.75; - if(fDebug) cout << " thetacer in photoninband max " << thetacer << endl; + AliDebug(1,Form("thetacer in photoninband max %f",thetacer)); FindThetaAtQuartz(thetacer); @@ -606,7 +253,7 @@ Int_t AliRICHRecon::PhotonInBand() Float_t padradius = sqrt(TMath::Power(xtoentr,2)+TMath::Power(ytoentr,2)); - if(fDebug) printf(" rmin %f r %f rmax %f \n",innerRadius,padradius,outerRadius); + AliDebug(1,Form("rmin %f r %f rmax %f",innerRadius,padradius,outerRadius)); if(padradius>=innerRadius && padradius<=outerRadius) return 1; return 0; @@ -628,8 +275,6 @@ void AliRICHRecon::FindThetaAtQuartz(Float_t thetaCerenkov) if(trackTheta == 0) { - if(fDebug) cout << " Theta sol unique " << thetaCerenkov << endl; - thetaAtQuartz = thetaCerenkov; SetThetaAtQuartz(thetaAtQuartz); return; @@ -649,17 +294,7 @@ void AliRICHRecon::FindThetaAtQuartz(Float_t thetaCerenkov) Double_t underSqrt = 1 + b*b - c*c; - if(fDebug) - { - cout << " trackTheta " << trackTheta << endl; - cout << " TrackPhi " << trackPhi << endl; - cout << " PhiPoint " << phiPoint << endl; - cout << " ThetaCerenkov " << thetaCerenkov << endl; - cout << " den b c " << den << " b " << b << " c " << c << endl; - } - if(underSqrt < 0) { - if(fDebug) cout << " sqrt negative !!!!" << underSqrt << endl; SetThetaAtQuartz(999.); return; } @@ -670,12 +305,11 @@ void AliRICHRecon::FindThetaAtQuartz(Float_t thetaCerenkov) Double_t thetaSol1 = 2*TMath::ATan(sol1); Double_t thetaSol2 = 2*TMath::ATan(sol2); - if(fDebug) cout << " Theta sol 1 " << thetaSol1 - << " Theta sol 2 " << thetaSol2 << endl; - if(thetaSol1>0 && thetaSol1 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol1; if(thetaSol2>0 && thetaSol2 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol2; +// AliDebug(1,Form(" Theta @ quartz window %f ",thetaAtQuartz)); + SetThetaAtQuartz(thetaAtQuartz); } @@ -692,9 +326,6 @@ void AliRICHRecon::FindThetaPhotonCerenkov() const Float_t kTollerance = 0.05; - // Float_t pmod = GetTrackMomentum(); - // Float_t trackTheta = GetTrackTheta(); - // Float_t trackPhi = GetTrackPhi(); Float_t phiPoint = GetPhiPoint(); @@ -703,9 +334,9 @@ void AliRICHRecon::FindThetaPhotonCerenkov() Float_t xPoint = GetEntranceX(); Float_t yPoint = GetEntranceY(); - Float_t distPoint = sqrt(xPoint*xPoint + yPoint*yPoint); + Float_t distPoint = TMath::Sqrt(xPoint*xPoint + yPoint*yPoint); - if(fDebug) cout << " DistPoint " << distPoint << endl; + AliDebug(1,Form(" DistPoint %f ",distPoint)); // Star minimization... @@ -756,8 +387,7 @@ void AliRICHRecon::FindThetaPhotonCerenkov() radiusMean = FromEmissionToCathode(); } - if(fDebug) cout << " r1 " << radiusMin << " rmean " - << radiusMean << " r2 " << radiusMax << endl; + AliDebug(1,Form(" r1 %f rmean %f r2 %f",radiusMin,radiusMean,radiusMax)); while (TMath::Abs(radiusMean-distPoint) > kTollerance) { @@ -784,12 +414,13 @@ void AliRICHRecon::FindThetaPhotonCerenkov() nIteration++; if(nIteration>=50) { - if(fDebug) printf(" max iterations in FindPhotonCerenkov\n"); + AliDebug(1,Form(" max iterations in FindPhotonCerenkov ",nIteration)); SetThetaPhotonCerenkov(999.); return; } } + AliDebug(1,Form(" distpoint %f radius %f ",distPoint,radiusMean)); SetThetaPhotonCerenkov(thetaCerMean); } @@ -811,9 +442,6 @@ void AliRICHRecon::FindAreaAndPortionOfRing() Float_t x0 = xemiss + shiftX; Float_t y0 = yemiss + shiftY; - // Float_t pmod = GetTrackMomentum(); - // Float_t trackTheta = GetTrackTheta(); - // Float_t trackPhi = GetTrackPhi(); SetPhotonEnergy(6.85); SetFreonRefractiveIndex(); @@ -1001,31 +629,6 @@ Int_t AliRICHRecon::CheckDetectorAcceptance() const return 999; } //__________________________________________________________________________________________________ -Float_t AliRICHRecon::PhotonPositionOnCathode() -{ - // find the photon position on the CsI - // Float_t massOfParticle; - Float_t beta; - Float_t nfreon; - - - SetPhotonEnergy(6.85); - SetEmissionPoint(fRadiatorWidth/2.); - SetMassHypotesis(0.139567); - - SetBetaOfParticle(); - SetFreonRefractiveIndex(); - - beta = GetBetaOfParticle(); - nfreon = GetFreonRefractiveIndex(); - - - Float_t radius = FromEmissionToCathode(); - if (radius == 999.) return 999.; - - return 0; -} - void AliRICHRecon::FindPhotonAnglesInDRS() { // Setup the rotation matrix of the track... @@ -1143,9 +746,9 @@ void AliRICHRecon::FindPhiPoint() Float_t argY = ytoentr - emissionPoint*tan(trackTheta)*sin(trackPhi); Float_t argX = xtoentr - emissionPoint*tan(trackTheta)*cos(trackPhi); - Float_t phipad = atan2(argY,argX); + Float_t phi = atan2(argY,argX); - SetPhiPoint(phipad); + SetPhiPoint(phi); } @@ -1406,7 +1009,7 @@ void AliRICHRecon::FlagPhotons() Int_t nPhotonHough = 0; Float_t thetaCerenkov = GetThetaCerenkov(); - if(fDebug) cout << " fThetaCerenkov " << thetaCerenkov << endl; + AliDebug(1,Form(" fThetaCerenkov ",thetaCerenkov)); Float_t thetaDist= thetaCerenkov - fThetaMin; Int_t steps = (Int_t)(thetaDist / fDTheta); @@ -1418,9 +1021,6 @@ void AliRICHRecon::FlagPhotons() tmin = tavg - 0.5*fWindowWidth; tmax = tavg + 0.5*fWindowWidth; - if(fDebug) cout << " tmin " << tmin << " tmax " << tmax << endl; - if(fDebug) cout << " thetac " << thetaCerenkov << endl; - // Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber(); Int_t nPhotons = GetPhotonsNumber(); @@ -1444,305 +1044,3 @@ void AliRICHRecon::FlagPhotons() SetHoughPhotons(nPhotonHough); } -void AliRICHRecon::DrawEvent(Int_t flag) const -{ - // draw event with rings - - flag=1; // dummy to be removed... -} - -Float_t AliRICHRecon::FindMassOfParticle() -{ - // find mass of the particle from theta cerenkov - - Float_t pmod = GetTrackMomentum(); - - SetPhotonEnergy(6.85); - SetFreonRefractiveIndex(); - - Float_t thetaCerenkov = GetThetaCerenkov(); - FindBetaFromTheta(thetaCerenkov); - - Double_t beta = (Double_t)(GetBetaOfParticle()); - Double_t den = 1. - beta*beta; - if(den<=0.) return 999.; - - Double_t gamma = 1./TMath::Sqrt(den); - - Float_t mass = pmod/(beta*(Float_t)gamma); - - return mass; -} - - -void AliRICHRecon::FillHistograms() -{ - // fill histograms.. - - Float_t fittedTrackTheta, fittedTrackPhi; - - Float_t thetaCerenkov = GetThetaCerenkov(); - if(thetaCerenkov == 999.) return; - - Float_t vertZ = GetEventVertexZ(); - - Float_t trackTheta = GetTrackTheta(); - Float_t trackPhi = GetTrackPhi(); - Float_t pmod = GetTrackMomentum(); - Float_t pt = GetTrackPt(); - Float_t trackEta = GetTrackEta(); - Int_t q = GetTrackCharge(); - Float_t tPCLastZ = GetTrackTPCLastZ(); - Float_t minDist = GetMinDist(); - - fittedTrackTheta = GetFittedTrackTheta(); - fittedTrackPhi = GetFittedTrackPhi(); - Int_t fittednPhotonHough = GetFittedHoughPhotons(); - - if(fDebug) - { - cout << " p " << pmod << " ThetaC " << thetaCerenkov - << " rings " << fNrings << endl; - } - - Int_t nPhotonHough = GetHoughPhotons(); -// Float_t nPhotonHoughNorm = GetHoughPhotonsNorm(); - Float_t inRing = GetPortionOfRing(); - - Float_t massOfParticle = FindMassOfParticle(); - - Float_t houghArea = GetHoughArea(); - Float_t multiplicity = GetEventMultiplicity(); - - - Float_t var[20]; - - var[0] = 0; - var[1] = 0; - var[2] = vertZ; - var[3] = pmod; - var[4] = pt; - var[5] = trackEta; - var[6] = trackTheta; - var[7] = trackPhi; - var[8] = fittedTrackTheta; - var[9] = fittedTrackPhi; - var[10] = q; - var[11] = thetaCerenkov; - var[12] = (Float_t)nPhotonHough; - var[13] = (Float_t)fittednPhotonHough; - var[14] = inRing; - var[15] = massOfParticle; - var[16] = houghArea; - var[17] = multiplicity; - var[18] = tPCLastZ; - var[19] = minDist; - - fNtuple->Fill(var); - - - fittedTrackTheta = GetFittedTrackTheta(); - fittedTrackPhi = GetFittedTrackPhi(); - - - - if(thetaCerenkov > 0.505 && thetaCerenkov < 0.605) { - SetPhotonEnergy(6.85); - SetFreonRefractiveIndex(); - } - - Int_t nPhotons = GetPhotonsNumber(); - - for (Int_t j=0; j < nPhotons;j++) - SetPhotonIndex(j); -}//FillHistograms() -//__________________________________________________________________________________________________ -void AliRICHRecon::Minimization() -{ - // minimization to find the best theta and phi of the track - - Double_t arglist; - Int_t ierflag = 0; - - static Double_t vstart[2]; - static Double_t lower[2], upper[2]; - static Double_t step[2]={0.001,0.001}; - - Double_t trackThetaNew,trackPhiNew; - TString chname; - Double_t eps, b1, b2; - Int_t ierflg; - - gAliRICHminuit = new TMinuit(2); - gAliRICHminuit->SetObjectFit((TObject *)this); - gAliRICHminuit->SetFCN(fcnrecon); - gAliRICHminuit->mninit(5,10,7); - - vstart[0] = (Double_t)GetTrackTheta(); - vstart[1] = (Double_t)GetTrackPhi(); - - lower[0] = vstart[0] - 0.03; - if(lower[0] < 0) lower[0] = 0.; - upper[0] = vstart[0] + 0.03; - lower[1] = vstart[1] - 0.03; - upper[1] = vstart[1] + 0.03; - - - gAliRICHminuit->mnparm(0,"theta",vstart[0],step[0],lower[0],upper[0],ierflag); - gAliRICHminuit->mnparm(1," phi ",vstart[1],step[1],lower[1],upper[1],ierflag); - - arglist = -1; - - // gAliRICHminuit->FixParameter(0); - - gAliRICHminuit->SetPrintLevel(-1); -// gAliRICHminuit->mnexcm("SET PRI",&arglist, 1, ierflag); - gAliRICHminuit->mnexcm("SET NOGR",&arglist, 1, ierflag); - gAliRICHminuit->mnexcm("SET NOW",&arglist, 1, ierflag); - arglist = 1; - gAliRICHminuit->mnexcm("SET ERR", &arglist, 1,ierflg); - arglist = -1; - - // gAliRICHminuit->mnscan(); - -// gAliRICHminuit->mnexcm("SIMPLEX",&arglist, 0, ierflag); - gAliRICHminuit->mnexcm("MIGRAD",&arglist, 0, ierflag); - gAliRICHminuit->mnexcm("EXIT" ,&arglist, 0, ierflag); - - gAliRICHminuit->mnpout(0,chname, trackThetaNew, eps , b1, b2, ierflg); - gAliRICHminuit->mnpout(1,chname, trackPhiNew, eps , b1, b2, ierflg); - - //values after the fit... - SetFittedTrackTheta((Float_t)trackThetaNew); - SetFittedTrackPhi((Float_t)trackPhiNew); - - delete gAliRICHminuit; - -} - -void AliRICHRecon::EstimationOfTheta() -{ - // theta estimate - - Int_t nPhotons = 0; - - Float_t shiftX = GetShiftX(); - Float_t shiftY = GetShiftY(); - - Float_t *candidatePhotonX = GetCandidatePhotonX(); - Float_t *candidatePhotonY = GetCandidatePhotonY(); - - Int_t nPhotonsCandidates = GetCandidatePhotonsNumber(); - - // cout << "MINIM: Nphotons " << nPhotonsCandidates << endl; - - for (Int_t j=0; j < nPhotonsCandidates; j++) - { - - SetPhotonIndex(j); - - if(!GetPhotonFlag()) continue; - - Float_t xtoentr = candidatePhotonX[j] - shiftX; - Float_t ytoentr = candidatePhotonY[j] - shiftY; - - SetEntranceX(xtoentr); - SetEntranceY(ytoentr); - - FindPhiPoint(); - - FindThetaPhotonCerenkov(); - - Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov(); - - // cout << " ACCEPTED!!! " << thetaPhotonCerenkov << endl; - - SetPhotonEta(thetaPhotonCerenkov); - - nPhotons++; - - } - - Float_t xmean = 0.; - Float_t x2mean = 0.; - Int_t nev = 0; - - for (Int_t j=0; j < nPhotonsCandidates;j++) - { - SetPhotonIndex(j); - - Float_t eta = GetPhotonEta(); - - if(eta != -999.) - { - if(GetPhotonFlag() == 2) - { - xmean += eta; - x2mean += eta*eta; - nev++; - } - } - } - - if(nev > 0) - { - xmean /=(Float_t)nev; - x2mean /=(Float_t)nev; - } else { - xmean = 0.; - x2mean = 0.; - } - - Float_t vRMS = sqrt(x2mean - xmean*xmean); - - // cout << " RMS " << vRMS; - - SetEstimationOfTheta(xmean); - SetEstimationOfThetaRMS(vRMS); -} - -void fcnrecon(Int_t& /*npar*/, Double_t* /*gin*/, Double_t &f, Double_t *par, Int_t) -{ - // function to be minimized - AliRICHRecon *gMyRecon = (AliRICHRecon*)gAliRICHminuit->GetObjectFit(); - - Float_t p0 = (Float_t)par[0]; - Float_t p1 = (Float_t)par[1]; - - gMyRecon->SetTrackTheta(p0); - gMyRecon->SetTrackPhi(p1); - - gMyRecon->EstimationOfTheta(); - Float_t vRMS = gMyRecon->GetEstimationOfThetaRMS(); - - Int_t houghPhotons = gMyRecon->GetHoughPhotons(); - - - f = (Double_t)(1000*vRMS/(Float_t)houghPhotons); - -// if(fDebug) cout << " f " << f -// << " theta " << par[0] << " phi " << par[1] -// << " HoughPhotons " << houghPhotons << endl; -// -// if(fDebug&&iflag == 3) -// { -// cout << " --- end convergence...summary --- " << endl; -// cout << " theta " << par[0] << endl; -// cout << " phi " << par[1] << endl; -// } -} - -void AliRICHRecon::Waiting() -{ - // wait, wait.... - if(!fIsDISPLAY) return; - cout << " Press any key to continue..."; - -// gSystem->ProcessEvents(); - getchar(); - - cout << endl; - - return; -} - diff --git a/RICH/AliRICHRecon.h b/RICH/AliRICHRecon.h index 03132c51be2..419ce8170a2 100644 --- a/RICH/AliRICHRecon.h +++ b/RICH/AliRICHRecon.h @@ -15,24 +15,15 @@ #include -class AliRICH; -class TFile; -class TNtuple; -class TCanvas; +class AliRICHHelix; class AliRICHRecon : public TTask { public : - AliRICHRecon(const char*, const char*); - ~AliRICHRecon(){EndProcessEvent();} + AliRICHRecon(AliRICHHelix *pHelix,TClonesArray *pClusters,Int_t iMipId); + virtual ~AliRICHRecon(){;} - AliRICH* Rich() {return fRich;} //main pointer to RICH - void StartProcessEvent(); // - void EndProcessEvent(); // - //PH void InitRecon(); // - void PatRec(); // - void Minimization(); // - void FillHistograms(); // + Double_t ThetaCerenkov(); // it returns reconstructed Theta Cerenkov void FindThetaPhotonCerenkov(); // void FindAreaAndPortionOfRing(); // void FindEmissionPoint(); // @@ -45,12 +36,8 @@ public : void FindWeightThetaCerenkov(); // void EstimationOfTheta(); // void FindIntersectionWithDetector(); // - void Waiting(); // - Float_t FindMassOfParticle(); // Float_t Cerenkovangle(Float_t n, Float_t b);// - Float_t PhotonPositionOnCathode(); // Int_t PhotonInBand(); // - void DrawEvent(Int_t flag) const; // Int_t CheckDetectorAcceptance() const; // Int_t GetFittedHoughPhotons() const{ return fFittedHoughPhotons;} // Int_t GetPhotonFlag() const{ return fPhotonFlag[fPhotonIndex];} // @@ -61,13 +48,10 @@ public : Int_t GetTrackIndex() const{ return fTrackIndex;} // Int_t GetCandidatePhotonsNumber() const{ return fCandidatePhotonsNumber;} // Int_t GetHoughPhotons() const{ return fHoughPhotons;} // - Float_t GetEventVertexZ() const{ return fEventVertZ;} // - Float_t GetEventMultiplicity() const{ return fEventMultiplicity;} // Float_t GetPhotonEnergy() const{ return fPhotonEnergy;} // Float_t GetFreonRefractiveIndex() const{ return fFreonRefractiveIndex;} // Float_t GetQuartzRefractiveIndex() const{ return fQuartzRefractiveIndex;} // Float_t GetGasRefractiveIndex() const{ return fGasRefractiveIndex;} // - Float_t GetFreonScaleFactor() const{ return fFreonScaleFactor;} // Float_t GetEmissionPoint() const{ return fLengthEmissionPoint;} // Float_t GetMassHypotesis() const{ return fMassHypotesis;} // Float_t GetBetaOfParticle() const{ return fTrackBeta;} // @@ -75,13 +59,8 @@ public : Float_t GetEntranceY() const{ return fYtoentr;} // Float_t GetThetaCerenkov() const{ return fThetaCerenkov;} // Float_t GetThetaPhotonCerenkov() const{ return fThetaPhotonCerenkov;} // - Float_t GetTrackMomentum() const{ return fTrackMomentum;} // - Float_t GetTrackEta() const{ return fTrackEta;} // Float_t GetTrackTheta() const{ return fTrackTheta;} // Float_t GetTrackPhi() const{ return fTrackPhi;} // - Float_t GetTrackPt() const{ return fTrackPt;} // - Float_t GetTrackTPCLastZ() const{ return fTrackTPCLastZ;} // - Float_t GetMinDist() const{ return fMinDist;} // Float_t GetXPointOnCathode() const{ return fPhotonLimitX;} // Float_t GetYPointOnCathode() const{ return fPhotonLimitY;} // Float_t GetThetaPhotonInDRS() const{ return fThetaPhotonInDRS;} // @@ -111,24 +90,17 @@ public : Float_t GetPhotonEta() const{ return fPhotonEta[fPhotonIndex];} // Float_t GetPhotonWeight() const{ return fPhotonWeight[fPhotonIndex];} // Float_t GetHoughRMS() const{ return fHoughRMS;} // - Float_t* GetCandidatePhotonX() const{ return fCandidatePhotonX;} // - Float_t* GetCandidatePhotonY() const{ return fCandidatePhotonY;} // - Float_t GetHoughPhotonsNorm() const{ return fHoughPhotonsNorm;} // Float_t GetFittedTrackTheta() const{ return fFittedTrackTheta;} // Float_t GetFittedTrackPhi() const{ return fFittedTrackPhi;} // Float_t GetFittedThetaCerenkov() const{ return fFittedThetaCerenkov;} // Float_t GetEstimationOfTheta() const{ return fEstimationOfTheta;} // Float_t GetEstimationOfThetaRMS() const{ return fEstimationOfThetaRMS;} // - void SetEventVertexZ(Float_t EventVertZ) { fEventVertZ = EventVertZ;} // - void SetEventMultiplicity(Float_t EventMultiplicity) { fEventMultiplicity = EventMultiplicity;} // void SetPhotonEnergy(Float_t PhotonEnergy) { fPhotonEnergy = PhotonEnergy;} // void SetFreonRefractiveIndex() {fFreonRefractiveIndex = fFreonScaleFactor*(1.177+0.0172*fPhotonEnergy);}// void SetQuartzRefractiveIndex() {fQuartzRefractiveIndex = sqrt(1+(46.411/(113.763556-TMath::Power(fPhotonEnergy,2)))+(228.71/(328.51563-TMath::Power(fPhotonEnergy,2))));}// void SetGasRefractiveIndex() { fGasRefractiveIndex = 1.;} // void SetFreonScaleFactor(Float_t FreonScaleFactor) {fFreonScaleFactor = FreonScaleFactor;} // void SetEmissionPoint(Float_t LengthEmissionPoint) { fLengthEmissionPoint = LengthEmissionPoint;} // - void SetMassHypotesis(Float_t mass) {fMassHypotesis = mass;} // - void SetBetaOfParticle() { fTrackBeta = fTrackMomentum/sqrt(TMath::Power(fTrackMomentum,2)+TMath::Power(fMassHypotesis,2));}// void SetEntranceX(Float_t Xtoentr) { fXtoentr = Xtoentr;} // void SetEntranceY(Float_t Ytoentr) { fYtoentr = Ytoentr;} // void SetThetaPhotonInTRS(Float_t Theta) {fThetaPhotonInTRS = Theta;} // @@ -149,14 +121,9 @@ public : void SetRadiusOuterRing(Float_t OuterRadius) {fOuterRadius = OuterRadius;} // void SetThetaCerenkov(Float_t ThetaCer) {fThetaCerenkov = ThetaCer;} // void SetThetaPhotonCerenkov(Float_t ThetaPhotCer) {fThetaPhotonCerenkov = ThetaPhotCer;} // - void SetTrackMomentum(Float_t TrackMomentum) {fTrackMomentum = TrackMomentum;} // - void SetTrackEta(Float_t TrackEta) {fTrackEta = TrackEta;} // void SetTrackTheta(Float_t TrackTheta) { fTrackTheta = TrackTheta;} // void SetTrackPhi(Float_t TrackPhi) { fTrackPhi = TrackPhi;} // - void SetTrackPt(Float_t TrackPt) { fTrackPt = TrackPt;} // void SetTrackCharge(Int_t TrackCharge) { fTrackCharge = TrackCharge;} // - void SetTrackTPCLastZ(Float_t TrackTPCLastZ) { fTrackTPCLastZ = TrackTPCLastZ;} // - void SetMinDist(Float_t MinDist) { fMinDist = MinDist;} // void SetShiftX(Float_t ShiftX) { fShiftX = ShiftX;} // void SetShiftY(Float_t ShiftY) { fShiftY = ShiftY;} // void SetDetectorWhereX(Float_t Xcoord) { fXcoord = Xcoord;} // @@ -175,9 +142,6 @@ public : void SetHoughRMS(Float_t HoughRMS) { fHoughRMS = HoughRMS;} // void SetMipIndex(Int_t MipIndex) { fMipIndex = MipIndex;} // void SetTrackIndex(Int_t TrackIndex) { fTrackIndex = TrackIndex;} // - void SetCandidatePhotonX(Float_t *CandidatePhotonX) { fCandidatePhotonX = CandidatePhotonX;} // - void SetCandidatePhotonY(Float_t *CandidatePhotonY) { fCandidatePhotonY = CandidatePhotonY;} // - void SetCandidatePhotonsNumber(Int_t CandidatePhotonsNumber) { fCandidatePhotonsNumber = CandidatePhotonsNumber;}// void SetHoughPhotons(Int_t HoughPhotons) { fHoughPhotons = HoughPhotons;} // void SetHoughPhotonsNorm(Float_t HoughPhotonsNorm) { fHoughPhotonsNorm = HoughPhotonsNorm;} // void SetFittedTrackTheta(Float_t FittedTrackTheta) { fFittedTrackTheta = FittedTrackTheta;} // @@ -191,7 +155,7 @@ public : Float_t FromEmissionToCathode(); // protected: - AliRICH* fRich; // main poiter to RICH + TClonesArray *fpClusters; // poiter to clusters Int_t fTrackCharge; // charge track Int_t fMipIndex; // index for Mip Int_t fTrackIndex; // index for track @@ -201,14 +165,9 @@ protected: Int_t fCandidatePhotonsNumber; // number of candidate photons Int_t fHoughPhotons; // n. photons after Hough Int_t fFittedHoughPhotons; // n. photons after Hough and after minimization - Float_t fEventVertZ; // z coord. of the primary vertex - Float_t fEventMultiplicity; // event primary multiplicity + Float_t fTrackTheta; // Theta of track at RICH Float_t fTrackPhi; // Phi of track at RICH - Float_t fTrackMomentum; // track momentum - Float_t fTrackEta; // track pseudorapidity - Float_t fTrackPt; // pt of track - Float_t fTrackTPCLastZ; // z last point of the TPC Float_t fMinDist; // min distance between extrapolated track and MIP Float_t fTrackBeta; // beta of the track Float_t fXtoentr; // X entrance to RICH @@ -236,7 +195,6 @@ protected: Float_t fPhotonLimitX; // X phys limit for photon Float_t fPhotonLimitY; // Y phys limit for photon Float_t fDistanceFromCluster; // distance from cluster - Float_t fMassHypotesis; // reconstructed mass Float_t fCerenkovAnglePad; // cherenkov angle of pad Float_t fThetaPhotonCerenkov; // theta cerenkov for photon Float_t fShiftX; // x shift to entrance in radiator @@ -245,6 +203,7 @@ protected: Float_t fYcoord; // .. Float_t fIntersectionX; // .. Float_t fIntersectionY; // .. + Float_t fMassHypotesis; // Float_t fThetaOfRing; // theta of ring Float_t fAreaOfRing; // area of the ring Float_t fPortionOfRing; // fraction of the accepted ring @@ -263,19 +222,9 @@ protected: Int_t fThetaBin; // bin in theta Float_t fThetaMin,fThetaMax; // min max Float_t fXmin,fXmax,fYmin,fYmax; // xy min max - TFile *fOutFile; // histogram output - TNtuple *fNtuple; // output ntuple - TCanvas *fDisplay; // for display Int_t fNrings; //current number of reconstructed rings - Bool_t fDebug; // flag for debug - Bool_t fIsDISPLAY; // flag for display Bool_t fIsWEIGHT; // flag to consider weight procedure Bool_t fIsBACKGROUND; // flag to simulate bkg - Bool_t fIsMINIMIZER; // flag to start with (theta,phi) minimization - Int_t fNpadX; // npadx - Int_t fNpadY; // npady - Float_t fPadSizeX; // sizepad x - Float_t fPadSizeY; // sizepad y Float_t fRadiatorWidth; // radiator width Float_t fQuartzWidth; // quartz width Float_t fGapWidth; // gap width diff --git a/RICH/AliRICHReconstructor.cxx b/RICH/AliRICHReconstructor.cxx index 581478b2988..47424a94887 100644 --- a/RICH/AliRICHReconstructor.cxx +++ b/RICH/AliRICHReconstructor.cxx @@ -13,7 +13,6 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* $Id$ */ /////////////////////////////////////////////////////////////////////////////// // // @@ -23,48 +22,70 @@ #include "AliRICHReconstructor.h" -#include "AliRunLoader.h" -#include "AliRun.h" #include "AliRICHClusterFinder.h" - +#include "AliRICHHelix.h" +#include +#include +#include ClassImp(AliRICHReconstructor) - -//_____________________________________________________________________________ -void AliRICHReconstructor::Reconstruct(AliRunLoader* runLoader) const +//__________________________________________________________________________________________________ +void AliRICHReconstructor::Reconstruct(AliRunLoader* pAL) const { -// reconstruct clusters - - AliRICH* rich = GetRICH(runLoader); - if (!rich) return; - AliRICHClusterFinder clusterer(rich); - clusterer.Exec(); +//Finds clusters out of digits + AliDebug(1,"Start cluster finder.");AliRICHClusterFinder clus(GetRICH(pAL)); clus.Exec(); } - -//_____________________________________________________________________________ -void AliRICHReconstructor::FillESD(AliRunLoader* /*runLoader*/, - AliESD* /*esd*/) const +//__________________________________________________________________________________________________ +void AliRICHReconstructor::FillESD(AliRunLoader* /*pAL*/, AliESD* pESD) const { -// nothing to be done - -} - +//This methode fills AliESDtrack with information from RICH + AliDebug(1,Form("Start with %i tracks",pESD->GetNumberOfTracks())); +/* const Double_t masses[5]={0.000511,0.105658,0.139567,0.493677,0.93828};//electron,muon,pion,kaon,proton + const Double_t refIndex = 1.29052; -//_____________________________________________________________________________ -AliRICH* AliRICHReconstructor::GetRICH(AliRunLoader* runLoader) const + Double_t thetaTh[5]; + Double_t sinThetaThNorm; + Double_t sigmaThetaTh[5]; + Double_t height[5]; + Double_t totalHeight=0; + Double_t x[3],p[3]; //tmp storage for track parameters + + for(Int_t iTrackN=0;iTrackNGetNumberOfTracks();iTrackN++){//ESD tracks loop + AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track +// if((pTrack->GetStatus()&AliESDtrack::kTOFout)==0) continue; //ignore tracks not recontructed by TOF + pTrack->GetXYZ(x); pTrack->GetPxPyPz(p); Double_t pmod=pTrack->GetP();//get running track parameters + continue; +// AliRICHHelix helix(x[0],x[1],x[2],p[0],p[1],p[2]); //construct helix from running track parameters + +// TVector rad(1,5); TVector pc(1,5); +// helix.RichIntersect(GetRICH(pAL)->P(),rad,pc); //returns cross point of track with RICH PC in LRS + + for(Int_t iPart=4;iPart>=0;iPart--){ + Double_t cosThetaTh = TMath::Sqrt(masses[iPart]*masses[iPart]+pmod*pmod)/(refIndex*pmod); + if(cosThetaTh>=1) {break;} + thetaTh[iPart] = TMath::ACos(cosThetaTh); + sinThetaThNorm = TMath::Sin(thetaTh[iPart])/TMath::Sqrt(1-1/(refIndex*refIndex)); + sigmaThetaTh[iPart] = (0.014*(1/sinThetaThNorm-1) + 0.0043)*1.25; + height[iPart] = TMath::Gaus(thetaExp,thetaTh[iPart],sigmaThetaTh[iPart]); + totalHeight +=height[iPart]; + } + + Double_t richPID[5]; + for(Int_t iPart=0;iPart<5;iPart++) richPID[iPart] = height[iPart]/totalHeight; + pTrack->SetRICHpid(richPID); + }//ESD tracks loop + */ +}//FillESD +//__________________________________________________________________________________________________ +AliRICH* AliRICHReconstructor::GetRICH(AliRunLoader* pAL) const { // get the RICH detector - if (!runLoader->GetAliRun()) runLoader->LoadgAlice(); - if (!runLoader->GetAliRun()) { - Error("GetRICH", "couldn't get AliRun object"); - return NULL; - } - AliRICH* rich = (AliRICH*) runLoader->GetAliRun()->GetDetector("RICH"); - if (!rich) { - Error("GetRICH", "couldn't get RICH detector"); - return NULL; - } - return rich; + if (!pAL->GetAliRun()) pAL->LoadgAlice(); + if (!pAL->GetAliRun()) {AliError("couldn't get AliRun object"); return NULL; } + AliRICH* pRich = (AliRICH*) pAL->GetAliRun()->GetDetector("RICH"); + if (!pRich) {AliError("couldn't get RICH detector"); return NULL; } + return pRich; } +//__________________________________________________________________________________________________ diff --git a/RICH/AliRICHReconstructor.h b/RICH/AliRICHReconstructor.h index a17c42d82b6..6bd045d8905 100644 --- a/RICH/AliRICHReconstructor.h +++ b/RICH/AliRICHReconstructor.h @@ -1,27 +1,26 @@ -#ifndef ALIRICHRECONSTRUCTOR_H -#define ALIRICHRECONSTRUCTOR_H +#ifndef AliRICHReconstructor_h +#define AliRICHReconstructor_h /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * See cxx source for full Copyright notice */ -/* $Id$ */ - -#include "AliReconstructor.h" +#include +#include +#include "AliRICHTracker.h" class AliRICH; - -class AliRICHReconstructor: public AliReconstructor { +class AliRICHReconstructor: public AliReconstructor +{ public: - AliRICHReconstructor(): AliReconstructor() {}; - virtual ~AliRICHReconstructor() {}; - - virtual void Reconstruct(AliRunLoader* runLoader) const; - virtual void FillESD(AliRunLoader* runLoader, AliESD* esd) const; - + AliRICHReconstructor(): AliReconstructor() {AliDebug(1,"Start.");} + virtual ~AliRICHReconstructor() {AliDebug(1,"Start.");} + virtual AliTracker* CreateTracker(AliRunLoader*)const {AliDebug(1,"Start.");return new AliRICHTracker;} //virtual from AliReconstructor + virtual void Reconstruct(AliRunLoader* pAL) const; //virtual from AliReconstructor + virtual void FillESD(AliRunLoader* pAL, AliESD* pESD) const; //virtual from AliReconstructor private: - AliRICH* GetRICH(AliRunLoader* runLoader) const; + AliRICH* GetRICH(AliRunLoader* pAL) const; - ClassDef(AliRICHReconstructor, 0) // class for the RICH reconstruction + ClassDef(AliRICHReconstructor, 0) //class for the RICH reconstruction }; #endif diff --git a/RICH/AliRICHTracker.cxx b/RICH/AliRICHTracker.cxx new file mode 100644 index 00000000000..0208b581133 --- /dev/null +++ b/RICH/AliRICHTracker.cxx @@ -0,0 +1,72 @@ +#include "AliRICHTracker.h" +#include +#include +#include +#include "AliRICH.h" +#include "AliRICHHelix.h" +#include +#include "AliRICHRecon.h" +#include +#include + +ClassImp(AliRICHTracker) + + +Int_t AliRICHTracker::PropagateBack(AliESD *pESD) +{ + Int_t iNtracks=pESD->GetNumberOfTracks(); + Double_t b=GetFieldMap()->SolenoidField()/10;// magnetic field in Tesla + AliDebug(1,Form("Start with %i tracks in %f Tesla field",iNtracks,b)); + Double_t xb[3],pb[3];//tmp storage for track parameters + TVector3 x0(0,0,0); TVector3 p0(0,0,0);//tmp storage for AliRICHHelix + + AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH")); + +// pRich->GetLoader()->GetRunLoader()->LoadHeader();pRich->GetLoader()->GetRunLoader()->LoadKinematics(); +// AliStack *pStack = pRich->GetLoader()->GetRunLoader()->Stack(); +// TParticle *pParticle = pStack->Particle(0); +// p0.SetMagThetaPhi(pParticle->P(),pParticle->Theta(),pParticle->Phi()); + + for(Int_t iTrackN=0;iTrackNGetTrack(iTrackN);// get next reconstructed track +// if((pTrack->GetStatus()&AliESDtrack::kTOFout)==0) continue; //ignore tracks not recontructed by TOF + pTrack->GetXYZ(xb); + pTrack->GetPxPyPz(pb); + Int_t status=pTrack->GetStatus()&AliESDtrack::kTOFout;//get running track parameters + AliDebug(1,Form("Track %i pmod=%f mass=%f stat=%i",iTrackN,pTrack->GetP(),pTrack->GetMass(),status)); +// x.Print();p.Print(); + x0.SetXYZ(xb[0],xb[1],xb[2]); p0.SetXYZ(xb[0],xb[1],xb[2]); + AliRICHHelix helix(x0,p0,pTrack->GetSign(),b); + Int_t iChamber=helix.RichIntersect(pRich->P()); + AliDebug(1,Form("intersection with %i chamber found",iChamber)); + if(!iChamber) continue;//intersection with no chamber found + + Double_t distMip=9999; //min distance between clusters and track position on PC + Int_t iMipId=0; //index of that min distance cluster + for(Int_t iClusN=0;iClusNClusters(iChamber)->GetEntries();iClusN++){//clusters loop for intersected chamber + AliRICHcluster *pClus=(AliRICHcluster*)pRich->Clusters(iChamber)->UncheckedAt(iClusN);//get pointer to current cluster + Double_t distCurrent=pClus->DistTo(helix.PosPc());//ditance between current cluster and helix intersection with PC + if(distCurrentDistTo(helix.PosPc()))); + }////clusters loop for intersected chamber + + AliDebug(1,Form("Min distance cluster: %i dist is %f",iMipId,distMip)); + + AliRICHRecon recon(&helix,pRich->Clusters(iChamber),iMipId); + Double_t thetaCerenkov=recon.ThetaCerenkov(); //search for mean Cerenkov angle for this track + AliDebug(1,Form("FINAL Theta Cerenkov=%f",thetaCerenkov)); + pTrack->SetRICHsignal(thetaCerenkov); + + + }//ESD tracks loop + AliDebug(1,"Stop."); + return 0; +} //pure virtual from AliTracker + +Int_t AliRICHTracker::LoadClusters(TTree *pTree) +{ +// Load clusters for RICH + AliDebug(1,"Start."); pTree->GetEntry(0); AliDebug(1,"Stop."); return 0; +} diff --git a/RICH/AliRICHTracker.h b/RICH/AliRICHTracker.h new file mode 100644 index 00000000000..7fd064dd1fc --- /dev/null +++ b/RICH/AliRICHTracker.h @@ -0,0 +1,28 @@ +#ifndef AliRICHTracker_h +#define AliRICHTracker_h + +#include +#include + +class AliCluster; +class AliESD; +class TTree; + +class AliRICHTracker : public AliTracker +{ +public: + AliRICHTracker() :AliTracker() {AliDebug(1,"Start.");} + virtual ~AliRICHTracker() {AliDebug(1,"Start.");} + Int_t Clusters2Tracks(AliESD *) {AliDebug(1,"Start.");return 0;} //pure virtual from AliTracker + Int_t RefitInward(AliESD *) {AliDebug(1,"Start.");return 0;} //pure virtual from AliTracker + void UnloadClusters() {AliDebug(1,"Start.");} //pure virtual from AliTracker + AliCluster *GetCluster(Int_t )const {AliDebug(1,"Start.");return 0;} //pure virtual from AliTracker + Int_t PropagateBack(AliESD *); //pure virtual from AliTracker + Int_t LoadClusters(TTree *); //pure virtual from AliTracker + +protected: + + ClassDef(AliRICHTracker,0) +};//class AliRICHTracker + +#endif//AliRICHTracker_h diff --git a/RICH/AliRICHv0.cxx b/RICH/AliRICHv0.cxx index 9bfebf42bae..06adbc84857 100644 --- a/RICH/AliRICHv0.cxx +++ b/RICH/AliRICHv0.cxx @@ -29,25 +29,21 @@ ClassImp(AliRICHv0) void AliRICHv0::StepManager() { //This StepManager is a provision for different test-learn activities on the current MC layer + static Int_t iStepN; const char *sParticle; switch(gMC->TrackPid()){ - case kProton: - sParticle="proton";break; - case kNeutron: - sParticle="neutron";break; - case kGamma: - sParticle="gamma";break; - case kCerenkov: - sParticle="photon";break; - case kPi0: - sParticle="Pi0";break; - case kElectron: - sParticle="electron";break; - default: - sParticle="not known";break; + case kProton: sParticle="proton" ;break; + case kNeutron: sParticle="neutron" ;break; + case kGamma: sParticle="gamma" ;break; + case kCerenkov: sParticle="photon" ;break; + case kPi0: sParticle="Pi0" ;break; + case kPiPlus: sParticle="Pi+" ;break; + case kPiMinus: sParticle="Pi-" ;break; + case kElectron: sParticle="electron" ;break; + default: sParticle="not known" ;break; } - Info("","event=%i hunt=%i tid=%i pid=%i(%s) m=%f q=%3.1f dEdX=%9.3f", + Info(Form("Step %i",iStepN),"event=%i hunt=%i tid=%i pid=%i(%s) m=%f q=%3.1f dEdX=%9.3f", gMC->CurrentEvent(), fIshunt, gAlice->GetMCApp()->GetCurrentTrackNumber(), @@ -55,8 +51,8 @@ void AliRICHv0::StepManager() sParticle, gMC->TrackMass(), gMC->TrackCharge(), - gMC->Edep()); - Info("","Flags:alive(%i) disap(%i) enter(%i) exit(%i) inside(%i) out(%i) stop(%i) new(%i)", + gMC->Edep()*1e9); + Info("Flags","alive(%i) disap(%i) enter(%i) exit(%i) inside(%i) out(%i) stop(%i) new(%i)", gMC->IsTrackAlive(), gMC->IsTrackDisappeared(), gMC->IsTrackEntering(), @@ -65,12 +61,12 @@ void AliRICHv0::StepManager() gMC->IsTrackOut(), gMC->IsTrackStop(), gMC->IsNewTrack()); - Int_t copy0,copy1,copy2,copy3; + Int_t copy0=0,copy1=0,copy2=0,copy3=0; Int_t vid0=gMC->CurrentVolID(copy0); Int_t vid1=gMC->CurrentVolOffID(1,copy1); Int_t vid2=gMC->CurrentVolOffID(2,copy2); Int_t vid3=gMC->CurrentVolOffID(3,copy3); - Info("","vid0=%i(%s)c%i vid1=%i(%s)c%i vid2=%i(%s)c%i vid3=%i(%s)c%i %s-%s-%s-%s", + Info("Volumes","vid0=%i(%s)c%i vid1=%i(%s)c%i vid2=%i(%s)c%i vid3=%i(%s)c%i %s-%s-%s-%s", vid0,gMC->VolName(vid0),copy0, vid1,gMC->VolName(vid1),copy1, vid2,gMC->VolName(vid2),copy2, @@ -82,72 +78,44 @@ void AliRICHv0::StepManager() Float_t a,z,den,rad,abs; a=z=den=rad=abs=kBad; Int_t mid=gMC->CurrentMaterial(a,z,den,rad,abs); - Info("","mid=%i a=%7.2f z=%7.2f den=%7.2f rad=%7.2f abs=%7.2f",mid,a,z,den,rad,abs); + Info("Material","id=%i a=%7.2f z=%7.2f den=%9.4f rad=%9.2f abs=%9.2f",mid,a,z,den,rad,abs); + + Int_t iTmedId=gMC->GetMedium(); + const char *sTmed; + switch(iTmedId){ + case kAir: sTmed="Air" ;break; + case kRoha: sTmed="Rohacell" ;break; + case kSiO2: sTmed="SiO2" ;break; + case kC6F14: sTmed="C6F14" ;break; + case kGridCu: sTmed="GridCu" ;break; + case kOpSiO2: sTmed="OpSiO2" ;break; + case kGap: sTmed="Gap" ;break; + case kGlass: sTmed="Glass" ;break; + case kCu: sTmed="Cu" ;break; + case kSteel: sTmed="Steel" ;break; + case kPerpex: sTmed="Perpex" ;break; + case kCH4: sTmed="CH4" ;break; + case kCsI: sTmed="CsI" ;break; + case kAl: sTmed="Al" ;break; + case kW: sTmed="W" ;break; + default: sTmed="not known";break; + } + Info("Medium","id=%i (%s)",iTmedId,sTmed); TLorentzVector x4; gMC->TrackPosition(x4); Float_t glo[3],loc[3]; glo[0]=x4.X();glo[1]=x4.Y();glo[2]=x4.Z(); gMC->Gmtod(glo,loc,1); - Info("","glo(%+8.3f,%+8.3f,%+8.3f) r=%8.3f theta=%8.3f phi=%8.3f", + Info("X4 glo","(x,y,z)=(%+8.3f,%+8.3f,%+8.3f) (r,theta,phi)=(%8.3f,%8.3f,%8.3f)", glo[0],glo[1],glo[2],x4.Rho(),x4.Theta()*TMath::RadToDeg(),x4.Phi()*TMath::RadToDeg()); - Info("","loc(%+8.3f,%+8.3f,%8.3f) by gMC->Gmtod()", loc[0],loc[1],loc[2]); - if(gMC->VolId("CSI ")==gMC->CurrentVolID(copy0)){ + Info("X4 loc","(x,y,z)=(%+8.3f,%+8.3f,%8.3f) by gMC->Gmtod()",loc[0],loc[1],loc[2]); + if(gMC->VolId("GAP ")==gMC->CurrentVolID(copy0)){ Int_t iChamber; gMC->CurrentVolOffID(2,iChamber); - TVector2 x2=C(iChamber)->Glob2Loc(x4); - Info("","loc(%+8.3f,%+8.3f) by Glob2Loc", x2.X(),x2.Y()); + TVector2 x2=C(iChamber)->Mrs2Pc(x4); + Info("X4 loc","(%+8.3f,%+8.3f) by Mrs2Pc", x2.X(),x2.Y()); } - Info("","end of current step\n"); + Info(Form("Step %i",iStepN++),"end of current step\n"); }//StepManager() -void AliRICHv0::CreateGeometry() -{ - if(GetDebug())Info("CreateGeometry","Start v0."); - - Double_t cm=1,mm=0.1; - //place radioactive source compartment if needed - Double_t containerLen=21.8*mm/2 ,containerR=38*mm/2; - Double_t screwLen=15*mm/2 ,screwR=5*mm/2; - Double_t srcLen=10*mm/2 ,srcR=2*mm/2; - Double_t perpexLen=20*mm/2 ,perpexR=34*mm/2; - Double_t perpexWholeLen=10*mm/2 ,perpexWholeR=4*mm/2; - Double_t alBottomLen=containerLen-perpexLen ,alWholeR=5*mm/2; -//volumes - Double_t par[3]; - -// Double_t anodWireD=20*mkm, cathWireD=100*mkm, collWireD=100*mkm; - - Double_t pcZ2=0.5*mm; - par[0]=68.8*cm; par[1]=70.86*cm; par[2]=13*cm; gMC->Gsvolu("RICH","BOX ",(*fIdtmed)[kCH4],par,3); //RICH - par[0]=P()->PcSizeX()/2; par[1]=P()->PcSizeY()/2; par[2]=pcZ2; gMC->Gsvolu("CSI ","BOX ",(*fIdtmed)[kCSI],par,3); //CSI - par[0]=P()->PcSizeX()/2; par[1]=P()->PcSizeY()/2; par[2]=P()->GapAmp(); gMC->Gsvolu("GAP ","BOX ",(*fIdtmed)[kCH4],par,3); //GAP -// par[0]=0; par[1]=; par[2]=P()->PcSizeX()/2;gMC->Gsvolu("WIAN","TUBE",(*fIdtmed)[kW],par,3);//Anod wire -// par[0]=0; par[1]=srcR; par[2]=P()->PcSizeX()/2;gMC->Gsvolu("WICA","TUBE",(*fIdtmed)[kCu],par,3);//Cathod wire -// par[0]=0; par[1]=srcR; par[2]=P()->PcSizeX()/2;gMC->Gsvolu("WICO","TUBE",(*fIdtmed)[kCu],par,3);//Collect wire - - par[0]=0 ;par[1]=containerR ;par[2]=containerLen; gMC->Gsvolu("CONT","TUBE",(*fIdtmed)[kCH4] ,par,3); //container - par[0]=perpexR;par[1]=containerR ;par[2]=containerLen; gMC->Gsvolu("ALWA","TUBE",(*fIdtmed)[kAl] ,par,3); //Al cylindric wall - par[0]=0 ;par[1]=perpexR ;par[2]=alBottomLen; gMC->Gsvolu("ALBO","TUBE",(*fIdtmed)[kAl] ,par,3); //Al bottom - par[0]=0 ;par[1]=alWholeR ;par[2]=alBottomLen; gMC->Gsvolu("ALWH","TUBE",(*fIdtmed)[kCH4] ,par,3); //Whole in Al bottom - par[0]=0 ;par[1]=perpexR ;par[2]=perpexLen; gMC->Gsvolu("PEPL","TUBE",(*fIdtmed)[kPerpex],par,3); //Perpex plug - par[0]=0 ;par[1]=perpexWholeR;par[2]=perpexWholeLen;gMC->Gsvolu("PEWH","TUBE",(*fIdtmed)[kCH4] ,par,3); //Whole in Perpex - par[0]=0 ;par[1]=screwR ;par[2]=screwLen; gMC->Gsvolu("SCRE","TUBE",(*fIdtmed)[kSteel] ,par,3); //Screw - par[0]=0 ;par[1]=srcR ;par[2]=srcLen; gMC->Gsvolu("SOUR","TUBE",(*fIdtmed)[kSteel] ,par,3); //Source itself -//nodes - gMC->Gspos("RICH",1,"ALIC",0,0,-9 ,0,"ONLY"); //RICH in ALIC - gMC->Gspos("CSI ",1,"RICH",0,0,-(P()->GapProx()+pcZ2) ,0,"ONLY"); //CsI in ALIC - gMC->Gspos("GAP ",1,"RICH",0,0,-(P()->GapProx()-P()->GapAmp()/2),0,"ONLY"); //GAP in ALIC - - gMC->Gspos("CONT",1,"RICH",0,0,1*cm ,0,"ONLY"); //Sr90 container in RICH - gMC->Gspos("ALWA",1,"CONT",0,0,0 ,0,"ONLY"); //Al wall - gMC->Gspos("ALBO",1,"CONT",0,0,-containerLen+alBottomLen ,0,"ONLY"); //Al bottom - gMC->Gspos("PEPL",1,"CONT",0,0,containerLen-perpexLen ,0,"ONLY"); //Perpex plug - - gMC->Gspos("ALWH",1,"ALBO",6*mm,0,0 ,0,"ONLY"); //Whole in Al bottom - - gMC->Gspos("PEWH",1,"PEPL",6*mm,0,-perpexLen+perpexWholeLen ,0,"ONLY"); //Whole in Perpex plug - gMC->Gspos("SOUR",1,"PEPL",6*mm,0, perpexLen-srcLen ,0,"ONLY"); //Source in Perpex plug - gMC->Gspos("SCRE",1,"PEPL",0,0,perpexLen-screwLen ,0,"ONLY"); //Screw in Perpex plug - if(GetDebug())Info("CreateGeometry","Stop v0."); -}//CreateGeometry() diff --git a/RICH/AliRICHv0.h b/RICH/AliRICHv0.h index 3594ff85b09..bf5ec4b6bb9 100644 --- a/RICH/AliRICHv0.h +++ b/RICH/AliRICHv0.h @@ -15,7 +15,6 @@ public: virtual void Init() {;} //interface from AliRICH virtual Int_t IsVersion() const{return 0;} //interface from AliRICH virtual void StepManager(); //interface from AliRICH - virtual void CreateGeometry(); //interface from AliRICH protected: ClassDef(AliRICHv0,1) //RICH coarse version for material budget study and debuging }; diff --git a/RICH/AliRICHv1.cxx b/RICH/AliRICHv1.cxx index 8d9efdababe..cca590608b9 100644 --- a/RICH/AliRICHv1.cxx +++ b/RICH/AliRICHv1.cxx @@ -32,42 +32,42 @@ ClassImp(AliRICHv1) void AliRICHv1::StepManager() { // Full Step Manager. -// 3- Ckovs absorbed in Collection electrods -// 5- Ckovs absorbed in Cathode wires -// 6- Ckovs absorbed in Anod wires +// 3- Ckovs absorbed on Collection electrods +// 5- Ckovs absorbed on Cathode wires +// 6- Ckovs absorbed on Anod wires Int_t copy; static Int_t iCurrentChamber; //history of Cerenkovs if(gMC->TrackPid()==kCerenkov){ - if(gMC->IsNewTrack() && (gMC->CurrentVolID(copy)==gMC->VolId("FRE1")||gMC->CurrentVolID(copy)==gMC->VolId("FRE2"))) fCounters(0)++;// 0- Ckovs produced in Freon - if(!gMC->IsTrackAlive() && (gMC->CurrentVolID(copy)==gMC->VolId("FRE1")||gMC->CurrentVolID(copy)==gMC->VolId("FRE2"))) fCounters(1)++;// 1- Ckovs absorbed in Freon - if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("QUAR")) fCounters(2)++;//2- Ckovs absorbed in Quartz - if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("META")) fCounters(4)++;//4- Ckovs absorbed in CH4 + if( gMC->IsNewTrack() && gMC->CurrentVolID(copy)==gMC->VolId("RRAD")) fCounters(0)++;// 0- Ckovs produced in radiator + if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("RRAD")) fCounters(1)++;// 1- Ckovs absorbed in radiator + if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("RRWI")) fCounters(2)++;// 2- Ckovs absorbed in radiator window + if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("RICH")) fCounters(4)++;// 4- Ckovs absorbed in CH4 } //Treat photons static TLorentzVector cerX4; - if((gMC->TrackPid()==kCerenkov||gMC->TrackPid()==kFeedback)&&gMC->CurrentVolID(copy)==gMC->VolId("CSI ")){//photon in CSI - if(gMC->Edep()>0){//CF+CSI+DE + if((gMC->TrackPid()==kCerenkov||gMC->TrackPid()==kFeedback)&&gMC->CurrentVolID(copy)==gMC->VolId("RPC ")){//photon in PC + if(gMC->Edep()>0){//photon in PC +DE if(IsLostByFresnel()){ - if(gMC->TrackPid()==kCerenkov) fCounters(7)++;// 7- Ckovs reflected on CsI + if(gMC->TrackPid()==kCerenkov) fCounters(7)++;// 7- Ckovs reflected from CsI gMC->StopTrack(); return; } - gMC->TrackPosition(cerX4); gMC->CurrentVolOffID(2,iCurrentChamber); + gMC->TrackPosition(cerX4); gMC->CurrentVolOffID(2,iCurrentChamber);//RICH-RPPF-RPC AddHit(iCurrentChamber,gAlice->GetMCApp()->GetCurrentTrackNumber(),cerX4.Vect(),cerX4.Vect());//HIT for PHOTON in conditions CF+CSI+DE - fCounters(8)++;//4- Ckovs converted in CsI + fCounters(8)++;//4- Ckovs converted to electron on CsI GenerateFeedbacks(iCurrentChamber); - }//CF+CSI+DE - }//CF in CSI + }//photon in PC and DE >0 + }//photon in PC //Treat charged particles static Float_t eloss; static TLorentzVector mipInX4,mipOutX4; - if(gMC->TrackCharge() && gMC->CurrentVolID(copy)==gMC->VolId("GAP ")){//MIP in GAP - gMC->CurrentVolOffID(3,iCurrentChamber); + if(gMC->TrackCharge() && gMC->CurrentVolID(copy)==gMC->VolId("RGAP")){//MIP in GAP + gMC->CurrentVolOffID(1,iCurrentChamber);//RICH-RGAP if(gMC->IsTrackEntering()||gMC->IsNewTrack()) {//MIP in GAP entering or newly created eloss=0; gMC->TrackPosition(mipInX4); @@ -93,7 +93,7 @@ Bool_t AliRICHv1::IsLostByFresnel() Double_t localTheta = TMath::ATan2(TMath::Sqrt(localTc),localMom[1]); Double_t cotheta = TMath::Abs(TMath::Cos(localTheta)); if(gMC->GetRandom()->Rndm() < Fresnel(p4.E()*1e9,cotheta,1)){ - if(GetDebug()) Info("IsLostByFresnel","Photon lost"); + AliDebug(1,"Photon lost"); return kTRUE; }else return kFALSE; @@ -106,11 +106,12 @@ void AliRICHv1::GenerateFeedbacks(Int_t iChamber,Float_t eloss) TLorentzVector x4; gMC->TrackPosition(x4); - TVector2 x2=C(iChamber)->Glob2Loc(x4); - Int_t sector=P()->Loc2Sec(x2); if(sector==kBad) return; //hit in dead zone, nothing to produce + TVector2 x2=C(iChamber)->Mrs2Pc(x4);//hit position on photocathode plane + TVector2 xspe=x2; + Int_t sector=P()->Loc2Sec(xspe); if(sector==kBad) return; //hit in dead zone, nothing to produce Int_t iTotQdc=P()->TotQdc(x2,eloss); - Int_t iNphotons=gMC->GetRandom()->Poisson(P()->AlphaFeedback(sector)*iTotQdc); - if(GetDebug())Info("GenerateFeedbacks","N photons=%i",iNphotons); + Int_t iNphotons=gMC->GetRandom()->Poisson(P()->C(iChamber)->AlphaFeedback(sector)*iTotQdc); + 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]; //Generate photons @@ -172,5 +173,5 @@ void AliRICHv1::GenerateFeedbacks(Int_t iChamber,Float_t eloss) outputNtracksStored, 1.0); }//feedbacks loop - if(GetDebug()) Info("GenerateFeedbacks","Stop."); + AliDebug(1,"Stop."); }//GenerateFeedbacks() diff --git a/RICH/CreateConfig.C b/RICH/CreateConfig.C index 1634afcb2c2..a1b5e53bac5 100644 --- a/RICH/CreateConfig.C +++ b/RICH/CreateConfig.C @@ -6,7 +6,7 @@ public: kPMD=8192,kDIPO=16384,kEMCAL=32768,kVZERO=65536,kMUON=131072,kZDC=262144,kSHILD=524288}; enum EProcesses {kDCAY=1,kPAIR=2,kCOMP=4,kPHOT=8,kPFIS=16,kDRAY=32,kANNI=64,kBREM=128,kMUNU=256,kCKOV=512,kHADR=1024,kLOSS=2048,kMULS=4096, kRAYL=8192}; - enum EGenTypes {kGun1,kGun7,kPythia7,kHijing,kHijingPara}; + enum EGenTypes {kGun1=100,kGun7,kPythia7,kHijing,kHijingPara}; KirConfig(const char*sFileName); ~KirConfig() {Info("ctor","");fMain->Cleanup(); delete fMain;} @@ -51,9 +51,12 @@ KirConfig::KirConfig(const char *sFileName) //Generator pVerFrame->AddFrame(pGenGrpFrm=new TGGroupFrame(pHorFrame,"Generator")); pGenGrpFrm->AddFrame(fGenTypeCombo=new TGComboBox(pGenGrpFrm,100)); - fGenTypeCombo->AddEntry("gun to single chamber",kGun1); fGenTypeCombo->AddEntry("gun to all chambers",kGun7); + fGenTypeCombo->AddEntry("gun to single chamber",kGun1); + fGenTypeCombo->AddEntry("gun to all chambers",kGun7); fGenTypeCombo->AddEntry("7 guns on top of Pythia",kPythia7); - fGenTypeCombo->AddEntry("HIJING",kHijing); fGenTypeCombo->AddEntry("parametrized HIJING",kHijingPara); + fGenTypeCombo->AddEntry("HIJING",kHijing); + fGenTypeCombo->AddEntry("parametrized HIJING",kHijingPara); + fGenTypeCombo->AddEntry("Sr90 source",kSr90); fGenTypeCombo->Select(kHijingPara); fGenTypeCombo->Resize(160,20); @@ -86,7 +89,7 @@ KirConfig::KirConfig(const char *sFileName) // Magnetic Field pVerFrame->AddFrame(pFldGrpFrm=new TGGroupFrame(pHorFrame,"Magnetic Field")); pFldGrpFrm->AddFrame(fMagFldChkBtn=new TGCheckButton(pFldGrpFrm,"On/Off")); - fMagFldChkBtn->SetState(kButtonUp); + fMagFldChkBtn->SetState(kButtonDown); //Detectors pHorFrame->AddFrame(fDetButGrp=new TGButtonGroup(pHorFrame,"Detectors")); fDetButGrp->Connect("Pressed(Int_t)","KirConfig",this,"AddDetector(Int_t)"); @@ -206,14 +209,15 @@ void KirConfig::CreateConfigFile() switch(fGenTypeCombo->GetSelected()){ case kHijingPara: fprintf(fp," AliGenHIJINGpara *pGen=new AliGenHIJINGpara(%i);\n",(int)fGenNprimEntry->GetNumber()); - fprintf(fp," pGen->SetMomentumRange(0,999); pGen->SetPhiRange(0,360); pGen->SetThetaRange(%f,%f);\n",Eta2Theta(8),Eta2Theta(-8)); + fprintf(fp," pGen->SetMomentumRange(0,999); pGen->SetThetaRange(%f,%f); pGen->SetPhiRange(0,360);\n",Eta2Theta(8),Eta2Theta(-8)); fprintf(fp," pGen->SetOrigin(0,0,0); pGen->SetSigma(0,0,0);\n"); fprintf(fp," pGen->Init();\n"); break; case kGun1: fprintf(fp," AliGenFixed *pGen=new AliGenFixed(1);\n"); - fprintf(fp," pGen->SetPart(%i); pGen->SetMomentum(%3.1f); pGen->SetOrigin(0,0,0);\n",fGenPartIdCombo->GetSelected(),float(fGenMomCombo->GetSelected())/10); - fprintf(fp," pGen->SetPhiRange(pRICH->C(%i)->PhiD()); pGen->SetThetaRange(pRICH->C(%i)->ThetaD()-2);\n",fGenChamberCombo->GetSelected(),fGenChamberCombo->GetSelected()); + fprintf(fp," pGen->SetPart(%i); pGen->SetOrigin(0,0,0);\n",fGenPartIdCombo->GetSelected()); + fprintf(fp," pGen->SetMomentum(%3.1f); pGen->SetTheta(pRICH->C(%i)->ThetaD()-2); pGen->SetPhi(pRICH->C(%i)->PhiD()); \n", + float(fGenMomCombo->GetSelected())/10,fGenChamberCombo->GetSelected(),fGenChamberCombo->GetSelected()); fprintf(fp," pGen->Init();\n"); break; case kGun7: @@ -221,7 +225,7 @@ void KirConfig::CreateConfigFile() fprintf(fp," for(int i=1;i<=7;i++){\n"); fprintf(fp," AliGenFixed *pFixed=new AliGenFixed(1);\n"); fprintf(fp," pFixed->SetPart(%i); pFixed->SetMomentum(2.5+i*0.4); pFixed->SetOrigin(0,0,0);\n",fGenPartIdCombo->GetSelected()); - fprintf(fp," pFixed->SetPhiRange(pRICH->C(i)->PhiD()); pFixed->SetThetaRange(pRICH->C(i)->ThetaD()-2);\n"); + fprintf(fp," pFixed->SetPhi(pRICH->C(i)->PhiD()); pFixed->SetTheta(pRICH->C(i)->ThetaD()-2);\n"); fprintf(fp," pCocktail->AddGenerator(pFixed,Form(\"Fixed %i\",i),1);\n }\n"); fprintf(fp," pCocktail->Init();\n"); break; @@ -238,6 +242,18 @@ void KirConfig::CreateConfigFile() fprintf(fp," pPythia->SetProcess(kPyMb); pPythia->SetEnergyCMS(14000);\n"); fprintf(fp," pCocktail->AddGenerator(pPythia,\"Pythia\",1);\n"); fprintf(fp," pCocktail->Init();\n"); + case kSr90: + fprintf(fp," AliGenRadioactive *pGen=new AliGenRadioactive(kSr90,%i);\n",(int)fGenNprimEntry->GetNumber()); + fprintf(fp," pGen->SetOrigin(0,0,0); pGen->SetSigma(0.1,0.1,0.05);\n"); + fprintf(fp," pGen->Init();\n"); + break; + case kHijing: + fprintf(fp," AliGenHijing *pGen=new AliGenHijing(-1); pGen->SetEnergyCMS(5500); pGen->SetReferenceFrame(\"CMS\");\n"); + fprintf(fp," pGen->SetProjectile(\"A\", 208, 82); pGen->SetTarget(\"A\", 208, 82);\n"); + fprintf(fp," pGen->SetJetQuenching(0); pGen->SetShadowing(0);\n"); + fprintf(fp," pGen->KeepFullEvent(); pGen->SetSelectAll(0); \n"); + fprintf(fp," pGen->SetImpactParameterRange(0., 5.);\n"); + fprintf(fp," pGen->Init();\n"); break; } //central before RICH detectors diff --git a/RICH/Geom.C b/RICH/Geom.C index d831c11c8b9..d0f46e45cfb 100644 --- a/RICH/Geom.C +++ b/RICH/Geom.C @@ -11,17 +11,20 @@ void Geom() Double_t dx,dy,dz,r1,r2,a=0,z=0,den=0,radlen=0,intlen=0; Double_t cm=1,m=100*cm,mm=0.1*cm,mkm=0.001*cm; //Medium - pAir=new TGeoMedium("Air" ,1, g->GetMaterial("Air")); - new TGeoMaterial("CH4" ,a=26.98 ,z=13,den=0.4224 ); pCH4=new TGeoMedium("CH4" ,2, g->GetMaterial("CH4")); - new TGeoMaterial("CsI" ,a=26.98 ,z=13,den=2.7 ,radlen=24.01,intlen=70.6); pCsI=new TGeoMedium("CsI" ,3, g->GetMaterial("CsI")); - new TGeoMaterial("Al" ,a=26.98 ,z=13,den=2.7 ,radlen=24.01,intlen=70.6); pAl =new TGeoMedium("Al" ,4, g->GetMaterial("Al")); - new TGeoMaterial("Fe" ,a=55.845,z=26,den=7.87,radlen=13.84,intlen=82.8); pFe =new TGeoMedium("Fe" ,5, g->GetMaterial("Fe")); - new TGeoMaterial("W" ,a=183.84,z=27,den=19.3,radlen= 9.59*cm,intlen=0.35*cm);pW =new TGeoMedium("W" ,6, g->GetMaterial("W")); - new TGeoMaterial("Cu" ,a=55.845,z=26,den=7.87,radlen=13.84,intlen=82.8); pCu =new TGeoMedium("Cu" ,7, g->GetMaterial("Cu")); - new TGeoMaterial("Perpex" ,a=55.845,z=26,den=7.87,radlen=13.84,intlen=82.8); new TGeoMedium("Perpex" ,8, g->GetMaterial("Perpex")); - new TGeoMaterial("Rohacell" ,a=12.01,z=6,den=0.1,radlen=18.8,intlen=0); new TGeoMedium("Rohacell" ,9, g->GetMaterial("Rohacell")); - new TGeoMaterial("SiO2" ,a=0,z=0,den=0); new TGeoMedium("SiO2" ,10, g->GetMaterial("SiO2")); - new TGeoMaterial("C6F14" ,a=0,z=0,den=0); new TGeoMedium("C6F14" ,11, g->GetMaterial("C6F14")); + pAir =new TGeoMedium("Air" ,1, g->GetMaterial("Air")); + new TGeoMaterial("CH4" ,a=26.98 ,z=13,den=0.4224 ); pCH4 =new TGeoMedium("CH4" ,2, g->GetMaterial("CH4")); + new TGeoMaterial("CsI" ,a=26.98 ,z=13,den=2.7 ,radlen=24.01*cm,intlen=70.6*cm); pCsI =new TGeoMedium("CsI" ,3, g->GetMaterial("CsI")); + new TGeoMaterial("Al" ,a=26.98 ,z=13,den=2.7 ,radlen=24.01*cm,intlen=70.6*cm); pAl =new TGeoMedium("Al" ,4, g->GetMaterial("Al")); + new TGeoMaterial("W" ,a=183.84,z=27,den=19.3,radlen= 9.59*cm,intlen=0.35*cm); pW =new TGeoMedium("W" ,5, g->GetMaterial("W")); + new TGeoMaterial("Cu" ,a=55.845,z=26,den=7.87,radlen=13.84*cm,intlen=82.8*cm); pCu =new TGeoMedium("Cu" ,6, g->GetMaterial("Cu")); + new TGeoMaterial("Rohacell",a=12.01 ,z=6 ,den=0.1,radlen=18.8,intlen=0); pRoha =new TGeoMedium("Rohacell" ,7, g->GetMaterial("Rohacell")); + new TGeoMaterial("SiO2" ,a=0 ,z=0 ,den=0); pSiO2= new TGeoMedium("SiO2" ,8,g->GetMaterial("SiO2")); + new TGeoMaterial("C6F14" ,a=0 ,z=0 ,den=0); pC6F14=new TGeoMedium("C6F14" ,9,g->GetMaterial("C6F14")); +//for Sr90 source + new TGeoMaterial("Perpex" ,a=55.845,z=26,den=7.87,radlen=13.84*cm,intlen=82.8*cm); pPerpex=new TGeoMedium("Perpex",10, g->GetMaterial("Perpex")); + new TGeoMaterial("Steel" ,a=55.845,z=26,den=7.87,radlen=13.84*cm,intlen=82.8*cm); pSteel=new TGeoMedium("Steel" ,11, g->GetMaterial("Steel")); + new TGeoMaterial("Mylar" ,a=55.845,z=26,den=7.87,radlen=13.84*cm,intlen=82.8*cm); pMylar=new TGeoMedium("Mylar" ,12, g->GetMaterial("Mylar")); + new TGeoMaterial("Sr" ,a=87.62 ,z=38,den=7.87,radlen=13.84*cm,intlen=82.8*cm); pSr =new TGeoMedium("Sr" ,13, g->GetMaterial("Sr")); //Mother TGeoVolume *pMother=g->MakeBox("Mother",pAir,dx=10*m/2,dy=10*m/2,dz=30*m/2); //arbitrary values //Rich @@ -43,7 +46,7 @@ void Geom() pRich->AddNode(pPpf,copy=5,new TGeoTranslation(-335*mm,+433*mm,8*cm+20*mm)); pRich->AddNode(pPpf,copy=6,new TGeoTranslation(+335*mm,+433*mm,8*cm+20*mm)); pPpf->AddNode( pPc ,copy=1,new TGeoTranslation( 0*mm, 0*mm,-19.15*mm));//PPF 2001P2 - pPpf->AddNode(g->GetVolume("PPFlarge"),copy=1,new TGeoTranslation(-224.5*mm,-151.875*mm, 0.85*mm)); + pPpf->AddNode(pPpfLarge,copy=1,new TGeoTranslation(-224.5*mm,-151.875*mm, 0.85*mm)); pPpf->AddNode(g->GetVolume("PPFlarge"),copy=2,new TGeoTranslation(-224.5*mm,- 50.625*mm, 0.85*mm)); pPpf->AddNode(g->GetVolume("PPFlarge"),copy=3,new TGeoTranslation(-224.5*mm,+ 50.625*mm, 0.85*mm)); pPpf->AddNode(g->GetVolume("PPFlarge"),copy=4,new TGeoTranslation(-224.5*mm,+151.875*mm, 0.85*mm)); @@ -99,52 +102,49 @@ void Geom() TGeoVolume *pRadShort =g->MakeBox( "RadShort" ,g->GetMedium("Neoceram") ,dx= 10*mm/2 ,dy= 403*mm/2 ,dz= 15*mm/2); TGeoVolume *pRadSpacer=g->MakeTube("RadSpacer",g->GetMedium("SiO2") ,r1= 0 ,r2=10*mm/2 ,dz= 15*mm/2); - pRich->AddNode(pRad,copy=1,new TGeoTranslation( 0*mm,-434*mm, -12*mm)); - pRich->AddNode(pRad,copy=2,new TGeoTranslation( 0*mm, 0*mm, -12*mm)); - pRich->AddNode(pRad,copy=3,new TGeoTranslation( 0*mm,+434*mm, -12*mm)); - pRad->AddNode(g->GetVolume("RadFront") ,copy=1,new TGeoTranslation( 0*mm, 0*mm, -10.0*mm)); - pRad->AddNode(g->GetVolume("RadWin") ,copy=1,new TGeoTranslation( 0*mm, 0*mm, 9.5*mm)); - pRad->AddNode(g->GetVolume("RadLong") ,copy=1,new TGeoTranslation( 0*mm,-204*mm, -0.5*mm)); - pRad->AddNode(g->GetVolume("RadLong") ,copy=2,new TGeoTranslation( 0*mm,+204*mm, -0.5*mm)); - pRad->AddNode(g->GetVolume("RadShort") ,copy=1,new TGeoTranslation(-660*mm, 0*mm, -0.5*mm)); - pRad->AddNode(g->GetVolume("RadShort") ,copy=2,new TGeoTranslation(+660*mm, 0*mm, -0.5*mm)); + pRich->AddNode(pRad ,copy=1,new TGeoTranslation( 0*mm,-434*mm, -12*mm)); + pRich->AddNode(pRad ,copy=2,new TGeoTranslation( 0*mm, 0*mm, -12*mm)); + pRich->AddNode(pRad ,copy=3,new TGeoTranslation( 0*mm,+434*mm, -12*mm)); + + pRad ->AddNode(pRadFront ,copy=1,new TGeoTranslation( 0*mm, 0*mm, -10.0*mm)); + pRad ->AddNode(pRadWin ,copy=1,new TGeoTranslation( 0*mm, 0*mm, 9.5*mm)); + pRad ->AddNode(pRadLong ,copy=1,new TGeoTranslation( 0*mm,-204*mm, -0.5*mm)); + pRad ->AddNode(pRadLong ,copy=2,new TGeoTranslation( 0*mm,+204*mm, -0.5*mm)); + pRad ->AddNode(pRadShort ,copy=1,new TGeoTranslation(-660*mm, 0*mm, -0.5*mm)); + pRad ->AddNode(pRadShort ,copy=2,new TGeoTranslation(+660*mm, 0*mm, -0.5*mm)); for(int i=0;i<3;i++) for(int j=0;j<10;j++) pRad->AddNode(pRadSpacer,copy=10*i+j,new TGeoTranslation(-1330*mm/2+116*mm+j*122*mm,(i-1)*105*mm,-0.5*mm)); }else{//Sr90 - Double_t containerZ=21.8*mm/2,containerR=38*mm/2; - Double_t screwLen=15*mm/2 ,screwR=5*mm/2; - Double_t srcLen=10*mm/2 ,srcR=2*mm/2; - Double_t perpexZ=20*mm/2 ,perpexR=34*mm/2; - Double_t perpexWholeLen=10*mm/2,perpexWholeR=4*mm/2; - Double_t alBottomLen=containerZ-perpexZ,alWholeR=5*mm/2; - g->MakeTube("Sr90" ,pAir ,0 ,containerR ,containerZ); - g->MakeTube("Sr90AlWall" ,pAl ,perpexR,containerR ,containerZ); - g->MakeTube("Sr90AlBottom" ,pAl ,0 ,perpexR ,alBottomLen); - g->MakeTube("Sr90AlWhole" ,pAir ,0 ,alWholeR ,alBottomLen); - g->MakeTube("Sr90PerpexPlug" ,g->GetMedium("Perpex"),0 ,perpexR ,perpexZ); - g->MakeTube("Sr90PerpexWhole",pAir ,0 ,perpexWholeR,perpexWholeLen); - g->MakeTube("Sr90Screw" ,g->GetMedium("Steel") ,0 ,screwR ,screwLen); - g->MakeTube("Sr90Source" ,g->GetMedium("Steel") ,0 ,srcR ,srcLen); + pSrc =g->MakeTube("Src" ,pCH4 , 0 , 70*mm/2 , 30*mm/2); //top container + pAlGlass =g->MakeTube("SrcAlGlass" ,pAl , 0 , 38*mm/2 ,21.8*mm/2); //Al glass wall + pPerpexPlug =g->MakeTube("SrcPerpex" ,pPerpex , 0 , 34*mm/2 , 20*mm/2); //Perpex plug + pSteelScrew =g->MakeTube("SrcSteelScrew" ,pSteel , 0 , 5*mm/2 , 15*mm/2); //Steel screw in the center + pSr90Screw =g->MakeTube("SrcSr90Screw" ,pSteel , 0 , 2*mm/2 , 10*mm/2); //Steel screw to support Sr90 + pSr90 =g->MakeTube("SrcSr90" ,pSr , 0 , 1*mm/2 , 1*mm/2); //Sr90 source + pHolePerpex =g->MakeTube("SrcHolePerpex" ,pAir , 0 , 4*mm/2 , 10*mm/2); //Air hole in perpex plug + pHoleAl =g->MakeTube("SrcHoleAl" ,pAir , 0 , 4*mm/2 , 1.8*mm/2); //Air hole in Al glass bottom + pMylarFoil =g->MakeTube("SrcMylarFoil" ,pMylar , 0 , 30*mm/2 , 50*mkm/2); //Mylar foil - pRich->AddNode(pContainer,1,new TGeoTranslation(30*cm,0,containerZ)); - pContainer ->AddNode(pAlWall,1); - pContainer ->AddNode(pAlBottom,1,new TGeoTranslation(0,0,-containerZ+alBottomLen)); - pContainer ->AddNode(pPerpexPlug,1,new TGeoTranslation(0,0,containerZ-perpexZ)); - pAlBottom ->AddNode(pAlWhole,1,new TGeoTranslation(6*mm,0,0)); - pPerpexPlug->AddNode(pPerpexWhole,1,new TGeoTranslation(6*mm,0,-perpexZ+perpexWholeLen)); - pPerpexPlug->AddNode(pSource,1 ,new TGeoTranslation(6*mm,0, perpexZ-srcLen)); - pPerpexPlug->AddNode(pScrew,1,new TGeoTranslation(0,0,perpexZ-screwLen)); + pRich->AddNode(pSrc,1,new TGeoTranslation(30*cm,0,alVesselZ)); + pSrc ->AddNode(pMylarFoil,1,new TGeoTranslation(0,0,21.8*mm/2+50*mkm/2)); + pSrc ->AddNode(pAlGlass,1,new TGeoTranslation(0,0,0));//Al glass to fake Src volume + pAlGlass->AddNode(pHoleAl,1,new TGeoTranslation(6*mm, 0*mm, 10.9*mm)); + pAlGlass->AddNode(pPerpexPlug,1,new TGeoTranslation(0*mm,0*mm,-0.9*mm)); + pPerpexPlug->AddNode(pSteelScrew,1,new TGeoTranslation(shift,0, perpexPlugZ-steelScrewZ)); + pPerpexPlug->AddNode(pHolePerpex,1,new TGeoTranslation(shift,0,-perpexPlugZ+airWholeZ)); + pPerpexPlug->AddNode(pSr90Screw ,1,new TGeoTranslation(shift,0, perpexPlugZ-sr90ScrewZ)); + pSr90Screw->AddNode(pSr90 ,1,new TGeoTranslation(0,0, -(sr90ScrewZ-sr90Z))); }//if(!isSr90) //Sandbox - TGeoVolume *pSandBox=g->MakeBox( "SandBox" ,g->GetMedium("Air") ,dx=1419*mm/2 ,dy=1378*mm/2 ,dz=50.5*mm/2); //2072P1 - g->MakeBox( "SandCover",g->GetMedium("Al") ,dx=1419*mm/2 ,dy=1378*mm/2 ,dz= 0.5*mm/2); - g->MakeBox( "SandComb" ,g->GetMedium("Rohacell") ,dx=1359*mm/2 ,dy=1318*mm/2 ,dz=49.5*mm/2); + TGeoVolume *pSandBox =g->MakeBox( "SandBox" ,pAir ,dx=1419*mm/2 ,dy=1378*mm/2 ,dz=50.5*mm/2); //2072P1 + TGeoVolume *pSandCover=g->MakeBox( "SandCover",pAl ,dx=1419*mm/2 ,dy=1378*mm/2 ,dz= 0.5*mm/2); + TGeoVolume *pSandComb =g->MakeBox( "SandComb" ,pRoha ,dx=1359*mm/2 ,dy=1318*mm/2 ,dz=49.5*mm/2); pRich->AddNode(pSandBox,copy=1,new TGeoTranslation( 0*mm,0*mm, -73.75*mm)); //p.84 TDR - pSandBox->AddNode(g->GetVolume("SandComb") ,copy=1,new TGeoTranslation( 0*mm,0*mm, 0*mm)); //2072P1 - pSandBox->AddNode(g->GetVolume("SandCover") ,copy=1,new TGeoTranslation( 0*mm,0*mm, +25*mm)); - pSandBox->AddNode(g->GetVolume("SandCover") ,copy=2,new TGeoTranslation( 0*mm,0*mm, -25*mm)); + pSandBox->AddNode(pSandComb ,copy=1,new TGeoTranslation( 0*mm,0*mm, 0*mm)); //2072P1 + pSandBox->AddNode(pSandCover ,copy=1,new TGeoTranslation( 0*mm,0*mm, +25*mm)); + pSandBox->AddNode(pSandCover ,copy=2,new TGeoTranslation( 0*mm,0*mm, -25*mm)); //Colors g->GetVolume("Pc") ->SetLineColor(kGreen); g->GetVolume("PPF") ->SetLineColor(33); diff --git a/RICH/Opticals.h b/RICH/Opticals.h index a5143b52ed3..669fad691ec 100644 --- a/RICH/Opticals.h +++ b/RICH/Opticals.h @@ -15,12 +15,11 @@ void Opticals() Float_t aIdxSiO2[kNbins]; Float_t aIdxOpSiO2[kNbins]; - Float_t e1=10.666,e2=18.125,f1=46.411,f2= 228.71; //RICH TDR page 35 for (i=0;iSetDebug(1); + gSystem->Exec("rm -rf *.root hlt hough ZZZ*"); + if(isDebug) AliLog::SetGlobalDebugLevel(AliLog::kDebug); - Info("my/RichBatch.C","%i event(s) requested, debug %i,config file %s",iNevents,isDebug,sConfigFileName); TStopwatch sw;TDatime time; - + + AliSimulation *pSim=new AliSimulation; pSim->Run(iNevents); delete pSim; - AliReconstruction *pRec=new AliReconstruction; pRec->SetRunLocalReconstruction("RICH"); pRec->SetFillESD("RICH"); pRec->Run(); delete pRec; + AliReconstruction *pRec=new AliReconstruction; + pRec->SetRunLocalReconstruction("ITS TPC TRD TOF RICH"); + pRec->SetFillESD("ITS TPC TRD TOF RICH"); + pRec->Run(); delete pRec; + + //pRec->SetRunLocalReconstruction("RICH"); pRec->SetFillESD("RICH"); + + cout<<"\n!!!!!!!!!!!!Info in : Start creating Control Plots \n"; - AliRunLoader* pAL = AliRunLoader::Open(); - pAL->LoadgAlice(); + AliRunLoader* pAL = AliRunLoader::Open(); pAL->LoadgAlice(); AliRICH *pRICH=(AliRICH*)pAL->GetAliRun()->GetDetector("RICH"); if(pRICH) pRICH->ControlPlots(); - cout<<"\nInfo in : Start time: ";time.Print(); - cout<<"Info in : Stop time: ";time.Set(); time.Print(); - cout<<"Info in : Time used: ";sw.Print(); + cout<<"\n!!!!!!!!!!!!Info in : Start time: ";time.Print(); + cout<<"!!!!!!!!!!!!Info in : Stop time: ";time.Set(); time.Print(); + cout<<"!!!!!!!!!!!!Info in : Time used: ";sw.Print(); gSystem->Exec("touch ZZZ______finished_______ZZZ"); +// gSystem->Exec("playwave .kde/my/end.wav"); } diff --git a/RICH/api.txt b/RICH/api.txt index 95c4933a2e9..0e059c75101 100644 --- a/RICH/api.txt +++ b/RICH/api.txt @@ -3,7 +3,7 @@ How to open session? How to retrive pointer to alice run loader: use pRICH->GetLoader()->GetRunLoader() (all detector classes inherit from AliDetector wich has GetLoader()) use methode AliRun::GetRunLoader for gAlice (depricated) -How to get pointers to deifferent root trees: +How to get pointers to different root trees: TreeE belongs to AliRunLoader, available after AliRunLoader::LoadHeader() TreeK belongs to AliRunLoader, available after AliRunLoader::LoadKinematics() TreeH belongs to AliLoader , available after AliLoader::LoadHits() @@ -13,19 +13,20 @@ How to get pointers to deifferent root trees: How to work with the stack of particles? - - pointer to the stack is returned by gAlice->Stack() (global gAlice of type AliRun) or AliRunLoader::Stack() but - before one needs to load event header by AliRunLoader::LoadHeader() otherwise both methods return 0. - Moreover loading header gives the information about number of particles only. - To retrive the list of particle one also needs to load kinematics by AliRunLoader::LoadKinematics() - - total amount of particles in stack for a given event: AliRunLoader::Stack()->GetNtrack() or AliRun::GetEvent() (after LoadHeader()) - - total amount of primiry particles in stack for a given event: AliRunLoader::Stack()->GetNprimary() or AliLoader::TreeH()->GetEntries() (after LoadHeader()) + - first of all, the stack includes primary as well as secondary particles + - pointer to the stack is returned by AliRun::Stack() (global gAlice of type AliRun - depricated way to do) + or by AliRunLoader::Stack() but before one needs to load event header by AliRunLoader::LoadHeader() otherwise both methods return 0. + Moreover loading header gives the information about number of particles only. + To retrive the list of particle one also needs to load kinematics by AliRunLoader::LoadKinematics() + - total amount of particles in stack for a given event: AliStack::GetNtrack() or AliRun::GetEvent() (after LoadHeader()) + - total amount of primiry particles in stack for a given event: AliStack::GetNprimary() or AliLoader::TreeH()->GetEntries() (after LoadHeader()) How to retrive hits: - Hits a stored on primiry by primiry basis. To retrieve all hits one needs to do: + Hits a stored on primary by primary basis. To retrieve all hits one needs to do: initialise the root tree and containers: AliLoader::LoadHits() - read number of primiries in current event: + read number of primaries in current event: loop on the list of primiries: @@ -36,13 +37,8 @@ How to retrive sdigits? One needs to say: pRich->GetLoader()->LoadSDigits(); this one open file, get the tree and invoke AliRICH::SetTreeAddress() - - -gAlice->GetMCApp()->GetCurrentTrackNumber() - - What are the debug methodes avail: - AliModule::GetDebug() + AliModule::GetDebug() now depricated, use AliDebug printout instead AliModule::SetDebug() AliRun::GetDebug() AliRun::SetDebug() @@ -61,5 +57,3 @@ How to deal with AliRunDigitizer? How to avoid using gAlice? Rich()->GetLoader()->GetRunLoader()->GetAliRun() returns gAlice global pointer. - - diff --git a/RICH/libRICH.pkg b/RICH/libRICH.pkg index b3a7f3460f4..4b1a09568f2 100644 --- a/RICH/libRICH.pkg +++ b/RICH/libRICH.pkg @@ -3,7 +3,7 @@ SRCS = AliRICH.cxx AliRICHv0.cxx AliRICHv1.cxx\ AliRICHClusterFinder.cxx AliRICHMap.cxx\ AliRICHChamber.cxx AliRICHRecon.cxx\ AliRICHDisplFast.cxx\ - AliRICHDigitizer.cxx AliRICHReconstructor.cxx + AliRICHDigitizer.cxx AliGenRadioactive.cxx AliRICHHelix.cxx AliRICHTracker.cxx AliRICHReconstructor.cxx HDRS = $(SRCS:.cxx=.h) DHDR= RICHLinkDef.h diff --git a/RICH/menu.C b/RICH/menu.C index c312014c493..10e3d182163 100644 --- a/RICH/menu.C +++ b/RICH/menu.C @@ -3,6 +3,7 @@ void ph(Int_t event=0) {R()->PrintHits(event);} //utility print hits for 'eve void ps(Int_t event=0) {R()->PrintSDigits(event);} //utility print sdigits void pd(Int_t event=0) {R()->PrintDigits(event);} //utility print digits void pc(Int_t event=0) {R()->PrintClusters(event);}//utility print clusters +void pt(Int_t event=0) {R()->PrintTracks(event);} //utility print tracks //__________________________________________________________________________________________________ void pp(int tid) @@ -19,6 +20,7 @@ void pp(int tid) //__________________________________________________________________________________________________ void PrintParticleInfo(int tid) { +// Prints particle info for a given TID TParticle *p=al->Stack()->Particle(tid); cout<GetName()<<"("<GetMother(0)!=-1){cout<<" from "; PrintParticleInfo(p->GetMother(0));} @@ -27,6 +29,7 @@ void PrintParticleInfo(int tid) //__________________________________________________________________________________________________ Int_t prim(Int_t tid) { +// Provides mother TID for the given TID R()->GetLoader()->GetRunLoader()->LoadHeader(); R()->GetLoader()->GetRunLoader()->LoadKinematics(); if(tid<0||tid>=R()->GetLoader()->GetRunLoader()->Stack()->GetNtrack()) @@ -147,13 +150,11 @@ AliRICH *r; Bool_t ReadAlice() { Info("ReadAlice","Tring to read ALICE from SIMULATED FILE."); - AliLoader::SetDebug(0); if(gAlice) delete gAlice; if(!(al=AliRunLoader::Open("galice.root","AlicE","update"))){ gSystem->Exec("rm -rf *.root *.dat"); Error("menu.C::ReadAlice","galice.root broken, removing all this garbage then init new one"); new AliRun("gAlice","Alice experiment system"); - gAlice->SetDebug(-1); gAlice->Init("Config.C"); r=(AliRICH*)gAlice->GetDetector("RICH"); return kFALSE; @@ -164,7 +165,6 @@ Bool_t ReadAlice() a->SetDebug(0); //RICH if(!(r=(AliRICH*)gAlice->GetDetector("RICH"))) Warning("RICH/menu.C::ReadAlice","No RICH in file"); - r->SetDebug(0); if(!(rl=al->GetLoader("RICHLoader"))) Warning("RICH/menu.C::ReadAlice","No RICH loader in file"); Info("ReadAlice","Run contains %i event(s)",gAlice->GetEventsPerRun()); @@ -353,22 +353,22 @@ void TestSeg() TPolyLine3D *pXaxis=new TPolyLine3D(2,X);pXaxis->SetLineColor(kRed); pXaxis->Draw(); TPolyLine3D *pYaxis=new TPolyLine3D(2,Y);pYaxis->SetLineColor(kGreen); pYaxis->Draw(); TPolyLine3D *pZaxis=new TPolyLine3D(2,Z);pZaxis->SetLineColor(kBlue); pZaxis->Draw(); -//chambers +//Draw PC for all chambers by trasfering Pc plane using Pc2Mrs methode Int_t iNpointsX=50,iNpointsY=50; for(Int_t iChamberN=1;iChamberN<=7;iChamberN++){//chamber loop TPolyMarker3D *pChamber=new TPolyMarker3D(iNpointsX*iNpointsY); Int_t i=0; for(Double_t x=0;xLoc2Glob(TVector2(x,y)); - pChamber->SetPoint(i++,v3.X(),v3.Y(),v3.Z()); - pChamber->SetMarkerSize(1); - pChamber->SetMarkerColor(iChamberN); + TVector3 v3=p.C(iChamberN)->Pc2Mrs(TVector2(x,y));//from regular grid of local PC points to MRS presentation + pChamber->SetPoint(i++,v3.X(),v3.Y(),v3.Z());//Pc plane poing in MRS }//step loop + pChamber->SetMarkerSize(1); + pChamber->SetMarkerColor(iChamberN); pChamber->Draw(); t.SetNDC();t.SetTextColor(iChamberN); t.DrawText(0.1,iChamberN*0.1,Form("Chamber %i",iChamberN)); }//chamber loop - gPad->GetView()->RotateView(94,45); +// gPad->GetView()->RotateView(94,45); }//void TestSeg() //__________________________________________________________________________________________________ void TestMenu() @@ -423,3 +423,92 @@ void ReadRaw() } fclose(pRawFile); } +//__________________________________________________________________________________________________ +void FillContribs(Int_t part,Int_t C, Bool_t print) +{ + static Int_t contrib[10][7][2]; + static Bool_t kEnter=kFALSE; + + C--; // chamber numbering from 1 to 7 + if(!kEnter) { + for(Int_t i=0;i<10;i++) { + for(Int_t j=0;j<7;j++) { + for(Int_t k=0;k<2;k++) contrib[i][j][k]=0; + } + } + kEnter=kTRUE; + } + + if(print) { + for(Int_t k=0;k<2;k++) {cout << "----------------Positives to RICH ---------------" << endl; + cout << " chamber 1 2 3 4 5 6 7 " << endl; + cout << " -------------------------------------------------" << endl; + for(Int_t i=0;i<4;i++) { + if(i==0) cout << " e"; + if(i==1) cout << "pi"; + if(i==2) cout << " K"; + if(i==3) cout << " p"; + if(k==0) cout << "+ "; + if(k==1) cout << "- "; + for(Int_t j=0;j<7;j++) { + cout << contrib[i][j][k] << " "; + } + cout << endl; + } + } + } + + // +ves + if(part==kPositron) contrib[0][C][0]++; + if(part==kPiPlus) contrib[1][C][0]++; + if(part==kKPlus) contrib[2][C][0]++; + if(part==kProton) contrib[3][C][0]++; + + // -ves + if(part==kElectron) contrib[0][C][1]++; + if(part==kPiMinus) contrib[1][C][1]++; + if(part==kKMinus) contrib[2][C][1]++; + if(part==kProtonBar) contrib[3][C][1]++; +} +//__________________________________________________________________________________________________ +void ParticleContribs() +{ + + TH1F *distH1 = new TH1F("distH1","distH1",100,0.,5.); + AliRICH *pRich = (AliRICH*)gAlice->GetDetector("RICH"); + Bool_t isHits =!pRich->GetLoader()->LoadHits(); + + pRich->GetLoader()->GetRunLoader()->LoadHeader(); pRich->GetLoader()->GetRunLoader()->LoadKinematics(); + + + if(!isHits){Error("Exec","No hits. No contribs to calculate.");return;} + + for(Int_t iEventN=0;iEventNGetEventsPerRun();iEventN++){//events Loop + pRich->GetLoader()->GetRunLoader()->GetEvent(iEventN); + cout << " event n. " << iEventN << endl; + Int_t nPrimaries = (Int_t)pRich->GetLoader()->TreeH()->GetEntries(); + TObjArray * Hits = new TObjArray[nPrimaries]; + + for(Int_t i=0;iGetLoader()->TreeH()->GetEntry(i); + Int_t nHits = pRich->Hits()->GetEntries(); + for(Int_t k=0;kHits()->At(k)); + + } +//deals with hits + for(Int_t i=0;iGetLoader()->TreeH()->GetEntry(i); + Int_t nHits = pRich->Hits()->GetEntries(); + for(Int_t j=0;jGetLoader()->GetRunLoader()->GetAliRun()->Stack()->Particle(pHit->GetTrack()); +// if(pParticle->P()>1) FillContribs(pParticle->GetPdgCode(),pHit->C(),kFALSE); + FillContribs(pParticle->GetPdgCode(),pHit->C(),kFALSE); + if(pParticle->GetPDG()->Charge()!=0) distH1->Fill(pHit->Length()); + }//hits loop + }//prims loop + }// event loop + FillContribs(0,0,kTRUE); + distH1->Draw(); + +}//ParticleContribs() -- 2.43.0