From e42a7b4622bb16919e0c2beb555bb52957968f22 Mon Sep 17 00:00:00 2001 From: kir Date: Mon, 31 May 2004 14:51:30 +0000 Subject: [PATCH] New convention on numbering and LRS --- RICH/AliRICH.cxx | 487 ++++++++++++++++++++++------------ RICH/AliRICH.h | 248 ++++++++--------- RICH/AliRICHChamber.cxx | 44 ++- RICH/AliRICHChamber.h | 77 ++---- RICH/AliRICHClusterFinder.cxx | 281 +++++++++----------- RICH/AliRICHClusterFinder.h | 43 ++- RICH/AliRICHDigitizer.cxx | 38 +-- RICH/AliRICHDigitizer.h | 4 +- RICH/AliRICHDisplFast.cxx | 97 ++----- RICH/AliRICHDisplFast.h | 8 +- RICH/AliRICHParam.cxx | 35 ++- RICH/AliRICHParam.h | 343 +++++++++++++----------- RICH/AliRICHv0.cxx | 58 +++- RICH/AliRICHv0.h | 13 +- RICH/AliRICHv1.cxx | 103 ++++++- RICH/AliRICHv1.h | 7 +- RICH/RICHLinkDef.h | 1 - 17 files changed, 1030 insertions(+), 857 deletions(-) diff --git a/RICH/AliRICH.cxx b/RICH/AliRICH.cxx index 2a8acbadda0..d57feb07398 100644 --- a/RICH/AliRICH.cxx +++ b/RICH/AliRICH.cxx @@ -16,6 +16,7 @@ #include "AliRICH.h" #include "AliRICHParam.h" #include "AliRICHChamber.h" +#include "AliRICHClusterFinder.h" #include #include #include @@ -29,7 +30,11 @@ #include #include #include +#include #include +#include +#include +#include ClassImp(AliRICHhit) //__________________________________________________________________________________________________ @@ -43,23 +48,27 @@ ClassImp(AliRICHdigit) //__________________________________________________________________________________________________ void AliRICHdigit::Print(Option_t*)const { - ::Info("digit","csxy=%6i, cfm=%9i, c=%2i, x=%3i, y=%3i, q=%8.3f, TID1=%5i, TID2=%5i, TID3=%5i", - Id(),fChFbMip,fChamber,fPadX,fPadY,fQdc,fTracks[0],fTracks[1],fTracks[2]); + ::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]); } //__________________________________________________________________________________________________ ClassImp(AliRICHcluster) //__________________________________________________________________________________________________ void AliRICHcluster::Print(Option_t*)const { - ::Info("cluster","CombiPid=%10i, c=%2i, size=%6i, dim=%5i, x=%7.3f, y=%7.3f, Q=%6i, st=%i", - fCombiPid,fChamber,fSize,fDimXY,fX,fY,fQdc,fStatus); -} -//__________________________________________________________________________________________________ -ClassImp(AliRICHreco) -//__________________________________________________________________________________________________ -void AliRICHreco::Print(Option_t*)const -{ - ::Info("reco","ThetaCherenkov=%9.6f, Nphotons=%4i, TID=%9i",fThetaCherenkov,fNphotons,fTid); + char *status=0; + switch(fStatus){ + case kRaw: status="raw" ;break; + case kResolved: status="resolved";break; + case kEmpty: status="empty" ;break; + } + if(fDigits) + ::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); + } //__________________________________________________________________________________________________ ClassImp(AliRICH) @@ -71,32 +80,23 @@ ClassImp(AliRICH) */ //END_HTML //__________________________________________________________________________________________________ -AliRICH::AliRICH() - :AliDetector() +AliRICH::AliRICH():AliDetector(),fpParam(0), fSdigits(0),fNsdigits(0),fDigitsNew(0),fClusters(0) { //Default ctor should not contain any new operators - fpParam =0; - fChambers =0; //AliDetector ctor deals with Hits and Digits - fSdigits =0; fNsdigits =0; - fDigitsNew =0; for(int i=0;iGetMCApp()->AddHitList(fHits); - fSdigits= 0; - fDigitsNew= 0; - fClusters= 0; - fRecos =0; + CreateHits(); gAlice->GetMCApp()->AddHitList(fHits); + fCounters.ResizeTo(20); fCounters.Zero(); if(GetDebug())Info("named ctor","Stop."); }//AliRICH::AliRICH(const char *name, const char *title) //__________________________________________________________________________________________________ @@ -106,59 +106,49 @@ AliRICH::~AliRICH() if(GetDebug()) Info("dtor","Start."); if(fpParam) delete fpParam; - if(fChambers) delete fChambers; if(fHits) delete fHits; if(fSdigits) delete fSdigits; if(fDigits) delete fDigits; if(fDigitsNew) {fDigitsNew->Delete(); delete fDigitsNew;} if(fClusters) {fClusters->Delete(); delete fClusters;} - if(fRecos) delete fRecos; if(GetDebug()) Info("dtor","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."); - - AliLoader * richLoader = GetLoader(); - AliRunLoader * runLoader = GetLoader()->GetRunLoader(); - for(Int_t iEventN=0;iEventNGetRunLoader()->GetAliRun()->GetEventsPerRun();iEventN++){//events loop - runLoader->GetEvent(iEventN); + GetLoader()->GetRunLoader()->GetEvent(iEventN);//get next event - if (!richLoader->TreeH()) richLoader->LoadHits(); - if (!runLoader->TreeE()) runLoader->LoadHeader(); - if (!runLoader->TreeK()) runLoader->LoadKinematics();//from - if (!richLoader->TreeS()) richLoader->MakeTree("S"); MakeBranch("S");//to + if(!GetLoader()->TreeH()) GetLoader()->LoadHits(); GetLoader()->GetRunLoader()->LoadHeader(); + GetLoader()->GetRunLoader()->LoadKinematics();//from + if(!GetLoader()->TreeS()) GetLoader()->MakeTree("S"); MakeBranch("S");//to for(Int_t iPrimN=0;iPrimNTreeH()->GetEntries();iPrimN++){//prims loop - richLoader->TreeH()->GetEntry(iPrimN); + GetLoader()->TreeH()->GetEntry(iPrimN); for(Int_t iHitN=0;iHitNGetEntries();iHitN++){//hits loop - AliRICHhit *pHit=(AliRICHhit*)Hits()->At(iHitN); - TVector2 x2 = P()->ShiftToWirePos(C(pHit->C())->Glob2Loc(pHit->OutX3())); - Int_t iTotQdc=P()->TotQdc(x2,pHit->Eloss()); + AliRICHhit *pHit=(AliRICHhit*)Hits()->At(iHitN);//get current hit + TVector2 x2 = C(pHit->C())->Glob2Loc(pHit->OutX3());//hit position in the chamber local system + Int_t iTotQdc=P()->TotQdc(x2,pHit->Eloss());//total charge produced by hit, 0 if hit in dead zone if(iTotQdc==0) continue; - Int_t iPadXmin,iPadXmax,iPadYmin,iPadYmax; - P()->Loc2Area(x2,iPadXmin,iPadYmin,iPadXmax,iPadYmax);//determine affected pads - if(GetDebug()) Info("Hits2SDigits","left-down=(%i,%i) right-up=(%i,%i)",iPadXmin,iPadYmin,iPadXmax,iPadYmax); - for(Int_t iPadY=iPadYmin;iPadY<=iPadYmax;iPadY++)//affected pads loop - for(Int_t iPadX=iPadXmin;iPadX<=iPadXmax;iPadX++){ - Double_t padQdc=iTotQdc*P()->FracQdc(x2,iPadX,iPadY); - if(padQdc>0.1) AddSDigit(pHit->C(),iPadX,iPadY,padQdc, - runLoader->Stack()->Particle(pHit->GetTrack())->GetPdgCode(),pHit->GetTrack()); + 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); + 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); + if(padQdc>0.1) AddSDigit(pHit->C(),pad,padQdc,GetLoader()->GetRunLoader()->Stack()->Particle(pHit->GetTrack())->GetPdgCode(),pHit->GetTrack()); }//affected pads loop }//hits loop }//prims loop - richLoader->TreeS()->Fill(); - richLoader->WriteSDigits("OVERWRITE"); + GetLoader()->TreeS()->Fill(); + GetLoader()->WriteSDigits("OVERWRITE"); ResetSDigits(); }//events loop - richLoader->UnloadHits(); - runLoader->UnloadHeader(); - runLoader->UnloadKinematics(); + GetLoader()->UnloadHits(); GetLoader()->GetRunLoader()->UnloadHeader(); GetLoader()->GetRunLoader()->UnloadKinematics(); GetLoader()->UnloadSDigits(); if(GetDebug()) Info("Hit2SDigits","Stop."); }//Hits2SDigits() @@ -177,7 +167,7 @@ void AliRICH::BuildGeometry() Float_t len=P()->SectorSizeY(); new TBRIK("PHOTO","PHOTO","void",wid/2,0.1,len/2); - for(int i=1;i<=kNCH;i++){ + 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->SetLineColor(kRed); @@ -208,18 +198,17 @@ void AliRICH::BuildGeometry() //______________________________________________________________________________ void AliRICH::CreateMaterials() { - // - // *** DEFINITION OF AVAILABLE RICH MATERIALS *** - +// Definition of available RICH materials #include "Opticals.h" Float_t a=0,z=0,den=0,radl=0,absl=0; 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( 1, "Air $",a=14.61,z=7.3, den=0.001205,radl=30420.0,absl=67500);//(Air) - AliMedium(1, "DEFAULT MEDIUM AIR$", 1, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); + 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); 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); @@ -230,11 +219,11 @@ void AliRICH::CreateMaterials() 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); - AliMaterial(50, "ALUM", a=26.98,z=13.0,den=2.7, radl=8.9, absl=0); //aluminium sheet (Al) + 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(31, "COPPER$", a=63.54,z=29.0,den=8.96, radl=1.4, absl=0); //(Cu) - AliMedium(12, "PCB_COPPER", 31, 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); 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 @@ -244,21 +233,30 @@ void AliRICH::CreateMaterials() AliMedium(8, "QUARTZO$", 21, 1, 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 (30, "FRE",aFreon,zFreon,den=1.7,-2,wmatFreon);//Freon (C6F14) - AliMedium(4, "FREON$", 30, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); + 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); Float_t aMethane[2]={12.01,1}; Float_t zMethane[2]={6,1}; Float_t wmatMethane[2]={1,4}; - AliMixture (40, "MET", aMethane, zMethane, den=7.17e-4,-2, wmatMethane);//methane (CH4) - AliMedium(5, "METHANE$", 40, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); + 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 + AliMedium(kSteel, "STEEL", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); - AliMixture (41, "METG", aMethane, zMethane, den=7.17e-4, -2, wmatMethane); - AliMedium(kGAP, "GAP$", 41, 1, isxfld, sxmgmx,tmaxfd, 0.1, -deemax, epsil, -stmin); + AliMaterial(material=46, "PERPEX$", a=63.54,z=29.0,den=8.96, radl=1.4, absl=0); //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 + 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 (32, "GLASS",aGlass, zGlass, den=1.74, 5, wGlass);//Glass 50%C+10.5%Si+35.5%O+3% + 1% - AliMedium(11, "GLASS", 32, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); + 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); @@ -387,14 +385,14 @@ void AliRICH::MakeBranch(Option_t* option) if(cD&&fLoader->TreeD()){//D CreateDigits(); - for(Int_t i=0;iTreeD(),Form("%s%d",GetName(),i+1),&((*fDigitsNew)[i]),kBufferSize,0); } }//D if(cR&&fLoader->TreeR()){//R CreateClusters(); - for(Int_t i=0;iTreeR(),Form("%sClusters%d",GetName(),i+1), &((*fClusters)[i]), kBufferSize, 0); }//R if(GetDebug())Info("MakeBranch","Stop."); @@ -420,7 +418,7 @@ void AliRICH::SetTreeAddress() if(fLoader->TreeD()){//D if(GetDebug())Info("SetTreeAddress","tree D is requested."); - for(int i=0;iTreeD()->GetBranch(Form("%s%d",GetName(),i+1)); if(branch){CreateDigits(); branch->SetAddress(&((*fDigitsNew)[i]));} } @@ -428,7 +426,7 @@ void AliRICH::SetTreeAddress() if(fLoader->TreeR()){//R if(GetDebug())Info("SetTreeAddress","tree R is requested."); - for(int i=0;iTreeR()->GetBranch(Form("%sClusters%d" ,GetName(),i+1)); if(branch){CreateClusters(); branch->SetAddress(&((*fClusters)[i]));} } @@ -440,14 +438,14 @@ void AliRICH::Print(Option_t *option)const { //Debug printout TObject::Print(option); - P()->Dump(); - fChambers->Print(option); + P()->Print(); + fCounters.Print(); }//void AliRICH::Print(Option_t *option)const //__________________________________________________________________________________________________ void AliRICH::CreateGeometry() { //Creates detailed geometry simulation (currently GEANT volumes tree) - if(GetDebug())Info("CreateGeometry","Start."); + if(GetDebug())Info("CreateGeometry","Start main."); //Opaque quartz thickness Float_t oquaThickness = .5; //CsI dimensions @@ -460,19 +458,16 @@ void AliRICH::CreateGeometry() Float_t zs; Int_t idrotm[1099]; Float_t par[3]; - -//External aluminium box - par[0]=68.8;par[1]=13;par[2]=70.86; gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3); -//Air - par[0]=66.3; par[1] = 13; par[2] = 68.35; gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], 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); +// 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); +// 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); -//Aluminium sheet - par[0]=66.3;par[1]=0.025;par[2]=68.35; gMC->Gsvolu("ALUM", "BOX ", idtmed[1009], 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; @@ -502,7 +497,7 @@ void AliRICH::CreateGeometry() //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()->ProximityGap()/2;par[2]=pcY/2;gMC->Gsvolu("GAP ","BOX ",(*fIdtmed)[kGAP],par,3); + 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 @@ -561,10 +556,10 @@ void AliRICH::CreateGeometry() 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("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"); @@ -621,11 +616,11 @@ void AliRICH::CreateGeometry() 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()->ProximityGap()/2 - 0.0001, 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()->ProximityGap()/2 - .1, 0., 0, "ONLY"); + 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 @@ -638,100 +633,244 @@ void AliRICH::CreateGeometry() //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"); - + //place chambers into mother volume ALIC - for(int i=1;i<=kNCH;i++){ + 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"); } - if(GetDebug())Info("CreateGeometry","Stop."); + + if(GetDebug())Info("CreateGeometry","Stop main."); }//void AliRICH::CreateGeometry() //__________________________________________________________________________________________________ -void AliRICH::CreateChambers() +void AliRICH::Reconstruct()const { -//create all RICH Chambers on each call. Previous chambers deleted - if(fChambers) delete fChambers; - if(GetDebug())Info("CreateChambers","Creating RICH chambers."); - fChambers=new TObjArray(kNCH); - fChambers->SetOwner(); - for(int i=0;iAddAt(new AliRICHChamber(i+1,P()),i); -}//void AliRICH::CreateChambers() + AliRICHClusterFinder finder(const_cast(this)); + finder.Exec(); +} //__________________________________________________________________________________________________ -void AliRICH::GenerateFeedbacks(Int_t iChamber,Float_t eloss) -{ -// Generate FeedBack photons -// eloss=0 means photon so only pulse height distribution is to be analysed. This one is done in AliRICHParam::TotQdc() - - TLorentzVector x4; - gMC->TrackPosition(x4); - TVector2 x2=C(iChamber)->Glob2Loc(x4); - Int_t sector=P()->Sector(x2); 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 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 - for(Int_t i=0;iGetRandom()->RndmArray(2,ranf); //Sample direction - cthf=ranf[0]*2-1.0; - if(cthf<0) continue; - sthf = TMath::Sqrt((1 - cthf) * (1 + cthf)); - phif = ranf[1] * 2 * TMath::Pi(); +void AliRICH::ControlPlots() +{ +//Creates a set of hists to control the results of simulation + + TH1F *pHxD=0,*pHyD=0,*pNumClusH1=0, + *pQdcH1=0, *pSizeH1=0, + *pPureMipQdcH1=0,*pPureMipSizeH1=0, + *pPureCerQdcH1=0,*pPureCerSizeH1=0, + *pPureFeeQdcH1=0,*pPureFeeSizeH1=0, + *pMipQdcH1=0, *pPhotQdcH1=0; + TH2F *pMapH2=0,*pPureMipMapH2=0,*pPureCerMapH2=0,*pPureFeeMapH2=0; + + Bool_t isDig =!GetLoader()->LoadDigits(); + Bool_t isClus=!GetLoader()->LoadRecPoints(); + + if(!isDig && !isClus){Error("ControlPlots","No digits and clusters! Nothing to do.");return;} + + TStopwatch sw;TDatime time; - if(Double_t randomNumber=gMC->GetRandom()->Rndm()<=0.57) - enfp = 7.5e-9; - else if(randomNumber<=0.7) - enfp = 6.4e-9; - else - enfp = 7.9e-9; + TFile *pFile = new TFile("$(HOME)/RCP.root","RECREATE"); + + if(isDig){ + cout<<"Digits available\n"; + 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 + + if(isClus){ + cout<<"Clusters available\n"; + pNumClusH1=new TH1F("NumClusPerEvent","Number of clusters per event;number",50,0,49); - - dir[0] = sthf * TMath::Sin(phif); dir[1] = cthf; dir[2] = sthf * TMath::Cos(phif); - gMC->Gdtom(dir, mom, 2); - mom[0]*=enfp; mom[1]*=enfp; mom[2]*=enfp; - mom[3] = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2]); + pQdcH1 =new TH1F("ClusQdc", "Cluster Charge all chambers;q [QDC]",R()->P()->MaxQdc(),0,R()->P()->MaxQdc()); + pSizeH1 =new TH1F("ClusSize", "Cluster size all chambers;size [number of pads in cluster]",100,0,100); + pMapH2 =new TH2F("ClusMap", "Cluster map;x [cm];y [cm]",1000,0,R()->P()->PcSizeX(),1000,0,R()->P()->PcSizeY()); + + pMipQdcH1 =new TH1F("QdcMip" ,"MIP Cluster Charge all chambers;q [QDC]",R()->P()->MaxQdc(),0,R()->P()->MaxQdc()); + pPhotQdcH1 =new TH1F("QdcPhot" ,"Cer+Fee Cluster Charge all chambers;q [QDC]",R()->P()->MaxQdc(),0,R()->P()->MaxQdc()); + + pPureMipQdcH1 =new TH1F("QdcPureMip" ,"MIP only Cluster Charge all chambers;q [QDC]",R()->P()->MaxQdc(),0,R()->P()->MaxQdc()); + pPureMipSizeH1=new TH1F("SizePureMip" ,"MIP only Cluster size all chambers;size [number of pads in cluster]",100,0,100); + pPureMipMapH2 =new TH2F("MapPureMip" ,"MIP only Cluster map;x [cm];y [cm]",1000,0,R()->P()->PcSizeX(),1000,0,R()->P()->PcSizeY()); + + pPureCerQdcH1 =new TH1F("QdcPureCer" ,"Cerenkov only Cluster Charge all chambers;q [QDC]",R()->P()->MaxQdc(),0,R()->P()->MaxQdc()); + pPureCerSizeH1=new TH1F("SizePureCer" ,"Cernekov only Cluster size all chambers;size [number of pads in cluster]",100,0,100); + pPureCerMapH2 =new TH2F("MapPureCer" ,"Cerenkov only Cluster map;x [cm];y [cm]",1000,0,R()->P()->PcSizeX(),1000,0,R()->P()->PcSizeY()); - // Polarisation - e1[0]= 0; e1[1]=-dir[2]; e1[2]= dir[1]; - e2[0]=-dir[1]; e2[1]= dir[0]; e2[2]= 0; - e3[0]= dir[1]; e3[1]= 0; e3[2]=-dir[0]; + pPureFeeQdcH1 =new TH1F("QdcPureFee" ,"Feedback only Cluster Charge all chambers;q [QDC]",R()->P()->MaxQdc(),0,R()->P()->MaxQdc()); + pPureFeeSizeH1=new TH1F("SizePureFee" ,"Feedback only Cluster size all chambers;size [number of pads in cluster]",100,0,100); + pPureFeeMapH2 =new TH2F("MapPureFee" ,"Feedback only Cluster map;x [cm];y [cm]",1000,0,R()->P()->PcSizeX(),1000,0,R()->P()->PcSizeY()); + }//isClus + + for(Int_t iEvtN=0;iEvtN < GetLoader()->GetRunLoader()->GetAliRun()->GetEventsPerRun();iEvtN++){//events loop + GetLoader()->GetRunLoader()->GetEvent(iEvtN); //gets current event - vmod=0; - for(j=0;j<3;j++) vmod+=e1[j]*e1[j]; - if (!vmod) for(j=0;j<3;j++) { - uswop=e1[j]; - e1[j]=e3[j]; - e3[j]=uswop; - } - vmod=0; - for(j=0;j<3;j++) vmod+=e2[j]*e2[j]; - if (!vmod) for(j=0;j<3;j++) { - uswop=e2[j]; - e2[j]=e3[j]; - e3[j]=uswop; + if(!GetLoader()->TreeH()) GetLoader()->LoadHits(); + for(Int_t iPrimN=0;iPrimN < GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop + GetLoader()->TreeH()->GetEntry(iPrimN); } - vmod=0; for(j=0;j<3;j++) vmod+=e1[j]*e1[j]; vmod=TMath::Sqrt(1/vmod); for(j=0;j<3;j++) e1[j]*=vmod; - vmod=0; for(j=0;j<3;j++) vmod+=e2[j]*e2[j]; vmod=TMath::Sqrt(1/vmod); for(j=0;j<3;j++) e2[j]*=vmod; + if(isClus) GetLoader()->TreeR()->GetEntry(0); + if(isDig) GetLoader()->TreeD()->GetEntry(0); - phi = gMC->GetRandom()->Rndm()* 2 * TMath::Pi(); - for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi); - gMC->Gdtom(pol, pol, 2); - Int_t outputNtracksStored; - gAlice->GetMCApp()->PushTrack(1, //do not transport - gAlice->GetMCApp()->GetCurrentTrackNumber(),//parent track - kFeedback, //PID - mom[0],mom[1],mom[2],mom[3], //track momentum - x4.X(),x4.Y(),x4.Z(),x4.T(), //track origin - pol[0],pol[1],pol[2], //polarization - kPFeedBackPhoton, - outputNtracksStored, - 1.0); - }//feedbacks loop -}//GenerateFeedbacks() + for(Int_t iChamN=1;iChamN<=7;iChamN++){//chambers loop + if(isClus){ + Int_t iNclusCham=Clusters(iChamN)->GetEntries(); if(iNclusCham) pNumClusH1->Fill(iNclusCham);//number of clusters per event + for(Int_t iClusN=0;iClusNAt(iClusN); + pQdcH1 ->Fill(pClus->Q()); + pSizeH1 ->Fill(pClus->Size()); + pMapH2 ->Fill(pClus->X(),pClus->Y()); //common + + if(pClus->IsSingleMip()) {pPureMipQdcH1 ->Fill(pClus->Q()); + pPureMipSizeH1->Fill(pClus->Size()); + pPureMipMapH2 ->Fill(pClus->X(),pClus->Y());}//Pure Mips + + if(pClus->IsSingleCerenkov()){pPureCerQdcH1 ->Fill(pClus->Q()); + pPureCerSizeH1->Fill(pClus->Size()); + pPureCerMapH2 ->Fill(pClus->X(),pClus->Y());}//Pure Cerenkovs + + if(pClus->IsSingleFeedback()){pPureFeeQdcH1 ->Fill(pClus->Q()); + pPureFeeSizeH1->Fill(pClus->Size()); + pPureFeeMapH2 ->Fill(pClus->X(),pClus->Y());}//Pure Feedbacks + + if(pClus->IsMip()) {pMipQdcH1 ->Fill(pClus->Q());} //MIP+ other contributions + if(!pClus->IsPureMip()) pPhotQdcH1->Fill(pClus->Q()); //not MIP + }//clusters loop + }//isClus + 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()); + pHxD->Fill(hitV2.X()-digV2.X()); pHyD->Fill(hitV2.Y()-digV2.Y()); + }//digits loop + }//isDig + }//chambers loop + Info("ControlPlots","Event %i processed.",iEvtN); + }//events loop + + if(isDig) GetLoader()->UnloadDigits(); + if(isClus) GetLoader()->UnloadRecPoints(); + + pFile->Write(); delete pFile; + sw.Print();time.Print(); +}//ControlPlots() +//__________________________________________________________________________________________________ +AliRICHhit* AliRICH::Hit(Int_t tid) +{ +//defines which hit provided by given tid for the currently loaded event + R()->GetLoader()->LoadHits(); + for(Int_t iPrimN=0;iPrimNGetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop + R()->GetLoader()->TreeH()->GetEntry(iPrimN); + for(Int_t iHitN=0;iHitNHits()->GetEntries();iHitN++){ + AliRICHhit *pHit=(AliRICHhit*)R()->Hits()->At(iHitN); + if(tid==pHit->Track()) {R()->GetLoader()->UnloadHits();return pHit;} + }//hits + }//prims loop + R()->GetLoader()->UnloadHits(); + return 0; +} +//__________________________________________________________________________________________________ +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); + R()->GetLoader()->GetRunLoader()->GetEvent(iEvtN); + if(R()->GetLoader()->LoadHits()) return; + + Int_t iTotalHits=0; + for(Int_t iPrimN=0;iPrimNGetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop + R()->GetLoader()->TreeH()->GetEntry(iPrimN); + R()->Hits()->Print(); + iTotalHits+=R()->Hits()->GetEntries(); + } + R()->GetLoader()->UnloadHits(); + Info("PrintHits","totally %i hits",iTotalHits); +} +//__________________________________________________________________________________________________ +void AliRICH::PrintSDigits(Int_t iEvtN) +{ +//prints a list of RICH sdigits for a given event + Info("PrintSDigits","List of RICH sdigits for event %i",iEvtN); + R()->GetLoader()->GetRunLoader()->GetEvent(iEvtN); + if(R()->GetLoader()->LoadSDigits()) return; + + R()->GetLoader()->TreeS()->GetEntry(0); + R()->SDigits()->Print(); + R()->GetLoader()->UnloadSDigits(); + Info("PrintSDigits","totally %i sdigits",R()->SDigits()->GetEntries()); +} +//__________________________________________________________________________________________________ +void AliRICH::PrintDigits(Int_t iEvtN) +{ +//prints a list of RICH digits for a given event + Info("PrintDigits","List of RICH digits for event %i",iEvtN); + R()->GetLoader()->GetRunLoader()->GetEvent(iEvtN); + if(R()->GetLoader()->LoadDigits()) return; + + Int_t iTotalDigits=0; + R()->GetLoader()->TreeD()->GetEntry(0); + for(Int_t iChamber=1;iChamber<=kNchambers;iChamber++){ + R()->Digits(iChamber)->Print(); + iTotalDigits+=R()->Digits(iChamber)->GetEntries(); + } + R()->GetLoader()->UnloadDigits(); + Info("PrintDigits","totally %i Digits",iTotalDigits); +} +//__________________________________________________________________________________________________ +void AliRICH::PrintClusters(Int_t iEvtN) +{ +//prints a list of RICH clusters for a given event + Info("PrintClusters","List of RICH clusters for event %i",iEvtN); + R()->GetLoader()->GetRunLoader()->GetEvent(iEvtN); + if(R()->GetLoader()->LoadRecPoints()) return; + + Int_t iTotalClusters=0; + R()->GetLoader()->TreeR()->GetEntry(0); + for(Int_t iChamber=1;iChamber<=kNchambers;iChamber++){ + R()->Clusters(iChamber)->Print(); + iTotalClusters+=R()->Clusters(iChamber)->GetEntries(); + } + R()->GetLoader()->UnloadRecPoints(); + Info("PrintClusters","totally %i clusters",iTotalClusters); +} +//__________________________________________________________________________________________________ +void AliRICH::FillESD(AliESD *pESD)const +{ +//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) {pTrack->ResetRICH(); 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 +} diff --git a/RICH/AliRICH.h b/RICH/AliRICH.h index 88a51e89090..fa06ed9ad04 100644 --- a/RICH/AliRICH.h +++ b/RICH/AliRICH.h @@ -18,7 +18,7 @@ class AliRICHhit : public AliHit { public: - AliRICHhit():AliHit() {fChamber=kBad; fEloss=kBad; fInX3.SetXYZ(0,0,0);fOutX3.SetXYZ(0,0,0);} + AliRICHhit():AliHit(),fChamber(kBad),fEloss(kBad) {fInX3.SetXYZ(0,0,0);fOutX3.SetXYZ(0,0,0);} AliRICHhit(Int_t c,Int_t tid,TVector3 in,TVector3 out,Double_t e):AliHit(0,tid) {fChamber=c;fInX3=in; fOutX3=out;fEloss=e; fX=out.X();fY=out.Y();fZ=out.Z();} virtual ~AliRICHhit() {;} @@ -42,26 +42,26 @@ protected: class AliRICHdigit :public AliDigit { public: - AliRICHdigit() {fChFbMip=fChamber=fPadX=fPadY=fTracks[0]=fTracks[1]=fTracks[2]=kBad;fQdc=kBad;} - AliRICHdigit(Int_t c,Int_t x,Int_t y,Double_t q,Int_t cpid,Int_t tid0,Int_t tid1,Int_t tid2) - {fPadX=x;fPadY=y;fQdc=q;fChamber=10*c+AliRICHParam::Sector(x,y);fChFbMip=cpid;fTracks[0]=tid0;fTracks[1]=tid1;fTracks[2]=tid2;} + AliRICHdigit():AliDigit(),fCFM(0),fChamber(0),fPadX(0),fPadY(0),fQdc(kBad){fTracks[0]=fTracks[1]=fTracks[2]=kBad;} + AliRICHdigit(Int_t c,TVector pad,Double_t q,Int_t cfm,Int_t tid0,Int_t tid1,Int_t tid2) + {fPadX=(Int_t)pad[0];fPadY=(Int_t)pad[1];fQdc=q;fChamber=10*c+AliRICHParam::Pad2Sec(pad);fCFM=cfm;fTracks[0]=tid0;fTracks[1]=tid1;fTracks[2]=tid2;} virtual ~AliRICHdigit() {;} Int_t Compare(const TObject *pObj) const {if(Id()==((AliRICHdigit*)pObj)->Id())return 0;else if(Id()>((AliRICHdigit*)pObj)->Id())return 1;else return -1;} //virtual - Bool_t IsSortable() const{return kTRUE;} //virtual - void Print(Option_t *option="") const; //virtual - Int_t ChFbMi() const{return fChFbMip;} //particle mixture for this digit + virtual Bool_t IsSortable() const{return kTRUE;} //sort interface + virtual void Print(Option_t *option="") const; //virtual + Int_t ChFbMi() const{return fCFM;} //particle mixture for this digit Int_t C() const{return fChamber/10;} //chamber number Int_t S() const{return fChamber-(fChamber/10)*10;} //sector number Int_t X() const{return fPadX;} //x position of the pad Int_t Y() const{return fPadY;} //y postion of the pad + TVector Pad() const{Float_t v[2]={fPadX,fPadY}; return TVector(2,v);} Int_t Id() const{return fChamber*10000000+fPadX*1000+fPadY;} //absolute id of this pad Double_t Q() const{return fQdc;} //charge in terms of ADC channels - Int_t Tid(Int_t i) const{return fTracks[i];} //track reference produced this digit void AddTidOffset(Int_t offset) {for (Int_t i=0; i<3; i++) if (fTracks[i]>0) fTracks[i]+=offset;}; protected: - Int_t fChFbMip; //1000000*Ncerenkovs+1000*Nfeedbacks+Nmips + Int_t fCFM; //1000000*Ncerenkovs+1000*Nfeedbacks+Nmips Int_t fChamber; //10*chamber number+ sector number Int_t fPadX; //pad number along X Int_t fPadY; //pad number along Y @@ -73,42 +73,48 @@ protected: class AliRICHcluster :public TObject { public: - enum ClusterStatus {kEdge,kShape,kSize,kRaw,kResolved}; - AliRICHcluster() {fSize=fQdc=fStatus=fChamber=fDimXY=0;fX=fY=kBad;fDigits=0;} - virtual ~AliRICHcluster() {delete fDigits;} - AliRICHcluster& operator=(const AliRICHcluster&) {return *this;} - Int_t Nlocals() const{return fSize - 10000*(fSize/10000);} // - Int_t Size() const{return fSize/10000;} // - Int_t Fsize() const{return fSize;} // - Int_t DimXY() const{return fDimXY;} // - Int_t C() const{return fChamber/10;} // - Int_t S() const{return fChamber-(fChamber/10)*10;} // - Int_t Fchamber() const{return fChamber;} // - Int_t Q() const{return fQdc;} // - Double_t X() const{return fX;} // - Double_t Y() const{return fY;} // - Int_t Status() const{return fStatus;} // - void SetStatus(Int_t status) {fStatus=status;} // - Int_t Nmips() const{return fCombiPid-1000000*Ncerenkovs()-1000*Nfeedbacks();} // - Int_t Ncerenkovs() const{return fCombiPid/1000000;} // - Int_t Nfeedbacks() const{return (fCombiPid-1000000*Ncerenkovs())/1000;} // - Bool_t IsPureMip() const{return fCombiPid<1000;} - Bool_t IsPureCerenkov() const{return Nmips()==0&&Nfeedbacks()==0;} // - Bool_t IsPureFeedback() const{return Nmips()==0&&Ncerenkovs()==0;} // - Int_t CombiPid() const{return fCombiPid;} // - void SetCombiPid(Int_t ckov,Int_t feeds,Int_t mips) {fCombiPid=1000000*ckov+1000*feeds+mips;} // - void Fill(AliRICHcluster *pRaw,Double_t x,Double_t y, Double_t q, Int_t combipid) - {fCombiPid=combipid;fChamber=pRaw->Fchamber();fSize=pRaw->Fsize(); - fQdc=(Int_t)(q*pRaw->Q());fX=x;fY=y;fStatus=kResolved;} // - TObjArray* Digits() const{return fDigits;} // - void Print(Option_t *option="")const; //virtual - inline void AddDigit(AliRICHdigit *pDig); // - inline void CoG(Int_t nLocals); // - void Reset() {fSize=fQdc=fStatus=fChamber=fDimXY=kBad;fX=fY=kBad;delete fDigits;fDigits=0;} // + 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();} + 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 + 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 + Int_t Fsize() const{return fSize;} // + Int_t Shape() const{return fShape;} //cluster shape rectangulare + Int_t C() const{return fChamber/10;} //chamber number + Int_t S() const{return fChamber-(fChamber/10)*10;} //sector number + Int_t Fchamber() const{return fChamber;} // + Int_t Q() const{return fQdc;} //cluster charge in QDC channels + Double_t X() const{return fX;} //cluster x position in LRS + Double_t Y() const{return fY;} //cluster y position in LRS + Int_t Status() const{return fStatus;} // + void SetStatus(Int_t status) {fStatus=status;} // + Int_t Nmips() const{return fCFM-1000000*Ncerenkovs()-1000*Nfeedbacks();} // + Int_t Ncerenkovs() const{return fCFM/1000000;} // + Int_t Nfeedbacks() const{return (fCFM-1000000*Ncerenkovs())/1000;} // + Bool_t IsPureMip() const{return fCFM<1000;} // + Bool_t IsPureCerenkov() const{return Nmips()==0&&Nfeedbacks()==0;} // + Bool_t IsPureFeedback() const{return Nmips()==0&&Ncerenkovs()==0;} // + Bool_t IsSingleMip() const{return Nmips()==1&&Ncerenkovs()==0&&Nfeedbacks()==0;} // + Bool_t IsSingleCerenkov() const{return Nmips()==0&&Ncerenkovs()==1&&Nfeedbacks()==0;} // + Bool_t IsSingleFeedback() const{return Nmips()==0&&Ncerenkovs()==0&&Nfeedbacks()==1;} // + Bool_t IsMip() const{return Nmips()!=0;} // + Bool_t IsCerenkov() const{return Ncerenkovs()!=0;} // + Bool_t IsFeedback() const{return Nfeedbacks()!=0;} // + 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;} // protected: - Int_t fCombiPid; //1000000*Ncerenkovs+1000*Nfeedbacks+Nmips + Int_t fCFM; //1000000*Ncerenkovs+1000*Nfeedbacks+Nmips Int_t fSize; //10000*(how many digits belong to this cluster) + nLocalMaxima - Int_t fDimXY; //100*xdim+ydim box containing the cluster + Int_t fShape; //100*xdim+ydim box containing the cluster Int_t fQdc; //QDC value Int_t fChamber; //10*module number+sector number Double_t fX; //local x postion @@ -119,53 +125,35 @@ protected: };//class AliRICHcluster //__________________________________________________________________________________________________ void AliRICHcluster::AddDigit(AliRICHdigit *pDig) -{// - if(!fDigits) {fQdc=fSize=fCombiPid=0;fDigits = new TObjArray;} +{ +// Adds a given digit to the list of digits belonging to this cluster + if(!fDigits) {fQdc=fSize=fCFM=0;fDigits = new TObjArray;} fQdc+=(Int_t)pDig->Q(); fDigits->Add(pDig); fChamber=10*pDig->C()+pDig->S(); fSize+=10000; + fStatus=kRaw; } //__________________________________________________________________________________________________ void AliRICHcluster::CoG(Int_t nLocals) -{// - Int_t xmin=999,ymin=999,xmax=0,ymax=0; +{ +// Calculates naive cluster position as a center of gravity of its digits. + Float_t xmin=999,ymin=999,xmax=0,ymax=0; fX=fY=0; for(Int_t iDig=0;iDigAt(iDig); - Int_t padX = pDig->X();Int_t padY = pDig->Y();Double_t q=pDig->Q(); - TVector2 x2=AliRICHParam::Pad2Loc(padX,padY); + TVector pad=pDig->Pad(); Double_t q=pDig->Q(); + TVector2 x2=AliRICHParam::Pad2Loc(pad); fX += x2.X()*q;fY +=x2.Y()*q; - if(padXxmax)xmax=padX;if(padYymax)ymax=padY; + if(pad[0]xmax)xmax=pad[0];if(pad[1]ymax)ymax=pad[1]; } fX/=fQdc;fY/=fQdc;//Center of Gravity - fDimXY = 100*(xmax-xmin+1)+ymax-ymin+1;//find box containing cluster + fShape=Int_t(100*(xmax-xmin+1)+ymax-ymin+1);//find box containing cluster fSize+=nLocals; fStatus=kRaw; }//CoG() - -class AliRICHreco: public TObject -{ -public: - AliRICHreco() {fTid=fNphotons=kBad; fThetaCherenkov=kBad;} - AliRICHreco(Int_t tid,Double_t thetaCherenkov,Int_t nPhotons) {fTid=tid;fThetaCherenkov=thetaCherenkov;fNphotons=nPhotons;} - - virtual ~AliRICHreco() {;} - - void Print(Option_t *option="")const; //virtual print - -protected: - Int_t fTid; // track Id reference - Int_t fNphotons; // number of photons contributed to the recontruction - Double_t fThetaCherenkov; // reconstructed Theta Cerenkov for a given charged track - - ClassDef(AliRICHreco,1) //RICH reco class - -};//class AliRICHreco - //__________________AliRICH_________________________________________________________________________ -class AliRICHParam; -class AliRICHChamber; +class AliESD; class AliRICH : public AliDetector { @@ -177,69 +165,66 @@ public: AliRICH& operator=(const AliRICH&) {return *this;} //framework part - virtual Int_t IsVersion() const =0; //virtual - virtual void StepManager() =0; //virtual - void Hits2SDigits(); //virtual - AliDigitizer* CreateDigitizer(AliRunDigitizer* man) const {return new AliRICHDigitizer(man);} //virtual - void Print(Option_t *option)const; //virtual - void SetTreeAddress(); //virtual - void MakeBranch(Option_t *opt=" "); //virtual - void CreateMaterials(); //virtual - virtual void BuildGeometry(); //virtual - virtual void CreateGeometry(); //virtual + virtual Int_t IsVersion() const =0; //interface from + 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 SetTreeAddress(); //interface from AliLoader + virtual void MakeBranch(Option_t *opt=" "); //interface from AliLoader + virtual void CreateMaterials(); //interface from AliMC + virtual void CreateGeometry(); //interface from AliMC + virtual void BuildGeometry(); //interface //private part Float_t AbsoCH4(Float_t x)const; //calculates absorption length for methane Float_t Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)const; //deals with Fresnel absorption - void GenerateFeedbacks(Int_t iChamber,Float_t eloss=0); //generates feedback photons; eloss=0 for photon - void CreateChambers(); //creates set of chambers inline void CreateHits(); //create hits container as a simple list inline void CreateSDigits(); //create sdigits container as a simple list inline void CreateDigits(); //create digits container as 7 lists, one per chamber inline void CreateClusters(); //create clusters container as 7 lists, one per chamber - inline void CreateRecos(); //create recos container // void ResetHits() {AliDetector::ResetHits();} //virtual void ResetSDigits() {fNsdigits=0; if(fSdigits) fSdigits ->Clear();} - void ResetDigits() {if(fDigitsNew)for(int i=0;iAt(i)->Clear();fNdigitsNew[i]=0;}} //virtual - void ResetClusters() {if(fClusters) for(int i=0;iAt(i)->Clear();fNclusters[i]=0;}} - void ResetRecos() {if(fRecos) fRecos->Clear();fNrecos=0;} + void ResetDigits() {if(fDigitsNew)for(int i=0;iAt(i)->Clear();fNdigitsNew[i]=0;}} //virtual + void ResetClusters() {if(fClusters) for(int i=0;iAt(i)->Clear();fNclusters[i]=0;}} TClonesArray* SDigits() const{return fSdigits;} TClonesArray* Digits(Int_t iC) const{if(fDigitsNew) return (TClonesArray *)fDigitsNew->At(iC-1);else return 0;} TClonesArray* Clusters(Int_t iC) const{if(fClusters) return (TClonesArray *)fClusters->At(iC-1);else return 0;} - TClonesArray* Recos() const{return fRecos;} - AliRICHChamber* C(Int_t iC) const{return (AliRICHChamber*)fChambers->At(iC-1);} - AliRICHParam* P() const{return fpParam;} - void PrintDigits() const{for(Int_t i=0;iAt(i)->Print();} - void PrintClusters() const{for(Int_t i=0;iAt(i)->Print();} + AliRICHChamber* C(Int_t iC) const{return fpParam->C(iC);} //provides pointer to a given chamber + AliRICHParam* P() const{return fpParam;} //provides pointer to a RICH params + AliRICH* R() {return this;} //provides pointer to RICH main object + TVector Counters() const{return fCounters;} //provides a set of counters + void ControlPlots(); //utility + 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 AddHit(Int_t chamber,Int_t tid,TVector3 iX3,TVector3 oX3,Double_t eloss=0) - {TClonesArray &tmp=*fHits;new(tmp[fNhits++])AliRICHhit(chamber,tid,iX3,oX3,eloss);} - inline void AddSDigit(Int_t c,Int_t x,Int_t y,Double_t q,Int_t pid,Int_t tid); - void AddDigit(int c,int x,int y,int q,int cfm,int *tid) - {TClonesArray &tmp=*((TClonesArray*)fDigitsNew->At(c-1));new(tmp[fNdigitsNew[c-1]++])AliRICHdigit(c,x,y,q,cfm,tid[0],tid[1],tid[2]);} + 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); + void AddDigit(int c,TVector pad,int q,int cfm,int *tid)//Add simulated digit + {TClonesArray &tmp=*((TClonesArray*)fDigitsNew->At(c-1));new(tmp[fNdigitsNew[c-1]++])AliRICHdigit(c,pad,q,cfm,tid[0],tid[1],tid[2]);} + void AddDigit(Int_t c,TVector pad,Int_t q)//for real data digits + {TClonesArray &tmp=*((TClonesArray*)fDigitsNew->At(0));new(tmp[fNdigitsNew[0]++])AliRICHdigit(c,pad,q,0,-1,-1,-1);} void AddCluster(AliRICHcluster &cl) - {Int_t c=cl.C()-1;/*cout<At(c));new(tmp[fNclusters[c]++])AliRICHcluster(cl);} - void AddReco(Int_t tid,Double_t thetaCherenkov,Int_t nPhotons) - {TClonesArray &tmp=*(TClonesArray*)fRecos;new(tmp[fNrecos++])AliRICHreco(tid,thetaCherenkov,nPhotons);} - + {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 {kCSI=6,kGAP=9}; + enum {kAir=1,kCSI=6,kGAP=9,kAl=10,kCH4=5,kSteel=15,kPerpex=16,kSr90=17}; AliRICHParam *fpParam; //main RICH parametrization - TObjArray *fChambers; //list of RICH chambers - //fHits and fDigits belong to AliDetector - TClonesArray *fSdigits; //! List of sdigits - Int_t fNsdigits; //! Current number of sdigits + TClonesArray *fSdigits; //! list of sdigits + Int_t fNsdigits; //! current number of sdigits - TObjArray *fDigitsNew; //! Each chamber holds it's one lists of digits - Int_t fNdigitsNew[kNCH]; //! Array of current numbers of digits + TObjArray *fDigitsNew; //! each chamber holds it's one lists of digits + Int_t fNdigitsNew[kNchambers]; //! array of current numbers of digits - TObjArray *fClusters; //! Each chamber holds it's one lists of clusters - Int_t fNclusters[kNCH]; //! Array of current numbers of raw clusters + TObjArray *fClusters; //! each chamber holds it's one lists of clusters + Int_t fNclusters[kNchambers]; //! array of current numbers of raw clusters - TClonesArray *fRecos; //! pointer to the list of recos - Int_t fNrecos; //! number of recos - - ClassDef(AliRICH,5) //Main RICH class + TVector fCounters; //Photon history conters, explanation in StepManager() + ClassDef(AliRICH,7) //Main RICH class };//class AliRICH //__________________________________________________________________________________________________ @@ -261,34 +246,27 @@ void AliRICH::CreateDigits() { if(fDigitsNew) return; if(GetDebug())Info("CreateDigits","creating digits containers."); - fDigitsNew = new TObjArray(kNCH); - for(Int_t i=0;iAddAt(new TClonesArray("AliRICHdigit",10000), i); fNdigitsNew[i]=0;} + fDigitsNew = new TObjArray(kNchambers); + for(Int_t i=0;iAddAt(new TClonesArray("AliRICHdigit",10000), i); fNdigitsNew[i]=0;} } //__________________________________________________________________________________________________ void AliRICH::CreateClusters() { if(fClusters) return; if(GetDebug())Info("CreateClusters","creating clusters containers."); - fClusters = new TObjArray(kNCH); - for(Int_t i=0;iAddAt(new TClonesArray("AliRICHcluster",10000), i); fNclusters[i]=0;} + fClusters = new TObjArray(kNchambers); + for(Int_t i=0;iAddAt(new TClonesArray("AliRICHcluster",10000), i); fNclusters[i]=0;} } //__________________________________________________________________________________________________ -void AliRICH::CreateRecos() -{ - if(fRecos) return; - if(GetDebug())Info("CreateRecos","creating recos containers."); - fRecos = new TClonesArray("AliRICHreco",1000);fNrecos=0; -} -//__________________________________________________________________________________________________ -void AliRICH::AddSDigit(Int_t c,Int_t x,Int_t y,Double_t q,Int_t pid,Int_t tid) -{ +void AliRICH::AddSDigit(Int_t c,TVector pad,Double_t q,Int_t pid,Int_t tid) +{ + Int_t cfm; switch(pid){ - case 50000050: pid=1000000;break;//cerenkov - case 50000051: pid=1000; break;//feedback - default: pid=1; break;//mip + case 50000050: cfm=1000000;break;//cerenkov + case 50000051: cfm=1000; break;//feedback + default: cfm=1; break;//mip } - TClonesArray &tmp=*fSdigits; - new(tmp[fNsdigits++])AliRICHdigit(c,x,y,q,pid,tid,kBad,kBad); -}//AddSDigit() + TClonesArray &tmp=*fSdigits; new(tmp[fNsdigits++])AliRICHdigit(c,pad,q,cfm,tid,kBad,kBad); +} //__________________________________________________________________________________________________ #endif//#ifndef AliRICH_h diff --git a/RICH/AliRICHChamber.cxx b/RICH/AliRICHChamber.cxx index 5fd0198cd6c..7f7f86f3205 100644 --- a/RICH/AliRICHChamber.cxx +++ b/RICH/AliRICHChamber.cxx @@ -14,67 +14,59 @@ // ************************************************************************** #include "AliRICHChamber.h" -#include "AliRICHParam.h" #include ClassImp(AliRICHChamber) -//______________________________________________________________________________ -AliRICHChamber::AliRICHChamber() -{ -//default ctor - fpParam=0; - fpRotMatrix=0; -} //______________________________________________________________________________ -AliRICHChamber::AliRICHChamber(Int_t iModuleN,AliRICHParam *pParam) +AliRICHChamber::AliRICHChamber(Int_t iChamber):TNamed() { //main ctor. Defines all geometry parameters for the given module. SetToZenith();//put to up position - switch(iModuleN){ + switch(iChamber){ case 1: - RotateX(-pParam->AngleYZ()); - RotateZ(-pParam->AngleXY()); + RotateX(-AliRICHParam::AngleYZ()); + RotateZ(-AliRICHParam::AngleXY()); fName="RICHc1";fTitle="RICH chamber 1"; break; case 2: - RotateZ(-pParam->AngleXY()); + RotateZ(-AliRICHParam::AngleXY()); fName="RICHc2";fTitle="RICH chamber 2"; break; case 3: - RotateX(-pParam->AngleYZ()); + RotateX(-AliRICHParam::AngleYZ()); fName="RICHc3";fTitle="RICH chamber 3"; break; case 4: fName="RICHc4";fTitle="RICH chamber 4"; //no turns break; case 5: - RotateX(pParam->AngleYZ()); + RotateX(AliRICHParam::AngleYZ()); fName="RICHc5";fTitle="RICH chamber 5"; break; case 6: - RotateZ(pParam->AngleXY()); + RotateZ(AliRICHParam::AngleXY()); fName="RICHc6";fTitle="RICH chamber 6"; break; case 7: - RotateX(pParam->AngleYZ()); - RotateZ(pParam->AngleXY()); + RotateX(AliRICHParam::AngleYZ()); + RotateZ(AliRICHParam::AngleXY()); fName="RICHc7";fTitle="RICH chamber 7"; break; default: - Fatal("named ctor","Wrong chamber number %i, check CreateChamber ctor",iModuleN); + Fatal("named ctor","Wrong chamber number %i, check CreateChamber ctor",iChamber); }//switch(iModuleN) - RotateZ(pParam->AngleRot());//apply common rotation + RotateZ(AliRICHParam::AngleRot());//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()); - fpParam=pParam; }//main ctor //__________________________________________________________________________________________________ -void AliRICHChamber::Print(Option_t *) const +void AliRICHChamber::Print(Option_t *opt) const { -//debug printout method - 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()); +// 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); }//Print() //__________________________________________________________________________________________________ diff --git a/RICH/AliRICHChamber.h b/RICH/AliRICHChamber.h index c6ff33359c5..08ebdb9db52 100644 --- a/RICH/AliRICHChamber.h +++ b/RICH/AliRICHChamber.h @@ -9,19 +9,14 @@ #include #include #include "AliRICHParam.h" -#include "AliSegmentation.h" -class AliRICHGeometry; -class AliRICHResponse; class TRotMatrix; -typedef enum {kMip, kPhoton} ResponseType; -class AliRICHParam; class AliRICHChamber : public TNamed { public: - AliRICHChamber(); - AliRICHChamber(Int_t iModuleN,AliRICHParam *pParam); + AliRICHChamber():TNamed(),fpRotMatrix(0) {;} + AliRICHChamber(Int_t iChamberN); AliRICHChamber(const AliRICHChamber &chamber):TNamed(chamber) {;} virtual ~AliRICHChamber() {;} AliRICHChamber& operator=(const AliRICHChamber&) {return *this;} @@ -29,65 +24,45 @@ public: TRotMatrix* RotMatrix() const{return fpRotMatrix;} TString RotMatrixName() const{return "rot"+fName;} TRotation Rot() const{return fRot;} - Double_t Rho() const{return fCenterV3.Mag();} - Double_t ThetaD() const{return fCenterV3.Theta()*TMath::RadToDeg();} - Double_t PhiD() const{return fCenterV3.Phi()*TMath::RadToDeg();} - 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);} - void RotateY(Double_t a) {fRot.RotateY(a);fCenterV3.RotateY(a);fPcX3.RotateY(a);} - void RotateZ(Double_t a) {fRot.RotateZ(a);fCenterV3.RotateZ(a);fPcX3.RotateZ(a);} + 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 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();} - TVector3 L2G(TVector3 x3) const{x3.Transform(fRot);x3+=fCenterV3;return x3;} - TVector3 G2L(TVector3 x3) const{x3-=fCenterV3;x3.Transform(fRot.Inverse()); return x3;} - inline TVector2 Glob2Loc(TVector3 x3, Bool_t isVector=kFALSE) const; - TVector2 Glob2Loc(TLorentzVector x4,Bool_t isVector=kFALSE) const{return Glob2Loc(x4.Vect(),isVector);} - TVector3 L2G(Double_t x,Double_t y,Double_t z) const{return L2G(TVector3(x,y,z));} - TVector3 G2L(TLorentzVector x4) const{return G2L(x4.Vect());} - Float_t G2Ly(TLorentzVector x4) const{TVector3 x3=G2L(x4.Vect()); return x3.Z();} - TVector3 G2L(Double_t x,Double_t y,Double_t z) const{return G2L(TVector3(x,y,z));} - Float_t G2Lx(Double_t x,Double_t y,Double_t z) const{TVector3 x3=G2L(x,y,z); return x3.X();} - Float_t G2Ly(Double_t x,Double_t y,Double_t z) const{TVector3 x3=G2L(x,y,z); return x3.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 - void LocaltoGlobal(Float_t pos[3],Float_t localpos[3]) { - //Transformation from local to global coordinates, chamber-dependant - TVector3 buf = L2G(localpos[0],localpos[1],localpos[2]); - pos[0]=buf.X();pos[1]=buf.Y();pos[2]=buf.Z(); - } - void GlobaltoLocal(Float_t pos[3],Float_t localpos[3]) { - //Transformation from Global to local coordinates, chamber-dependant - TVector3 buf = G2L(pos[0],pos[1],pos[2]); - localpos[0]=buf.X();localpos[1]=buf.Y();localpos[2]=buf.Z(); - } + inline void SetToZenith(); TRotMatrix *GetRotMatrix() const{return fpRotMatrix;} protected: - TVector3 fCenterV3; //chamber center position in MRS (cm) + TVector3 fCenterV3; //chamber 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 - AliRICHParam *fpParam; //main RICH parameters description - ClassDef(AliRICHChamber,5) //single RICH chamber description + ClassDef(AliRICHChamber,6) //single RICH chamber description };//class AliRICHChamber //__________________________________________________________________________________________________ void AliRICHChamber::SetToZenith() { - fCenterV3.SetXYZ(0,AliRICHParam::Offset()-AliRICHParam::GapThickness()/2,0); - fPcX3.SetXYZ(0,AliRICHParam::Offset()-AliRICHParam::GapThickness()/2+5.276+0.25,0); +//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); } //__________________________________________________________________________________________________ -TVector2 AliRICHChamber::Glob2Loc(TVector3 x3,Bool_t isVector)const -{ - if(!isVector) x3-=fPcX3; - x3.Transform(fRot.Inverse()); - return TVector2(x3.X(),x3.Z());//attention Y and Z are misplaced! -} -//__________________________________________________________________________________________________ #endif //AliRICHChamber_h diff --git a/RICH/AliRICHClusterFinder.cxx b/RICH/AliRICHClusterFinder.cxx index 5b4031576cd..c12ff0e9ad0 100644 --- a/RICH/AliRICHClusterFinder.cxx +++ b/RICH/AliRICHClusterFinder.cxx @@ -29,219 +29,209 @@ ClassImp(AliRICHClusterFinder) //__________________________________________________________________________________________________ AliRICHClusterFinder::AliRICHClusterFinder(AliRICH *pRICH) {//main ctor -// Info("main ctor","Start."); + fRICH = pRICH;//this goes first as GetDebug depends on it. + if(GetDebug())Info("main ctor","Start."); - fRICH = pRICH; - - fHitMap = 0; + fDigitMap = 0; fRawCluster.Reset(); fResolvedCluster.Reset(); - + if(GetDebug())Info("main ctor","Stop."); }//main ctor //__________________________________________________________________________________________________ -void AliRICHClusterFinder::FindLocalMaxima() -{ -// Find number of local maxima in the cluster -// Info("FindLocalMaxima","Start."); - - fNlocals = 0; - for(Int_t iDig1=0;iDig1At(iDig1); - Int_t padX1 = pDig1->X(); - Int_t padY1 = pDig1->Y(); - Int_t padQ1 = (Int_t)(pDig1->Q()+0.1); - Int_t padC1 = pDig1->ChFbMi(); - for(Int_t iDig2=0;iDig2At(iDig2); - Int_t padX2 = pDig2->X(); - Int_t padY2 = pDig2->Y(); - Int_t padQ2 = (Int_t)(pDig2->Q()+0.1); - if(iDig1==iDig2) continue; - Int_t diffx = TMath::Sign(padX1-padX2,1); - Int_t diffy = TMath::Sign(padY1-padY2,1); - if((diffx+diffy)<=1) { - if(padQ2>=padQ1) iNotMax++; - } - } - if(iNotMax==0) { - if (fNlocals >= fLocalX.GetSize()) { - fLocalX.Set(fNlocals+100); - fLocalY.Set(fNlocals+100); - fLocalQ.Set(fNlocals+100); - fLocalC.Set(fNlocals+100); - } - TVector2 x2=AliRICHParam::Pad2Loc(padX1,padY1); - fLocalX[fNlocals]=x2.X();fLocalY[fNlocals]=x2.Y(); - fLocalQ[fNlocals] = (Double_t)padQ1; - fLocalC[fNlocals] = padC1; - fNlocals++; - } - } -}//FindLocalMaxima() -//__________________________________________________________________________________________________ void AliRICHClusterFinder::Exec() { +//Main method of cluster finder. Loops on events and chambers, everything else is done in FindClusters() if(GetDebug()) Info("Exec","Start."); - - - Rich()->GetLoader()->LoadDigits(); - - Rich()->GetLoader()->GetRunLoader()->LoadHeader(); - Rich()->GetLoader()->GetRunLoader()->LoadKinematics(); + + R()->GetLoader() ->LoadDigits(); + R()->GetLoader()->GetRunLoader()->LoadHeader(); + R()->GetLoader()->GetRunLoader()->LoadKinematics(); for(Int_t iEventN=0;iEventNGetEventsPerRun();iEventN++){//events loop - if(GetDebug()) Info("Exec","Event %i processed.",iEventN+1); -// gAlice->GetRunLoader()->GetEvent(iEventN); - Rich()->GetLoader()->GetRunLoader()->GetEvent(iEventN); + if(GetDebug()) Info("Exec","Processing event %i...",iEventN); + R()->GetLoader()->GetRunLoader()->GetEvent(iEventN); - Rich()->GetLoader()->MakeTree("R"); Rich()->MakeBranch("R"); - Rich()->ResetDigits(); Rich()->ResetClusters(); + R()->GetLoader()->MakeTree("R"); R()->MakeBranch("R"); + R()->ResetDigits(); R()->ResetClusters(); - Rich()->GetLoader()->TreeD()->GetEntry(0); + R()->GetLoader()->TreeD()->GetEntry(0); for(Int_t iChamber=1;iChamber<=kNCH;iChamber++){//chambers loop FindClusters(iChamber); }//chambers loop - Rich()->GetLoader()->TreeR()->Fill(); - Rich()->GetLoader()->WriteRecPoints("OVERWRITE"); + R()->GetLoader()->TreeR()->Fill(); R()->GetLoader()->WriteRecPoints("OVERWRITE");//write out clusters for current event }//events loop - Rich()->GetLoader()->UnloadDigits(); Rich()->GetLoader()->UnloadRecPoints(); - Rich()->ResetDigits(); Rich()->ResetClusters(); - Rich()->GetLoader()->GetRunLoader()->UnloadHeader(); - Rich()->GetLoader()->GetRunLoader()->UnloadKinematics(); + R()->ResetDigits();//reset and unload everything + R()->ResetClusters(); + R()->GetLoader() ->UnloadDigits(); + R()->GetLoader() ->UnloadRecPoints(); + R()->GetLoader()->GetRunLoader()->UnloadHeader(); + R()->GetLoader()->GetRunLoader()->UnloadKinematics(); if(GetDebug()) Info("Exec","Stop."); }//Exec() //__________________________________________________________________________________________________ void AliRICHClusterFinder::FindClusters(Int_t iChamber) { - //finds neighbours and fill the tree with raw clusters - Int_t nDigits=Rich()->Digits(iChamber)->GetEntriesFast(); -// Info("FindClusters","Start for Chamber %i with %i digits.",iChamber,nDigits); - if(nDigits==0)return; +//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); + + if(iNdigits==0)return;//no digits for a given chamber, nothing to do - fHitMap=new AliRICHMap(Rich()->Digits(iChamber));//create digit map for the given chamber + fDigitMap=new AliRICHMap(R()->Digits(iChamber));//create digit map for the given chamber - for(Int_t iDig=0;iDigDigits(iChamber)->At(iDig); + for(Int_t iDigN=0;iDigNDigits(iChamber)->At(iDigN); Int_t i=dig->X(); Int_t j=dig->Y(); - if(fHitMap->TestHit(i,j)==kUsed) continue; + if(fDigitMap->TestHit(i,j)==kUsed) continue;//this digit is already taken, go after next digit - FormRawCluster(i,j); - FindLocalMaxima(); -// cout << " fNlocals in FindCluster " << fNlocals << endl; - fRawCluster.CoG(fNlocals); // first initial approxmation of the CoG...to start minimization. + FormRawCluster(i,j);//form raw cluster starting from (i,j) pad + if(GetDebug()) {Info("FindClusters","After FormRawCluster:");fRawCluster.Print();} + FindLocalMaxima(); //find number of local maxima and initial center of gravity + if(GetDebug()) {Info("FindClusters","After FindLocalMaxima:");fRawCluster.Print();} - if(AliRICHParam::IsResolveClusters()) { - ResolveCluster(); // ResolveCluster serialization will happen inside - } else { - WriteRawCluster(); // simply output of the RawCluster found without deconvolution + if(AliRICHParam::IsResolveClusters()&&fRawCluster.Size()>1&&fRawCluster.Size()<6){ + FitCoG(); //serialization of resolved clusters will happen inside + }else{//cluster size=1 or resolving is switched off + WriteRawCluster();//simply output the formed raw cluster without deconvolution } - fRawCluster.Reset(); - fResolvedCluster.Reset(); - }//digits loop + fRawCluster.Reset(); fResolvedCluster.Reset(); + }//digits loop for a given chamber - delete fHitMap; -// Info("FindClusters","Stop."); + delete fDigitMap; + if(GetDebug())Info("FindClusters","Stop."); }//FindClusters() //__________________________________________________________________________________________________ void AliRICHClusterFinder::FindClusterContribs(AliRICHcluster *pCluster) { -// finds CombiPid for a given cluster -// Info("FindClusterContribs","Start"); +//Finds cerenkov-feedback-mip mixture for a given cluster + if(GetDebug()) {Info("FindClusterContribs","Start.");pCluster->Print();} TObjArray *pDigits = pCluster->Digits(); + if(!pDigits) return; //?????????? Int_t iNmips=0,iNckovs=0,iNfeeds=0; TArrayI contribs(3*pCluster->Size()); Int_t *pindex = new Int_t[3*pCluster->Size()]; for(Int_t iDigN=0;iDigNSize();iDigN++) {//loop on digits of a given cluster - contribs[3*iDigN] =((AliRICHdigit*)pDigits->At(iDigN))->Tid(0); + contribs[3*iDigN] =((AliRICHdigit*)pDigits->At(iDigN))->GetTrack(0); if (contribs[3*iDigN] >= 10000000) contribs[3*iDigN] = 0; - contribs[3*iDigN+1]=((AliRICHdigit*)pDigits->At(iDigN))->Tid(1); + contribs[3*iDigN+1]=((AliRICHdigit*)pDigits->At(iDigN))->GetTrack(1); if (contribs[3*iDigN+1] >= 10000000) contribs[3*iDigN+1] = 0; - contribs[3*iDigN+2]=((AliRICHdigit*)pDigits->At(iDigN))->Tid(2); + contribs[3*iDigN+2]=((AliRICHdigit*)pDigits->At(iDigN))->GetTrack(2); if (contribs[3*iDigN+2] >= 10000000) contribs[3*iDigN+2] = 0; }//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 Tid + 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); if(contribs[pindex[iDigN]]!=contribs[pindex[iDigN+1]]) { - TParticle* particle = Rich()->GetLoader()->GetRunLoader()->Stack()->Particle(contribs[pindex[iDigN]]); + TParticle* particle = R()->GetLoader()->GetRunLoader()->Stack()->Particle(contribs[pindex[iDigN]]); Int_t code = particle->GetPdgCode(); Double_t charge = 0; - if (particle->GetPDG()) particle->GetPDG()->Charge(); + if(particle->GetPDG()) charge=particle->GetPDG()->Charge(); + if(GetDebug()) Info("FindClusterContribs"," charge of particle %f",charge); if(code==50000050) iNckovs++; - else if(code==50000051) iNfeeds++; - else if(charge!=0) iNmips++; - if (contribs[pindex[iDigN+1]]==kBad) break; + if(code==50000051) iNfeeds++; + if(charge!=0) iNmips++; } }//loop on digits to sort Tid - pCluster->SetCombiPid(iNckovs,iNfeeds,iNmips); -// pCluster->Print(); + + if (contribs[pindex[3*pCluster->Size()-1]]!=kBad) { + + TParticle* particle = R()->GetLoader()->GetRunLoader()->Stack()->Particle(contribs[pindex[3*pCluster->Size()-1]]); + 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); + if(code==50000050) iNckovs++; + if(code==50000051) iNfeeds++; + if(charge!=0) iNmips++; + } + + pCluster->CFM(iNckovs,iNfeeds,iNmips); +// delete [] pindex; + if(GetDebug()) {pCluster->Print();Info("FindClusterContribs","Stop.");} }//FindClusterContribs() //__________________________________________________________________________________________________ void AliRICHClusterFinder::FormRawCluster(Int_t i, Int_t j) { -// Builder of the final Raw Cluster (before deconvolution) - if(GetDebug()) Info("FormRawCluster","Start with digit(%i,%i)",i,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()); - fRawCluster.AddDigit((AliRICHdigit*) fHitMap->GetHit(i,j)); - fHitMap->FlagHit(i,j);// Flag hit as taken + fRawCluster.AddDigit((AliRICHdigit*) fDigitMap->GetHit(i,j));//take this pad in cluster + fDigitMap->FlagHit(i,j);//flag this pad as taken Int_t listX[4], listY[4]; // Now look recursively for all neighbours - for (Int_t iNeighbour=0;iNeighbourP()->PadNeighbours(i,j,listX,listY);iNeighbour++) - if(fHitMap->TestHit(listX[iNeighbour],listY[iNeighbour])==kUnused) - FormRawCluster(listX[iNeighbour],listY[iNeighbour]); + for (Int_t iNei=0;iNeiP()->PadNeighbours(i,j,listX,listY);iNei++) + if(fDigitMap->TestHit(listX[iNei],listY[iNei])==kUnused) FormRawCluster(listX[iNei],listY[iNei]); }//FormRawCluster() //__________________________________________________________________________________________________ -void AliRICHClusterFinder::ResolveCluster() -{// Decluster algorithm - if(GetDebug()) {Info("ResolveCluster","Start."); fRawCluster.Print();} - - switch (fRawCluster.Size()) { - - case 1: // nothing to decluster: cluster size = 1 - WriteRawCluster();break; - default: // cluster size > 1: if=2 FitCoG; if>2 Resolve and FitCoG - FitCoG();break; - } -}//ResolveCluster() +void AliRICHClusterFinder::FindLocalMaxima() +{ +//find number of local maxima in the current raw cluster and then calculates initial center of gravity + fNlocals=0; + for(Int_t iDig1=0;iDig1At(iDig1); + TVector pad1 = pDig1->Pad(); + Int_t padQ1 = (Int_t)(pDig1->Q()+0.1); + Int_t padC1 = pDig1->ChFbMi(); + for(Int_t iDig2=0;iDig2At(iDig2); + TVector pad2 = pDig2->Pad(); + Int_t padQ2 = (Int_t)(pDig2->Q()+0.1); + if(iDig1==iDig2) continue; + Int_t diffx = TMath::Sign(Int_t(pad1[0]-pad2[0]),1); + Int_t diffy = TMath::Sign(Int_t(pad1[1]-pad2[1]),1); + if((diffx+diffy)<=1) { + if(padQ2>=padQ1) iNotMax++; + } + } + if(iNotMax==0) { + TVector2 x2=AliRICHParam::Pad2Loc(pad1); + fLocalX[fNlocals]=x2.X();fLocalY[fNlocals]=x2.Y(); + fLocalQ[fNlocals] = (Double_t)padQ1; + fLocalC[fNlocals] = padC1; + fNlocals++; + } + } + fRawCluster.CoG(fNlocals); //first initial approximation of the CoG...to start minimization. +}//FindLocalMaxima() //__________________________________________________________________________________________________ void AliRICHClusterFinder::WriteRawCluster() { -// out the current raw cluster +//Add the current raw cluster to the list of clusters if(GetDebug()) Info("WriteRawCluster","Start."); - FindClusterContribs(&fRawCluster); - if(GetDebug()) fRawCluster.Dump(); - Rich()->AddCluster(fRawCluster); -// fRawCluster.Print(); + FindClusterContribs(&fRawCluster); + R()->AddCluster(fRawCluster); + + if(GetDebug()){fRawCluster.Print(); Info("WriteRawCluster","Stop.");} }//WriteRawCluster() //__________________________________________________________________________________________________ void AliRICHClusterFinder::WriteResolvedCluster() { -// out the current resolved cluster +//Add the current resolved cluster to the list of clusters if(GetDebug()) Info("WriteResolvedCluster","Start."); -// FindClusterContribs(&fResolvedCluster); - Rich()->AddCluster(fResolvedCluster); + FindClusterContribs(&fResolvedCluster); + R()->AddCluster(fResolvedCluster); + if(GetDebug()){fResolvedCluster.Print(); Info("WriteResolvedCluster","Stop.");} }//WriteResolvedCluster() //__________________________________________________________________________________________________ void AliRICHClusterFinder::FitCoG() -{// Fit cluster size 2 by single Mathieson - if(GetDebug()) Info("FitCoG","Start with size %3i and local maxima %3i",fRawCluster.Size(),fNlocals); +{ +//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();} Double_t arglist; Int_t ierflag = 0; -// fRawCluster.Print(); -// if(fNlocals==0) fRawCluster.Print(); - if(fNlocals==0||fNlocals>3) {WriteRawCluster();return;} - TMinuit *pMinuit = new TMinuit(3*fNlocals-1); pMinuit->mninit(5,10,7); @@ -265,19 +255,15 @@ void AliRICHClusterFinder::FitCoG() lower = vstart - 2*AliRICHParam::PadSizeX(); upper = vstart + 2*AliRICHParam::PadSizeX(); pMinuit->mnparm(3*i ,Form("xCoG %i",i),vstart,stepX,lower,upper,ierflag); -// cout << Form("xCoG %i",i) << "start " << vstart << "lower " << lower << "upper " << upper << endl; vstart = fLocalY[i]; lower = vstart - 2*AliRICHParam::PadSizeY(); upper = vstart + 2*AliRICHParam::PadSizeY(); pMinuit->mnparm(3*i+1,Form("yCoG %i",i),vstart,stepY,lower,upper,ierflag); -// cout << Form("yCoG %i",i) << "start " << vstart << "lower " << lower << "upper " << upper << endl; if(i==fNlocals-1) break; // last parameter is constrained vstart = fLocalQ[i]/fRawCluster.Q(); lower = 0; upper = 1; pMinuit->mnparm(3*i+2,Form("qfrac %i",i),vstart,stepQ,lower,upper,ierflag); -// cout << Form("qCoG %i",i) << "start " << vstart << "lower " << lower << "upper " << upper << endl; - } arglist = -1; @@ -303,23 +289,16 @@ void AliRICHClusterFinder::FitCoG() qfracCoG[fNlocals-1] = 1 - qfraclast; delete pMinuit; - for(Int_t i=0;iGetObjectFit())->GetRawCluster(); @@ -338,19 +317,13 @@ void RICHMinMathieson(Int_t &npar, Double_t *, Double_t &chi2, Double_t *par, In chi2 = 0; Int_t qtot = pRawCluster->Q(); for(Int_t i=0;iSize();i++) { - Int_t padX = ((AliRICHdigit *)pRawCluster->Digits()->At(i))->X(); - Int_t padY = ((AliRICHdigit *)pRawCluster->Digits()->At(i))->Y(); + TVector pad=((AliRICHdigit *)pRawCluster->Digits()->At(i))->Pad(); Double_t padQ = ((AliRICHdigit *)pRawCluster->Digits()->At(i))->Q(); Double_t qfracpar=0; for(Int_t j=0;j 2 - void FitCoG(); //Evaluate the CoG as the best - void WriteRawCluster(); //write in the list of cluster - void WriteResolvedCluster(); //write in the list of cluster - AliRICHcluster *GetRawCluster() {return &fRawCluster;} //Return pointer to the current raw cluster - Bool_t GetDebug() const {return fRICH->GetDebug();} //is debug printout needed? + AliRICH *R() {return fRICH;} //returns pointer to RICH + void Exec(); //loop on events and chambers + void FindClusters(Int_t iChamber); //find all clusters for a given chamber + void FindClusterContribs(AliRICHcluster *pCluster); //find CFM for the current cluster + void FormRawCluster(Int_t i, Int_t j); //form a raw cluster + void FindLocalMaxima(); //find local maxima in a cluster + void FitCoG(); //evaluate the CoG as the best + void WriteRawCluster(); //write in the list of cluster + void WriteResolvedCluster(); //write in the list of cluster + AliRICHcluster *GetRawCluster() {return &fRawCluster;} //returns pointer to the current raw cluster + Bool_t GetDebug() const{return fRICH->GetDebug();} //is debug printout needed? protected: - AliRICH *fRICH; //Pointer to RICH - AliHitMap *fHitMap; //Hit Map with digit positions - AliRICHcluster fRawCluster; //Current raw cluster before deconvolution - AliRICHcluster fResolvedCluster; //Current cluster after deconvolution - Int_t fNlocals; // number of local maxima - TArrayD fLocalX,fLocalY; // list of locals X,Y - TArrayD fLocalQ; // list of locals charge Q - TArrayI fLocalC; // list of locals CombiPid - ClassDef(AliRICHClusterFinder,0) //Finds raw clusters, trasfers them to resolved clusters through declustering. + AliRICH *fRICH; //pointer to RICH + AliHitMap *fDigitMap; //map of digits positions + AliRICHcluster fRawCluster; //current raw cluster before deconvolution + AliRICHcluster fResolvedCluster; //current cluster after deconvolution + Int_t fNlocals; //number of local maxima + Double_t fLocalX[100],fLocalY[100]; //list of locals X,Y + Double_t fLocalQ[100]; //list of locals charge Q + Int_t fLocalC[100]; //list of locals CombiPid + ClassDef(AliRICHClusterFinder,0) //finds raw clusters, trasfers them to resolved clusters through declustering. }; #endif diff --git a/RICH/AliRICHDigitizer.cxx b/RICH/AliRICHDigitizer.cxx index 3eaa2806e80..dbca6d8572b 100644 --- a/RICH/AliRICHDigitizer.cxx +++ b/RICH/AliRICHDigitizer.cxx @@ -36,16 +36,15 @@ Bool_t AliRICHDigitizer::Init() { //This methode is called from AliRunDigitizer after the corresponding file is open if(GetDebug())Info("Init","Start."); - AliRunLoader *pOutAL = - AliRunLoader::GetRunLoader(fManager->GetOutputFolderName()); + AliRunLoader *pOutAL = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName()); if (!pOutAL->GetAliRun()) pOutAL->LoadgAlice(); - fRich=(AliRICH*)pOutAL->GetAliRun()->GetDetector("RICH"); - Rich()->P()->GenSigmaThMap(); + fRich=(AliRICH*)pOutAL->GetAliRun()->GetDetector("RICH");//retrive RICH pointer from OUTPUT stream return kTRUE; }//Init() //__________________________________________________________________________________________________ void AliRICHDigitizer::Exec(Option_t*) { +//this method invoked if(GetDebug())Info("Exec","Start with %i input(s) for event %i",fManager->GetNinputs(),fManager->GetOutputEventNr()); AliRunLoader *pInAL=0, *pOutAL;//in and out Run loaders @@ -53,47 +52,52 @@ void AliRICHDigitizer::Exec(Option_t*) pOutAL = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName()); pOutRL = pOutAL->GetLoader("RICHLoader"); - pOutRL->MakeTree("D"); Rich()->MakeBranch("D"); //create TreeD with RICH branches in output stream + pOutRL->MakeTree("D"); R()->MakeBranch("D"); //create TreeD with RICH branches in output stream TClonesArray tmpCA("AliRICHdigit");//tmp storage for sdigits sum up from all input files Int_t total=0; for(Int_t inFileN=0;inFileNGetNinputs();inFileN++){//files loop - pInAL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inFileN)); + pInAL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inFileN));//get run loader from current input pInRL = pInAL->GetLoader("RICHLoader"); if(pInRL==0) continue;//no RICH in this input, check the next input if (!pInAL->GetAliRun()) pInAL->LoadgAlice(); AliRICH* rich=(AliRICH*)pInAL->GetAliRun()->GetDetector("RICH"); pInRL->LoadSDigits(); pInRL->TreeS()->GetEntry(0); - Info("Exec","input %i has %i sdigits",inFileN,rich->SDigits()->GetEntries()); - for(Int_t i=0;iSDigits()->GetEntries();i++) { + if(GetDebug())Info("Exec","input %i has %i sdigits",inFileN,rich->SDigits()->GetEntries()); + for(Int_t i=0;iSDigits()->GetEntries();i++){//collect sdigits from current input to tmpCA new(tmpCA[total++]) AliRICHdigit(*(AliRICHdigit*)rich->SDigits()->At(i)); ((AliRICHdigit*)tmpCA[total-1])->AddTidOffset(fManager->GetMask(inFileN)); } pInRL->UnloadSDigits(); rich->ResetSDigits(); }//files loop - tmpCA.Sort(); //sort them according to Id() methode + tmpCA.Sort(); //sort them according to Id() method - Int_t chFbMip=0,chamber=0,x=0,y=0,tid[3],id=0; Double_t q=0; + if(GetDebug()) {tmpCA.Print();Info("Exec","Totally %i sdigits in %i inputs",tmpCA.GetEntries(),fManager->GetNinputs());} + Int_t cfm=0,chamber=0,id=0; //cfm is cerenkov feedback mip mixture + TVector pad(2); pad[0]=0;pad[1]=0; + Double_t q=0; + Int_t tid[3]={0,0,0}; Int_t iNdigitsPerPad=0;//how many sdigits for a given pad for(Int_t i=0;iId()==id){//still the same pad - iNdigitsPerPad++; q+=pSdig->Q(); chFbMip+=pSdig->ChFbMi();//sum up charge and cfm - if(iNdigitsPerPad<=3) tid[iNdigitsPerPad-1]=pSdig->Tid(0); + iNdigitsPerPad++; q+=pSdig->Q(); cfm+=pSdig->ChFbMi();//sum up charge and cfm + if(pSdig->ChFbMi()==1) tid[0] = pSdig->GetTrack(0); // force the first tid to be mip's tid if it exists in the current pad + if(iNdigitsPerPad<=3) tid[iNdigitsPerPad-1]=pSdig->GetTrack(0); else if(GetDebug())Warning("Exec","More then 3 sdigits for the given pad"); }else{//new pad, add the pevious one - if(id!=kBad&&Rich()->P()->IsOverTh(chamber,x,y,q)) Rich()->AddDigit(chamber,x,y,(Int_t)q,chFbMip,tid); //add newly created dig - chFbMip=pSdig->ChFbMi(); chamber=pSdig->C(); id=pSdig->Id(); x=pSdig->X(); y=pSdig->Y(); q=pSdig->Q(); //init all values by current sdig - iNdigitsPerPad=1; tid[0]=pSdig->Tid(0); tid[1]=tid[2]=kBad; + if(id!=kBad&&R()->P()->IsOverTh(chamber,pad,q)) R()->AddDigit(chamber,pad,(Int_t)q,cfm,tid); //add newly created dig + cfm=pSdig->ChFbMi(); chamber=pSdig->C(); id=pSdig->Id(); pad=pSdig->Pad(); q=pSdig->Q(); //init all values by current sdig + iNdigitsPerPad=1; tid[0]=pSdig->GetTrack(0); tid[1]=tid[2]=kBad; } }//sdigits loop (sorted) - if(tmpCA.GetEntries()&&Rich()->P()->IsOverTh(chamber,x,y,q)) Rich()->AddDigit(chamber,x,y,(Int_t)q,chFbMip,tid);//add the last dig + if(tmpCA.GetEntries()&&R()->P()->IsOverTh(chamber,pad,q)) R()->AddDigit(chamber,pad,(Int_t)q,cfm,tid);//add the last dig pOutRL->TreeD()->Fill(); //fill the tree with the list of digits pOutRL->WriteDigits("OVERWRITE"); //serialize them to file tmpCA.Clear(); - pOutRL->UnloadDigits(); Rich()->ResetDigits(); + pOutRL->UnloadDigits(); R()->ResetDigits(); if(GetDebug())Info("Exec","Stop."); }//Exec() diff --git a/RICH/AliRICHDigitizer.h b/RICH/AliRICHDigitizer.h index c907b652b48..d6b24ffdc9b 100644 --- a/RICH/AliRICHDigitizer.h +++ b/RICH/AliRICHDigitizer.h @@ -13,7 +13,7 @@ class AliRICH; class AliRICHDigitizer : public AliDigitizer { public: - AliRICHDigitizer() {;} + AliRICHDigitizer():fRich(0) {;} AliRICHDigitizer(AliRunDigitizer * manager):AliDigitizer(manager) {if(GetDebug())Info("main ctor","Start.");} virtual ~AliRICHDigitizer() {if(GetDebug())Info("dtor","Start.");} @@ -21,7 +21,7 @@ public: void Exec(Option_t* option=0); //virtual Bool_t Init(); //virtual Bool_t GetDebug() const {return ((gAlice) ? gAlice->GetDebug() : kFALSE);} - AliRICH* Rich() const {return fRich;} + AliRICH* R() const {return fRich;} //returns pointer to RICH protected: AliRICH* fRich; //pointer to main RICH object ClassDef(AliRICHDigitizer,0) diff --git a/RICH/AliRICHDisplFast.cxx b/RICH/AliRICHDisplFast.cxx index 610e55f3311..656872f5339 100644 --- a/RICH/AliRICHDisplFast.cxx +++ b/RICH/AliRICHDisplFast.cxx @@ -28,7 +28,7 @@ ClassImp(AliRICHDisplFast) //__________________________________________________________________________________________________ -void AliRICHDisplFast::Exec() +void AliRICHDisplFast::Exec(Option_t *) { AliRICH *pRich = (AliRICH*)gAlice->GetDetector("RICH"); Bool_t isHits =!pRich->GetLoader()->LoadHits(); @@ -42,14 +42,14 @@ void AliRICHDisplFast::Exec() TH2F *pHitsH2=0,*pDigitsH2=0,*pClustersH2=0; - if(isHits) pHitsH2 = new TH2F("pHitsH2" , "Event Display",165,-AliRICHParam::PcSizeX()/2,AliRICHParam::PcSizeX()/2, - 144,-AliRICHParam::PcSizeY()/2,AliRICHParam::PcSizeY()/2); + if(isHits) pHitsH2 = new TH2F("pHitsH2" , "Event Display;x,cm;y,cm",165,0,AliRICHParam::PcSizeX(), + 144,0,AliRICHParam::PcSizeY()); if(pHitsH2) pHitsH2->SetStats(kFALSE); - if(isDigits) pDigitsH2 = new TH2F("pDigitsH2" ,"Event Display",165,-AliRICHParam::PcSizeX()/2,AliRICHParam::PcSizeX()/2, - 144,-AliRICHParam::PcSizeY()/2,AliRICHParam::PcSizeY()/2); - if(isClusters) pClustersH2 = new TH2F("pClustersH2","Event Display",165,-AliRICHParam::PcSizeX()/2,AliRICHParam::PcSizeX()/2, - 144,-AliRICHParam::PcSizeY()/2,AliRICHParam::PcSizeY()/2); + if(isDigits) pDigitsH2 = new TH2F("pDigitsH2" ,"Event Display",165,0,AliRICHParam::PcSizeX(), + 144,0,AliRICHParam::PcSizeY()); + if(isClusters) pClustersH2 = new TH2F("pClustersH2","Event Display",165,0,AliRICHParam::PcSizeX(), + 144,0,AliRICHParam::PcSizeY()); @@ -101,7 +101,7 @@ void AliRICHDisplFast::Exec() pRich->GetLoader()->TreeD()->GetEntry(0); for(Int_t j=0;jDigits(iChamber)->GetEntries();j++){//digits loop AliRICHdigit *pDig = (AliRICHdigit*)pRich->Digits(iChamber)->At(j); - TVector2 x2=AliRICHParam::Pad2Loc(pDig->X(),pDig->Y()); + TVector2 x2=AliRICHParam::Pad2Loc(pDig->Pad()); pDigitsH2->Fill(x2.X(),x2.Y(),100); }//digits loop pDigitsH2->SetMarkerColor(kGreen); pDigitsH2->SetMarkerStyle(29); pDigitsH2->SetMarkerSize(0.4); @@ -133,70 +133,19 @@ void AliRICHDisplFast::Exec() //__________________________________________________________________________________________________ void AliRICHDisplFast::DrawSectors() { - Double_t x1[5] = {-AliRICHParam::PcSizeX()/2, - -AliRICHParam::SectorSizeX()/2-AliRICHParam::DeadZone(), - -AliRICHParam::SectorSizeX()/2-AliRICHParam::DeadZone(), - -AliRICHParam::PcSizeX()/2, - -AliRICHParam::PcSizeX()/2}; - Double_t y1[5] = {AliRICHParam::DeadZone()/2, - AliRICHParam::DeadZone()/2, - AliRICHParam::PcSizeY()/2, - AliRICHParam::PcSizeY()/2, - AliRICHParam::DeadZone()/2}; - Double_t x2[5] = {-AliRICHParam::SectorSizeX()/2, - AliRICHParam::SectorSizeX()/2, - AliRICHParam::SectorSizeX()/2, - -AliRICHParam::SectorSizeX()/2, - -AliRICHParam::SectorSizeX()/2}; - Double_t y2[5] = {AliRICHParam::DeadZone()/2, - AliRICHParam::DeadZone()/2, - AliRICHParam::PcSizeY()/2, - AliRICHParam::PcSizeY()/2, - AliRICHParam::DeadZone()/2}; - Double_t x3[5] = { AliRICHParam::SectorSizeX()/2+AliRICHParam::DeadZone(), - AliRICHParam::PcSizeX()/2, - AliRICHParam::PcSizeX()/2, - AliRICHParam::SectorSizeX()/2+AliRICHParam::DeadZone(), - AliRICHParam::SectorSizeX()/2+AliRICHParam::DeadZone()}; - Double_t y3[5] = {AliRICHParam::DeadZone()/2, - AliRICHParam::DeadZone()/2, - AliRICHParam::PcSizeY()/2, - AliRICHParam::PcSizeY()/2, - AliRICHParam::DeadZone()/2}; - Double_t x4[5] = {-AliRICHParam::PcSizeX()/2, - -AliRICHParam::SectorSizeX()/2-AliRICHParam::DeadZone(), - -AliRICHParam::SectorSizeX()/2-AliRICHParam::DeadZone(), - -AliRICHParam::PcSizeX()/2, - -AliRICHParam::PcSizeX()/2}; - Double_t y4[5] = {-AliRICHParam::PcSizeY()/2, - -AliRICHParam::PcSizeY()/2, - -AliRICHParam::DeadZone()/2, - -AliRICHParam::DeadZone()/2, - -AliRICHParam::PcSizeY()/2}; - Double_t x5[5] = {-AliRICHParam::SectorSizeX()/2, - AliRICHParam::SectorSizeX()/2, - AliRICHParam::SectorSizeX()/2, - -AliRICHParam::SectorSizeX()/2, - -AliRICHParam::SectorSizeX()/2}; - Double_t y5[5] = {-AliRICHParam::PcSizeY()/2, - -AliRICHParam::PcSizeY()/2, - -AliRICHParam::DeadZone()/2, - -AliRICHParam::DeadZone()/2, - -AliRICHParam::PcSizeY()/2}; - Double_t x6[5] = { AliRICHParam::SectorSizeX()/2+AliRICHParam::DeadZone(), - AliRICHParam::PcSizeX()/2, - AliRICHParam::PcSizeX()/2, - AliRICHParam::SectorSizeX()/2+AliRICHParam::DeadZone(), - AliRICHParam::SectorSizeX()/2+AliRICHParam::DeadZone()}; - Double_t y6[5] = {-AliRICHParam::PcSizeY()/2, - -AliRICHParam::PcSizeY()/2, - -AliRICHParam::DeadZone()/2, - -AliRICHParam::DeadZone()/2, - -AliRICHParam::PcSizeY()/2}; - TPolyLine *sector1 = new TPolyLine(5,x1,y1); sector1->SetLineColor(21); sector1->Draw(); - TPolyLine *sector2 = new TPolyLine(5,x2,y2); sector2->SetLineColor(21); sector2->Draw(); - TPolyLine *sector3 = new TPolyLine(5,x3,y3); sector3->SetLineColor(21); sector3->Draw(); - TPolyLine *sector4 = new TPolyLine(5,x4,y4); sector4->SetLineColor(21); sector4->Draw(); - TPolyLine *sector5 = new TPolyLine(5,x5,y5); sector5->SetLineColor(21); sector5->Draw(); - TPolyLine *sector6 = new TPolyLine(5,x6,y6); sector6->SetLineColor(21); sector6->Draw(); + AliRICHParam p; + Double_t xLeft[5] = {0,0,p.SectorSizeX(),p.SectorSizeX(),0}; + Double_t xRight[5] = {p.SectorSizeX()+p.DeadZone(),p.SectorSizeX()+p.DeadZone(),p.PcSizeX(),p.PcSizeX(),p.SectorSizeX()+p.DeadZone()}; + + Double_t yDown[5] = {0,p.SectorSizeY(),p.SectorSizeY(),0,0}; + Double_t yCenter[5] = { p.SectorSizeY()+p.DeadZone(),2*p.SectorSizeY()+p.DeadZone(),2*p.SectorSizeY()+p.DeadZone(), + p.SectorSizeY()+p.DeadZone(),p.SectorSizeY()+p.DeadZone()}; + Double_t yUp[5] = {2*p.SectorSizeY()+2*p.DeadZone(),p.PcSizeY(),p.PcSizeY(),2*p.SectorSizeY()+2*p.DeadZone(),2*p.SectorSizeY()+2*p.DeadZone()}; + + TPolyLine *sec1 = new TPolyLine(5,xLeft ,yDown); sec1->SetLineColor(21); sec1->Draw(); + TPolyLine *sec2 = new TPolyLine(5,xRight,yDown); sec2->SetLineColor(21); sec2->Draw(); + TPolyLine *sec3 = new TPolyLine(5,xLeft ,yCenter); sec3->SetLineColor(21); sec3->Draw(); + TPolyLine *sec4 = new TPolyLine(5,xRight,yCenter); sec4->SetLineColor(21); sec4->Draw(); + TPolyLine *sec5 = new TPolyLine(5,xLeft, yUp); sec5->SetLineColor(21); sec5->Draw(); + TPolyLine *sec6 = new TPolyLine(5,xRight,yUp); sec6->SetLineColor(21); sec6->Draw(); }//DrawSectors() diff --git a/RICH/AliRICHDisplFast.h b/RICH/AliRICHDisplFast.h index 9d4ef00876a..f46ccd6afcb 100644 --- a/RICH/AliRICHDisplFast.h +++ b/RICH/AliRICHDisplFast.h @@ -13,10 +13,10 @@ class AliRICH; class AliRICHDisplFast : public TTask { public : - AliRICHDisplFast() {;} - virtual ~AliRICHDisplFast() {;} - static void DrawSectors(); //Draw sectors in plot - void Exec(); //virtual do the main job + AliRICHDisplFast() {;} + virtual ~AliRICHDisplFast() {;} + static void DrawSectors(); //Draw sectors in plot + virtual void Exec(Option_t *opt=0); //virtual do the main job protected: ClassDef(AliRICHDisplFast,0) //Utility class to draw the current event topology }; diff --git a/RICH/AliRICHParam.cxx b/RICH/AliRICHParam.cxx index f8d0a245cf6..ad1679fb9c8 100644 --- a/RICH/AliRICHParam.cxx +++ b/RICH/AliRICHParam.cxx @@ -13,33 +13,30 @@ // * provided "as is" without express or implied warranty. * // ************************************************************************** #include "AliRICHParam.h" -#include +#include "AliRICHChamber.h" ClassImp(AliRICHParam) Bool_t AliRICHParam::fgIsWireSag =kTRUE; Bool_t AliRICHParam::fgIsResolveClusters =kTRUE; +Bool_t AliRICHParam::fgIsRadioSrc =kFALSE; Double_t AliRICHParam::fgAngleRot =-60; -Int_t AliRICHParam::fgHV[kNsectors] ={2150,2100,2050,2000,2150,2150}; +Int_t AliRICHParam::fgHV[kNsectors] ={2050,2050,2050,2050,2050,2050}; Int_t AliRICHParam::fgNsigmaTh =4; -Float_t AliRICHParam::fgSigmaThMean =1.5; -Float_t AliRICHParam::fgSigmaThSpread =0.5; -Float_t AliRICHParam::fSigmaThMap[kNCH][kNpadsX][kNpadsY]; +Float_t AliRICHParam::fgSigmaThMean =1.132; //QDC +Float_t AliRICHParam::fgSigmaThSpread =0.035; // -void AliRICHParam::GenSigmaThMap() -{ -// Generate the map of thresholds sigmas for all pads of all chambers - for(Int_t iChamber=0;iChamberRndm())*SigmaThSpread(); - // Info("GenSigmaThMap"," Threshold map generated for all RICH chambers"); -} //__________________________________________________________________________________________________ -void AliRICHParam::Print() +void AliRICHParam::Print(Option_t*) { - cout<<"\nPads in chamber ("<Print(); } //__________________________________________________________________________________________________ +void AliRICHParam::CreateChambers() +{ +//Create all RICH Chambers on each call. Previous chambers deleted. + if(fpChambers) delete fpChambers; + fpChambers=new TObjArray(kNchambers); + fpChambers->SetOwner(); + for(int iChamberN=0;iChamberNAddAt(new AliRICHChamber(iChamberN+1),iChamberN); +}//void AliRICH::CreateChambers() diff --git a/RICH/AliRICHParam.h b/RICH/AliRICHParam.h index 9d0ac822eac..958e90c84a6 100644 --- a/RICH/AliRICHParam.h +++ b/RICH/AliRICHParam.h @@ -4,38 +4,43 @@ #include #include #include +#include #include #include +#include -static const int kNCH=7; //number of RICH chambers -static const int kNpadsX = 144; //number of pads along X in single chamber -static const int kNpadsY = 160; //number of pads along Y in single chamber +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; // nb. of sectors per chamber +static const int kNsectors=6; //number of sectors per chamber -static const int kadc_satm = 4096; //dynamic range (10 bits) +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 +class AliRICHChamber; class AliRICHParam :public TObject { public: - AliRICHParam() {;} - virtual ~AliRICHParam() {;} - static const Int_t NpadsX() {return kNpadsX;} //pads along X in chamber - static const Int_t NpadsY() {return kNpadsY;} //pads along Y in chamber - static Int_t NpadsXsec() {return NpadsX()/3;} //pads along X in sector - static Int_t NpadsYsec() {return NpadsY()/2;} //pads along Y in sector + 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 + 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 + static Int_t NpadsYsec() {return NpadsY()/3;} //pads along Y in sector static Double_t DeadZone() {return 2.6;} //dead zone size in cm - static Double_t PadSizeX() {return 0.84;} //pad size x in cm - static Double_t PadSizeY() {return 0.8;} //pad size y in cm - static Double_t SectorSizeX() {return NpadsX()*PadSizeX()/3;} //sector size x in cm - static Double_t SectorSizeY() {return NpadsY()*PadSizeY()/2;} //sector size y in cm - static Double_t PcSizeX() {return NpadsX()*PadSizeX()+2*DeadZone();} //photocathode size x in cm - static Double_t PcSizeY() {return NpadsY()*PadSizeY()+DeadZone();} //photocathode size y in cm - static Double_t WirePitch() {return PadSizeX()/2;} //distance between anode wires + 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;} @@ -45,23 +50,32 @@ public: 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 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 GapThickness() {return 8.0;} static Double_t RadiatorToPads() {return FreonThickness()+QuartzThickness()+GapThickness();} - static Double_t ProximityGap() {return 0.445;} - static Double_t AnodeCathodeGap() {return 0.2;} + 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;} - static TVector2 MathiesonDelta() {return TVector2(5*0.18,5*0.18);} - static Int_t MaxQdc() {return 4095;} - static Double_t AlphaFeedback(Int_t sec) {HV(sec);return 0.036;} - + 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? + static Bool_t IsRadioSrc() {return fgIsRadioSrc;} //add radioactive source inside CH4? static Int_t HV(Int_t sector) { if (sector>=1 && sector <=6) return fgHV[sector-1]; @@ -70,188 +84,193 @@ public: return kBad; } } //high voltage for this sector - static void IsResolveClusters(Bool_t a) {fgIsResolveClusters=a;} + static void SetDeclustering(Bool_t a) {fgIsResolveClusters=a;} + static void SetRadioSrc(Bool_t a) {fgIsRadioSrc=a;} 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;} - inline static void Loc2Area(TVector2 x2,Int_t &padxMin,Int_t &padyMin,Int_t &padxMax,Int_t &padyMax); // - inline static Int_t Loc2Pad(TVector2 x2,Int_t &padx,Int_t &pady); //return sector and pad - inline static TVector2 Pad2Loc(Int_t padx,Int_t pady); //return center of the pad - static Int_t Sector(Int_t padx,Int_t pady) {return Pad2Sec(padx,pady);} //sector of this pad - static Int_t Sector(TVector2 x2) {int x,y;return Loc2Pad(x2,x,y);} //sector of this point + inline static TVector Loc2Area(TVector2 x2); //return area affected by hit x2 + inline static TVector Loc2Pad(TVector2 x2); //return pad containing given position + inline static TVector2 Pad2Loc(TVector pad); //return center of the pad + static TVector2 Pad2Loc(Int_t x,Int_t y) {TVector pad(2);pad[0]=x;pad[1]=y;return Pad2Loc(pad);} inline static Int_t PadNeighbours(Int_t iPadX,Int_t iPadY,Int_t aListX[4],Int_t aListY[4]); //number of neighbours for this pad - inline static TVector2 ShiftToWirePos(TVector2 x2); //shift to the nearest wire - inline static Double_t Mathieson(Double_t lx1,Double_t lx2,Double_t ly1,Double_t ly2); //Mathienson integral over these limits - inline static Double_t GainSag(Double_t y,Int_t sector); //gain variations in % - inline static Double_t QdcSlope(Int_t sec); //weight of electon in QDC channels - inline static Double_t Gain(TVector2 x2); //gain for point in ChRS - inline static Double_t FracQdc(TVector2 x2,Int_t padx,Int_t pady); //charge fraction to pad from hit + inline static Double_t Mathieson(Double_t x1,Double_t x2,Double_t y1,Double_t y2); //Mathienson integral over these limits + inline static Double_t GainSag(Double_t x,Int_t sector); //gain variations in % + static Double_t QdcSlope(Int_t sec){switch(sec){case kBad: return 0; default: return 33;}} //weight of electon in QDC channels + static Double_t Gain(TVector2 x2){if(IsWireSag()) return QdcSlope(Loc2Sec(x2))*(1+GainSag(x2.X(),Loc2Sec(x2))/100);else return QdcSlope(Loc2Sec(x2));}//gain for point in chamber RS + inline static Double_t FracQdc(TVector2 x2,TVector pad); //charge fraction to pad from hit inline static Int_t TotQdc(TVector2 x2,Double_t eloss); //total charge for hit eloss=0 for photons - inline Bool_t IsOverTh(Int_t iChamber, Int_t x, Int_t y, Double_t q); // - static Int_t NsigmaTh() {return fgNsigmaTh;} // - static Float_t SigmaThMean() {return fgSigmaThMean;} // - static Float_t SigmaThSpread() {return fgSigmaThSpread;} // - void GenSigmaThMap(); //generate pedestal map - static void Print(); -protected: + inline Bool_t IsOverTh(Int_t c,TVector pad,Double_t q); //is QDC of the pad registered by FEE + static Int_t NsigmaTh() {return fgNsigmaTh;} // + 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(Int_t &padx,Int_t &pady); //return sector, (padx,pady)->Sector RS - static Bool_t fgIsWireSag; //is wire sagitta taken into account - static Bool_t fgIsResolveClusters; //performs declustering or not - static Int_t fgHV[6]; //HV applied to anod wires - static Double_t fgAngleRot; //rotation of RICH from up postion (0,0,490)cm - static Float_t fSigmaThMap[kNCH][kNpadsX][kNpadsY]; //sigma of the pedestal distributions for all pads - static Int_t fgNsigmaTh; //n. of sigmas to cut for zero suppression - static Float_t fgSigmaThMean; //sigma threshold value - static Float_t fgSigmaThSpread; //spread of sigma - ClassDef(AliRICHParam,4) //RICH main parameters + inline static Int_t Pad2Sec(TVector pad); //return sector +protected: + TObjArray *fpChambers; //list of chambers + static Bool_t fgIsWireSag; //wire sagitta ON/OFF flag + static Bool_t fgIsResolveClusters; //declustering ON/OFF flag + static Bool_t fgIsRadioSrc; //radioactive source ON/OFF flag + static Int_t fgHV[6]; //HV applied to anod wires + static Double_t fgAngleRot; //module rotation from up postion (0,0,490)cm + static Int_t fgNsigmaTh; //n. of sigmas to cut for zero suppression + static Float_t fgSigmaThMean; //sigma threshold value + static Float_t fgSigmaThSpread; //spread of sigma + ClassDef(AliRICHParam,5) //RICH main parameters class }; //__________________________________________________________________________________________________ Int_t AliRICHParam::PadNeighbours(Int_t iPadX,Int_t iPadY,Int_t listX[4],Int_t listY[4]) { -// Determines all the neighbouring pads for the given one. Returns total amount of these pads. +// Determines all the neighbouring pads for the given one (iPadX,iPadY). Returns total number of these pads. // Dead zones are taken into account. +// 1 +// 2 3 +// 4 Int_t nPads=0; - if(iPadY!=NpadsY()&&iPadY!=NpadsYsec()) {listX[nPads]=iPadX; listY[nPads]=iPadY+1; nPads++;} - if(iPadX!=NpadsXsec()&&iPadX!=2*NpadsXsec()&&iPadX!=NpadsX()){listX[nPads]=iPadX+1; listY[nPads]=iPadY; nPads++;} - if(iPadY!=1&&iPadY!=NpadsYsec()+1) {listX[nPads]=iPadX; listY[nPads]=iPadY-1; nPads++;} - if(iPadX!=1&&iPadX!=NpadsXsec()+1&&iPadX!=2*NpadsXsec()+1) {listX[nPads]=iPadX-1; listY[nPads]=iPadY; nPads++;} + if(iPadY!=NpadsY()&&iPadY!=2*NpadsYsec()&&iPadY!=NpadsYsec()){listX[nPads]=iPadX; listY[nPads]=iPadY+1; nPads++;} //1 + if(iPadX!=1&&iPadX!=NpadsXsec()+1) {listX[nPads]=iPadX-1; listY[nPads]=iPadY; nPads++;} //2 + if(iPadX!=NpadsXsec()&&iPadX!=NpadsX()) {listX[nPads]=iPadX+1; listY[nPads]=iPadY; nPads++;} //3 + if(iPadY!=1&&iPadY!=NpadsYsec()+1&&2*NpadsYsec()+1) {listX[nPads]=iPadX; listY[nPads]=iPadY-1; nPads++;} //4 return nPads; }//Pad2ClosePads() //__________________________________________________________________________________________________ -Int_t AliRICHParam::Loc2Sec(TVector2 &x2) +Int_t AliRICHParam::Loc2Sec(TVector2 &v2) { // Determines sector containing the given point and trasform this point to the local system of that sector. -// Returns sector code: 1 2 3 -// 4 5 6 +// Returns sector code: +//y ^ 5 6 +// | 3 4 +// | 1 2 +// -------> x + Double_t x0=0; Double_t x1=SectorSizeX(); Double_t x2=SectorSizeX()+DeadZone(); Double_t x3=PcSizeX(); + Double_t y0=0; Double_t y1=SectorSizeY(); Double_t y2=SectorSizeY()+DeadZone(); Double_t y3=2*SectorSizeY()+DeadZone(); + Double_t y4=PcSizeY()-SectorSizeY(); Double_t y5=PcSizeY(); + Int_t sector=kBad; - Double_t p1=-0.5*PcSizeX(); Double_t p2=-0.5*SectorSizeX()-DeadZone(); Double_t p3=-0.5*SectorSizeX(); - Double_t p4= 0.5*SectorSizeX(); Double_t p5= 0.5*SectorSizeX()+DeadZone(); Double_t p6= 0.5*PcSizeX(); - Double_t x,y; - if (x2.X()>=p1&&x2.X()<=p2) {sector=1;x=x2.X()+0.5*PcSizeX();} - else if(x2.X()>=p3&&x2.X()<=p4) {sector=2;x=x2.X()+0.5*SectorSizeX();} - else if(x2.X()>=p5&&x2.X()<=p6) {sector=3;x=x2.X()-0.5*SectorSizeX()-DeadZone();} - else {return kBad;} //in dead zone or out of chamber + 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;} - if (x2.Y()>=-0.5*PcSizeY() &&x2.Y()<=-0.5*DeadZone()) {y=x2.Y()+0.5*PcSizeY();sector+=3;} //sectors 4,5,6 - else if(x2.Y()> -0.5*DeadZone()&&x2.Y()< 0.5*DeadZone()) {return kBad;} //in dead zone - else if(x2.Y()>= 0.5*DeadZone()&&x2.Y()<= 0.5*PcSizeY()) {y=x2.Y()-0.5*DeadZone();} //sectors 1,2,3 - else {return kBad;} //out of chamber - x2.Set(x,y); + 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;} + v2.Set(x,y); return sector; }//Loc2Sec(Double_t x, Double_t y) //__________________________________________________________________________________________________ -Int_t AliRICHParam::Loc2Pad(TVector2 x2,Int_t &padx,Int_t &pady) +TVector AliRICHParam::Loc2Pad(TVector2 x2) { -// Determines pad number (padx,pady) containing the given point x2 defined the chamber RS. +// Determines pad number TVector(padx,pady) containing the given point x2 defined the chamber RS. // Pad count starts in lower left corner from 1,1 to 144,160 in upper right corner of a chamber. // Returns sector number of the determined pad. +//y ^ 5 6 +// | 3 4 +// | 1 2 +// -------> x + TVector pad(2); Int_t sector=Loc2Sec(x2);//trasforms x2 to sector reference system - if(sector==kBad) {padx=pady=kBad; return sector;} + if(sector==kBad) {pad[0]=pad[1]=kBad; return pad;} - padx=Int_t(x2.X()/PadSizeX())+1; if(padx>NpadsXsec()) padx= NpadsXsec(); - if(sector==2||sector==5) padx+= NpadsXsec(); // 1 2 3 - if(sector==3||sector==6) padx+=2*NpadsXsec(); // 4 5 6 + pad[0]=Int_t(x2.X()/PadSizeX())+1; if(pad[0]>NpadsXsec()) pad[0]= NpadsXsec(); + if(sector==2||sector==4||sector==6) pad[0]+= NpadsXsec(); - pady=Int_t(x2.Y()/PadSizeY())+1; if(pady>NpadsYsec()) pady= NpadsYsec(); - if(sector<4) pady+=NpadsYsec(); - return sector; + pad[1]=Int_t(x2.Y()/PadSizeY())+1; if(pad[1]>NpadsYsec()) pad[1]= NpadsYsec(); + if(sector==3||sector==4) pad[1]+=NpadsYsec(); + if(sector==5||sector==6) pad[1]+=2*NpadsYsec(); + return pad; } //__________________________________________________________________________________________________ -Int_t AliRICHParam::Pad2Sec(Int_t &padx, Int_t &pady) +Int_t AliRICHParam::Pad2Sec(TVector pad) { -// Determines sector containing the given pad (padx,pady) and trasform it to the local RS of that sector. +// Determines sector containing the given pad. Int_t sector=kBad; - if (padx>=1 &&padx<=NpadsXsec()) {sector=1;} - else if(padx> NpadsXsec() &&padx<=NpadsXsec()*2) {sector=2;padx-=NpadsXsec();} - else if(padx> NpadsXsec()*2&&padx<=NpadsX()) {sector=3;padx-=NpadsXsec()*2;} - else {return 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]); + + 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]); - if (pady>=1 &&pady<=NpadsYsec()) {return sector+3;} - else if(pady>NpadsYsec() &&pady<=NpadsY()) {pady-=NpadsYsec();return sector;} - else {return kBad;} + return sector; }//Pad2Sec() //__________________________________________________________________________________________________ -TVector2 AliRICHParam::Pad2Loc(Int_t padx,Int_t pady) +TVector2 AliRICHParam::Pad2Loc(TVector pad) { -// Returns position of the center of the given pad (padx,pady) in local RS of the chamber - Int_t sector=Pad2Sec(padx,pady);//shifts to sector RS - if(sector==kBad) return TVector2(-101,-101); - Double_t x,y; - if(sector<=3) - y=0.5*DeadZone()+pady*PadSizeY()-0.5*PadSizeY(); // 1 2 3 - else{ // 4 5 6 - y=-0.5*PcSizeY()+pady*PadSizeY()-0.5*PadSizeY(); - } - if(sector==1||sector==4) - x=-0.5*PcSizeX()+padx*PadSizeX()-0.5*PadSizeX(); - else if(sector==2||sector==5) - x=-0.5*SectorSizeX()+padx*PadSizeX()-0.5*PadSizeX(); +// Returns position of the center of the given pad in local system of the chamber +// y ^ 5 6 +// | 3 4 chamber structure +// | 1 2 +// -------> x + Double_t x=kBad,y=kBad; + if(pad[0] > 0 && pad[0] <= NpadsXsec())//it's 1 or 3 or 5 + x=(pad[0]-0.5)*PadSizeX(); + else if(pad[0] > NpadsXsec() && pad[0] <= NpadsX())//it's 2 or 4 or 6 + x=(pad[0]-0.5)*PadSizeX()+DeadZone(); else - x= 0.5*SectorSizeX()+DeadZone()+padx*PadSizeX()-0.5*PadSizeX(); + ::Error("Pad2Loc","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(); + else if(pad[1] > NpadsYsec() && pad[1] <= 2*NpadsYsec())//it's 3 or 4 + y=(pad[1]-0.5)*PadSizeY()+DeadZone(); + 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]); + return TVector2(x,y); } //__________________________________________________________________________________________________ -Double_t AliRICHParam::GainSag(Double_t y,Int_t sector) +Double_t AliRICHParam::GainSag(Double_t x,Int_t sector) { // Returns % of gain variation due to wire sagita. -// All cureves are parametrized per sector basis, so y must be scaled to the Sector RS. - if(y>0) y-=SectorSizeY()/2; else y+=SectorSizeY()/2; +// All curves are parametrized as per sector basis, so y must be apriory transformed to the Sector RS. + x-=SectorSizeX()/2; + if(x>SectorSizeX()) x-=SectorSizeX(); switch(HV(sector)){ - case 2150: return 9e-6*TMath::Power(y,4)+2e-7*TMath::Power(y,3)-0.0316*TMath::Power(y,2)-3e-4*y+25.367;//% - case 2100: return 8e-6*TMath::Power(y,4)+2e-7*TMath::Power(y,3)-0.0283*TMath::Power(y,2)-2e-4*y+23.015; - case 2050: return 7e-6*TMath::Power(y,4)+1e-7*TMath::Power(y,3)-0.0254*TMath::Power(y,2)-2e-4*y+20.888; - case 2000: return 6e-6*TMath::Power(y,4)+8e-8*TMath::Power(y,3)-0.0227*TMath::Power(y,2)-1e-4*y+18.961; + case 2150: return 9e-6*TMath::Power(x,4)+2e-7*TMath::Power(x,3)-0.0316*TMath::Power(x,2)-3e-4*x+25.367;//% + case 2100: return 8e-6*TMath::Power(x,4)+2e-7*TMath::Power(x,3)-0.0283*TMath::Power(x,2)-2e-4*x+23.015; + case 2050: return 7e-6*TMath::Power(x,4)+1e-7*TMath::Power(x,3)-0.0254*TMath::Power(x,2)-2e-4*x+20.888; + case 2000: return 6e-6*TMath::Power(x,4)+8e-8*TMath::Power(x,3)-0.0227*TMath::Power(x,2)-1e-4*x+18.961; default: return 0; } } //__________________________________________________________________________________________________ -Double_t AliRICHParam::QdcSlope(Int_t sec) -{ -// Returns number of QDC channels per single electron at the unknown yet ???? point for a given sector - switch(sec){ - case kBad: return 0; - default: return 27; - } -} -//__________________________________________________________________________________________________ -Double_t AliRICHParam::Gain(TVector2 x2) -{ -// - if(IsWireSag()) - return QdcSlope(Sector(x2))*(1+GainSag(x2.Y(),Sector(x2))/100); - else - return QdcSlope(Sector(x2)); -} -//__________________________________________________________________________________________________ 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. +// 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 - if(Sector(x2)==kBad) return 0; //hit in the dead zone + 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; for(Int_t i=1;i<=iNelectrons;i++) qdc+=-Gain(x2)*TMath::Log(gRandom->Rndm()); return Int_t(qdc); } //__________________________________________________________________________________________________ -Double_t AliRICHParam::FracQdc(TVector2 x2,Int_t padx,Int_t pady) +Double_t AliRICHParam::FracQdc(TVector2 x2,TVector pad) { -// Calculates the charge fraction for a given pad (padx,pady) from the given hit point. -// Mathieson distribution integrated is used. - TVector2 center2=Pad2Loc(padx,pady);//gives center of requested 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(Sector(x2)!=Sector(padx,pady)) return 0;//requested pad does not belong to the sector of given point - else return Mathieson(normXmin, normYmin, normXmax, normYmax); + 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 AliRICHParam::Mathieson(Double_t xMin,Double_t yMin,Double_t xMax,Double_t yMax) @@ -268,32 +287,40 @@ Double_t AliRICHParam::Mathieson(Double_t xMin,Double_t yMin,Double_t xMax,Doubl return 4*kX4*(TMath::ATan(ux2)-TMath::ATan(ux1))*kY4*(TMath::ATan(uy2)-TMath::ATan(uy1)); } //__________________________________________________________________________________________________ -void AliRICHParam::Loc2Area(TVector2 x2,Int_t &iPadXmin,Int_t &iPadYmin,Int_t &iPadXmax,Int_t &iPadYmax) +TVector AliRICHParam::Loc2Area(TVector2 x2) { // Calculates the area of disintegration for a given point. It's assumed here that this points lays on anode wire. // Area is a rectangulare set of pads defined by its left-down and right-up coners. - Loc2Pad(x2-MathiesonDelta(),iPadXmin,iPadYmin); - Loc2Pad(x2+MathiesonDelta(),iPadXmax,iPadYmax); + TVector area(4); + TVector pad=Loc2Pad(x2); + area[0]=area[2]=pad[0]; area[1]=area[3]=pad[1];//area is just a pad fired + if(pad[0]!=1 && pad[0]!= NpadsXsec()+1 ) area[0]--; //left down coner X + if(pad[1]!=1 && pad[1]!= NpadsYsec()+1 && pad[1]!= 2*NpadsYsec()+1) area[1]--; //left down coner Y + if(pad[0]!=NpadsXsec() && pad[0]!= NpadsX() ) area[2]++; //right up coner X + if(pad[1]!=NpadsYsec() && pad[1]!= 2*NpadsYsec() && pad[1]!= NpadsY() ) area[3]++; //right up coner Y + return area; } //__________________________________________________________________________________________________ -Bool_t AliRICHParam::IsOverTh(Int_t c,Int_t x,Int_t y,Double_t q) +Bool_t AliRICHParam::IsOverTh(Int_t ,TVector ,Double_t q) { -// Calculate the new charge subtracting pedestal and if the current digit is over threshold - if (c>0 && x>0 && y>0 && cNsigmaTh()*fSigmaThMap[c-1][x-1][y-1]) return kTRUE; - return kFALSE; +// 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())); } //__________________________________________________________________________________________________ -TVector2 AliRICHParam::ShiftToWirePos(TVector2 x2) +void AliRICHParam::PropogateHelix(TVector3 x0,TVector3 p0,Double_t s,TVector3 *x,TVector3 *p) { -// Calculate the position of the wire nearest to the hit - Int_t padx,pady; - Loc2Pad(x2,padx,pady); - Double_t x; - TVector2 center2=Pad2Loc(padx,pady); - if(x2.X()>center2.X()) x=center2.X()+0.5*WirePitch(); - else x=center2.X()-0.5*WirePitch(); - x2.Set(x,x2.Y()); - return x2; +// 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/AliRICHv0.cxx b/RICH/AliRICHv0.cxx index b4b9fd8ba6d..9bfebf42bae 100644 --- a/RICH/AliRICHv0.cxx +++ b/RICH/AliRICHv0.cxx @@ -22,6 +22,8 @@ #include #include +#include + ClassImp(AliRICHv0) void AliRICHv0::StepManager() @@ -93,11 +95,59 @@ void AliRICHv0::StepManager() if(gMC->VolId("CSI ")==gMC->CurrentVolID(copy0)){ Int_t iChamber; gMC->CurrentVolOffID(2,iChamber); - TVector3 x3=C(iChamber)->G2L(x4); - Info("","loc(%+8.3f,%+8.3f,%8.3f) by G2L", x3.X(),x3.Y(),x3.Z()); TVector2 x2=C(iChamber)->Glob2Loc(x4); - Info("","loc(%+8.3f,%+8.3f) by Global2Local", x2.X(),x2.Y()); + Info("","loc(%+8.3f,%+8.3f) by Glob2Loc", x2.X(),x2.Y()); } Info("","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 47f0a4813ac..3594ff85b09 100644 --- a/RICH/AliRICHv0.h +++ b/RICH/AliRICHv0.h @@ -9,12 +9,13 @@ class AliRICHv0 : public AliRICH { public: - AliRICHv0():AliRICH() {;} - AliRICHv0(const char *name, const char *title):AliRICH(name,title) {;} - virtual ~AliRICHv0() {;} - virtual void Init() {;} - virtual Int_t IsVersion() const{return 0;} - virtual void StepManager(); + AliRICHv0():AliRICH() {;} //default ctor + AliRICHv0(const char *name, const char *title):AliRICH(name,title) {;} //named ctor + virtual ~AliRICHv0() {;} //dtor + 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 a00b589e5ff..8d9efdababe 100644 --- a/RICH/AliRICHv1.cxx +++ b/RICH/AliRICHv1.cxx @@ -31,19 +31,34 @@ ClassImp(AliRICHv1) //__________________________________________________________________________________________________ void AliRICHv1::StepManager() { -//Full Step Manager - +// Full Step Manager. +// 3- Ckovs absorbed in Collection electrods +// 5- Ckovs absorbed in Cathode wires +// 6- Ckovs absorbed in 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 + } + //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(IsLostByFresnel()){ gMC->StopTrack(); return;} + if(gMC->Edep()>0){//CF+CSI+DE + if(IsLostByFresnel()){ + if(gMC->TrackPid()==kCerenkov) fCounters(7)++;// 7- Ckovs reflected on CsI + gMC->StopTrack(); + return; + } gMC->TrackPosition(cerX4); gMC->CurrentVolOffID(2,iCurrentChamber); AddHit(iCurrentChamber,gAlice->GetMCApp()->GetCurrentTrackNumber(),cerX4.Vect(),cerX4.Vect());//HIT for PHOTON in conditions CF+CSI+DE + fCounters(8)++;//4- Ckovs converted in CsI GenerateFeedbacks(iCurrentChamber); }//CF+CSI+DE }//CF in CSI @@ -68,6 +83,7 @@ void AliRICHv1::StepManager() //__________________________________________________________________________________________________ Bool_t AliRICHv1::IsLostByFresnel() { +// Calculate probability for the photon to be lost by Fresnel reflection. TLorentzVector p4; Double_t mom[3],localMom[3]; gMC->TrackMomentum(p4); mom[0]=p4(1); mom[1]=p4(2); mom[2]=p4(3); @@ -77,9 +93,84 @@ 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",""); + if(GetDebug()) Info("IsLostByFresnel","Photon lost"); return kTRUE; }else return kFALSE; }//IsLostByFresnel() //__________________________________________________________________________________________________ +void AliRICHv1::GenerateFeedbacks(Int_t iChamber,Float_t eloss) +{ +// Generate FeedBack photons for the current particle. To be invoked from StepManager(). +// eloss=0 means photon so only pulse height distribution is to be analysed. This one is done in AliRICHParam::TotQdc() + + 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 + 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 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 + for(Int_t i=0;iGetRandom()->RndmArray(2,ranf); //Sample direction + cthf=ranf[0]*2-1.0; + if(cthf<0) continue; + sthf = TMath::Sqrt((1 - cthf) * (1 + cthf)); + phif = ranf[1] * 2 * TMath::Pi(); + + if(Double_t randomNumber=gMC->GetRandom()->Rndm()<=0.57) + enfp = 7.5e-9; + else if(randomNumber<=0.7) + enfp = 6.4e-9; + else + enfp = 7.9e-9; + + + dir[0] = sthf * TMath::Sin(phif); dir[1] = cthf; dir[2] = sthf * TMath::Cos(phif); + gMC->Gdtom(dir, mom, 2); + mom[0]*=enfp; mom[1]*=enfp; mom[2]*=enfp; + mom[3] = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2]); + + // Polarisation + e1[0]= 0; e1[1]=-dir[2]; e1[2]= dir[1]; + e2[0]=-dir[1]; e2[1]= dir[0]; e2[2]= 0; + e3[0]= dir[1]; e3[1]= 0; e3[2]=-dir[0]; + + vmod=0; + for(j=0;j<3;j++) vmod+=e1[j]*e1[j]; + if (!vmod) for(j=0;j<3;j++) { + uswop=e1[j]; + e1[j]=e3[j]; + e3[j]=uswop; + } + vmod=0; + for(j=0;j<3;j++) vmod+=e2[j]*e2[j]; + if (!vmod) for(j=0;j<3;j++) { + uswop=e2[j]; + e2[j]=e3[j]; + e3[j]=uswop; + } + + vmod=0; for(j=0;j<3;j++) vmod+=e1[j]*e1[j]; vmod=TMath::Sqrt(1/vmod); for(j=0;j<3;j++) e1[j]*=vmod; + vmod=0; for(j=0;j<3;j++) vmod+=e2[j]*e2[j]; vmod=TMath::Sqrt(1/vmod); for(j=0;j<3;j++) e2[j]*=vmod; + + phi = gMC->GetRandom()->Rndm()* 2 * TMath::Pi(); + for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi); + gMC->Gdtom(pol, pol, 2); + Int_t outputNtracksStored; + gAlice->GetMCApp()->PushTrack(1, //do not transport + gAlice->GetMCApp()->GetCurrentTrackNumber(),//parent track + kFeedback, //PID + mom[0],mom[1],mom[2],mom[3], //track momentum + x4.X(),x4.Y(),x4.Z(),x4.T(), //track origin + pol[0],pol[1],pol[2], //polarization + kPFeedBackPhoton, + outputNtracksStored, + 1.0); + }//feedbacks loop + if(GetDebug()) Info("GenerateFeedbacks","Stop."); +}//GenerateFeedbacks() diff --git a/RICH/AliRICHv1.h b/RICH/AliRICHv1.h index f848de9bff1..41a1ff6d467 100644 --- a/RICH/AliRICHv1.h +++ b/RICH/AliRICHv1.h @@ -15,10 +15,11 @@ public: virtual void Init() {;} virtual Int_t IsVersion() const{return 1;} - virtual void StepManager(); - Bool_t IsLostByFresnel(); + virtual void StepManager(); //full slow step manager + Bool_t IsLostByFresnel(); //checks if the photon lost on Fresnel reflection + void GenerateFeedbacks(Int_t iChamber,Float_t eloss=0); //generates feedback photons; eloss=0 for photon private: - ClassDef(AliRICHv1,1)//RICH full version for simulation + ClassDef(AliRICHv1,1) //RICH full version for simulation }; #endif diff --git a/RICH/RICHLinkDef.h b/RICH/RICHLinkDef.h index 4768f4892aa..87321e3505b 100644 --- a/RICH/RICHLinkDef.h +++ b/RICH/RICHLinkDef.h @@ -9,7 +9,6 @@ #pragma link C++ class AliRICHhit+; #pragma link C++ class AliRICHdigit+; #pragma link C++ class AliRICHcluster+; -#pragma link C++ class AliRICHreco+; #pragma link C++ class AliRICHChamber+; #pragma link C++ class AliRICHMap+; #pragma link C++ class AliRICHClusterFinder+; -- 2.39.3