From d3eb6079b2ff0d28ed93f221ab88b7ce58552ecc Mon Sep 17 00:00:00 2001 From: kir Date: Wed, 2 Mar 2005 06:26:55 +0000 Subject: [PATCH] RAW DDL first attempt, also virtual hidded problem addressed --- RICH/AliRICH.cxx | 240 ++++++++++++++++++++++------------ RICH/AliRICH.h | 28 ++-- RICH/AliRICHChamber.cxx | 14 +- RICH/AliRICHChamber.h | 59 +++++---- RICH/AliRICHClusterFinder.cxx | 4 +- RICH/AliRICHClusterFinder.h | 2 +- RICH/AliRICHDisplFast.cxx | 55 ++++---- RICH/AliRICHDisplFast.h | 16 +-- RICH/AliRICHHelix.cxx | 16 ++- RICH/AliRICHHelix.h | 67 +++++----- RICH/AliRICHMap.h | 4 +- RICH/AliRICHParam.cxx | 163 ++++++++++------------- RICH/AliRICHParam.h | 14 +- RICH/AliRICHReconstructor.h | 2 + RICH/AliRICHTracker.cxx | 41 +++--- RICH/RichMake | 10 +- RICH/RichMenu.C | 238 +++++++++++++++++++-------------- RICH/api.txt | 58 ++++---- RICH/libRICHsim.pkg | 2 +- 19 files changed, 583 insertions(+), 450 deletions(-) diff --git a/RICH/AliRICH.cxx b/RICH/AliRICH.cxx index 7601834d333..73ffe932e5b 100644 --- a/RICH/AliRICH.cxx +++ b/RICH/AliRICH.cxx @@ -16,8 +16,8 @@ #include "AliRICH.h" #include "AliRICHParam.h" #include "AliRICHChamber.h" -#include "AliRICHClusterFinder.h" -#include +#include "AliRICHTracker.h" +//#include #include #include #include @@ -34,21 +34,28 @@ #include #include #include -#include +#include #include +#include +#include +#include +#include + ClassImp(AliRICHhit) //__________________________________________________________________________________________________ void AliRICHhit::Print(Option_t*)const { - AliInfo(Form("Ch=%1i, TID=%6i, eloss=%9.3f eV, in-out dist=%9.4f, OUT(%7.2f,%7.2f,%7.2f)" - ,fChamber,fTrack,fEloss*1e9,Length(),fOutX3.X(),fOutX3.Y(),fOutX3.Z())); + AliInfo(Form("Ch=%1i,TID=%6i,Elos=%9.3f eV,IN(%6.2f,%6.2f,%6.2f)-OUT(%6.2f,%6.2f,%6.2f)=%9.4f" + ,fChamber,fTrack,fEloss*1e9,fInX3.X() ,fInX3.Y() ,fInX3.Z(), + fOutX3.X(),fOutX3.Y(),fOutX3.Z(),Length())); } //__________________________________________________________________________________________________ ClassImp(AliRICHdigit) //__________________________________________________________________________________________________ void AliRICHdigit::Print(Option_t*)const { +//Print current digit AliInfo(Form("cfm=%9i, cs=%2i, x=%3i, y=%3i, q=%8.3f, TID1=%5i, TID2=%5i, TID3=%5i", fCFM,fChamber,fPadX,fPadY,fQdc,fTracks[0],fTracks[1],fTracks[2])); } @@ -57,6 +64,7 @@ ClassImp(AliRICHcluster) //__________________________________________________________________________________________________ void AliRICHcluster::Print(Option_t*)const { +//Print current cluster const char *status=0; switch(fStatus){ case kRaw: status="raw" ;break; @@ -124,14 +132,14 @@ void AliRICH::Hits2SDigits() GetLoader()->GetRunLoader()->GetEvent(iEventN);//get next event if(!GetLoader()->TreeH()) GetLoader()->LoadHits(); GetLoader()->GetRunLoader()->LoadHeader(); - GetLoader()->GetRunLoader()->LoadKinematics();//from + if(!GetLoader()->GetRunLoader()->TreeK()) GetLoader()->GetRunLoader()->LoadKinematics();//from if(!GetLoader()->TreeS()) GetLoader()->MakeTree("S"); MakeBranch("S");//to for(Int_t iPrimN=0;iPrimNTreeH()->GetEntries();iPrimN++){//prims loop GetLoader()->TreeH()->GetEntry(iPrimN); for(Int_t iHitN=0;iHitNGetEntries();iHitN++){//hits loop AliRICHhit *pHit=(AliRICHhit*)Hits()->At(iHitN);//get current hit - TVector2 x2 = C(pHit->C())->Mrs2Pc(0.5*(pHit->InX3()+pHit->OutX3()));//hit position in the anod plane + TVector2 x2 = C(pHit->C())->Mrs2Anod(0.5*(pHit->InX3()+pHit->OutX3()));//hit position in the anod plane Int_t iTotQdc=P()->TotQdc(x2,pHit->Eloss());//total charge produced by hit, 0 if hit in dead zone if(iTotQdc==0) continue; // @@ -163,6 +171,51 @@ void AliRICH::Hits2SDigits() AliDebug(1,"Stop."); }//Hits2SDigits() //__________________________________________________________________________________________________ +void AliRICH::Digits2Raw() +{ +//Loops over all digits and creates raw data files in DDL format. GetEvent() is done outside (AliSimulation) + AliDebug(1,"Start."); + GetLoader()->LoadDigits(); + GetLoader()->TreeD()->GetEntry(0); + + Int_t kRichOffset=0x700; //currently one DDL per 3 sectors + + ofstream file135,file246;//output streams 2 DDL per chamber + AliRawDataHeader header;//empty DDL miniheader + UInt_t word32=1; //32 bits data word + + for(Int_t iChN=1;iChN<=kNchambers;iChN++){ //2 DDL per chamber open both in parallel + file135.open(Form("RICH_%4i.ddl",kRichOffset+2*iChN-1)); //left part of chamber; sectors 1-3-5 + file246.open(Form("RICH_%4i.ddl",kRichOffset+2*iChN)); //right part of chamber; sectors 2-4-6 +//common DDL header defined by standart, now dummy as the total number of bytes is not yet known + file135.write((char*)&header,sizeof(header)); //dummy header just place holder + file246.write((char*)&header,sizeof(header)); //actual will be written later +//now coming RICH private header 16 32-bits wors currently signifying nothing + for(Int_t i=0;i<16;i++){ + word32='t'; + file135.write((char*)&word32,sizeof(word32)); //dummy header just place holder + file246.write((char*)&word32,sizeof(word32)); //actual will be written later + } + + Int_t counter135=0,counter246=0;//counts total number of records per DDL + + for(Int_t iDigN=0;iDigNGetEntriesFast();iDigN++){//digits loop for a given chamber + AliRICHdigit *pDig=(AliRICHdigit*)Digits(iChN)->At(iDigN); + word32=UInt_t (pDig->Q()+0x400*pDig->X()+0x4000*pDig->Y()); //arbitrary structure + switch(pDig->S()){//readout by vertical sectors: 1,3,5-left DDL 2,4,6-right DDL + case 1: case 3: case 5: file135.write((char*)&word32,sizeof(word32)); counter135++; break; + case 2: case 4: case 6: file246.write((char*)&word32,sizeof(word32)); counter246++; break; + }//switch sector + }//digits loop for a given chamber +//now count total byte number for each DDL file and rewrite actual header + header.fSize=sizeof(header)+counter135*sizeof(word32); header.SetAttribute(0); file135.seekp(0); file135.write((char*)&header,sizeof(header)); + header.fSize=sizeof(header)+counter246*sizeof(word32); header.SetAttribute(0); file246.seekp(0); file246.write((char*)&header,sizeof(header)); + file135.close(); file246.close(); + }//chambers loop + GetLoader()->UnloadDigits(); + AliDebug(1,"Stop."); +}//Digits2Raw() +//__________________________________________________________________________________________________ void AliRICH::BuildGeometry() { //Builds a TNode geometry for event display @@ -441,38 +494,48 @@ void AliRICH::Print(Option_t *option)const //__________________________________________________________________________________________________ void AliRICH::ControlPlots() { -// Creates a set of hists to control the results of simulation. Hists are in file $HOME/RCP.root - +//Creates a set of QA hists to control the results of simulation. Hists are in file $HOME/RCP.root TH1F *pElecP=0 ,*pMuonP=0 ,*pPionP=0 ,*pKaonP=0 ,*pProtP=0, //stack particles - *pHxD=0,*pHyD=0,*pNumClusH1=0, + *pHxD=0,*pHyD=0,*pHxSd=0,*pHySd=0, //diff hit position - digit sdigit position + *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; - - GetLoader()->GetRunLoader()->LoadHeader(); - GetLoader()->GetRunLoader()->LoadKinematics(); - - Bool_t isDig =!GetLoader()->LoadDigits(); + TH1F *pelecRadius=0,*pprotRadius=0,*pprotbarRadius=0; +//load all information + GetLoader()->GetRunLoader()->LoadHeader(); + GetLoader()->GetRunLoader()->LoadKinematics(); + GetLoader()->LoadHits(); + Bool_t isSdig=0;//!GetLoader()->LoadSDigits(); + Bool_t isDig =0;//!GetLoader()->LoadDigits(); Bool_t isClus=!GetLoader()->LoadRecPoints(); -// if(!isDig && !isClus){AliError("No digits and clusters! Nothing to do.");return;} - - TStopwatch sw;TDatime time; + gBenchmark->Start("ControlPlots"); TFile *pFile = new TFile("$(HOME)/RCP.root","RECREATE"); - pElecP=new TH1F("Pelec","e versus momentum;P [GeV]",1000,-10,10); - pMuonP=new TH1F("Pmuon","mu versus momentum;P [GeV]",1000,-10,10); - pPionP=new TH1F("Ppion","pi versus momentum;P [GeV]",1000,-10,10); - pKaonP=new TH1F("Pkaon","K versus momentum;P [GeV]",1000,-10,10); - pProtP=new TH1F("Pprot","p versus momentum;P [GeV]",1000,-10,10); + + pElecP=new TH1F("Pelec","Electrons made hit in RICH;p [GeV]",1000,-30,30); + pMuonP=new TH1F("Pmuon","Muons made hit in RICH;p [GeV]",1000,-30,30); + pPionP=new TH1F("Ppion","Pions made hit in RICH;p [GeV]",1000,-30,30); + pKaonP=new TH1F("Pkaon","Kaon made hit in RICH;p [GeV]",1000,-30,30); + pProtP=new TH1F("Pprot","Protons made hit in RICH;p [GeV]",1000,-30,30); + pelecRadius=new TH1F("elecRadius","elec",600,0.,600.); + pprotRadius=new TH1F("protRadius","elec",600,0.,600.); + pprotbarRadius=new TH1F("protbarRadius","elec",600,0.,600.); + + if(isSdig){ + AliInfo("SDigits available"); + pHxSd=new TH1F("DiffHitSDigitX","Hit-SDigit diff X all chambers;diff [cm]",300,-10,10); + pHySd=new TH1F("DiffHitSDigitY","Hit-SDigit diff Y all chambers;diff [cm]",300,-10,10); + }//isSdig if(isDig){ AliInfo("Digits available"); - pHxD=new TH1F("HitDigitDiffX","Hit-Digits diff X all chambers;diff [cm]",100,-10,10); - pHyD=new TH1F("HitDigitDiffY","Hit-Digits diff Y all chambers;diff [cm]",100,-10,10); + pHxD=new TH1F("DiffHitDigitX","Hit-Digit diff X all chambers;diff [cm]",300,-10,10); + pHyD=new TH1F("DiffHitDigitY","Hit-Digit diff Y all chambers;diff [cm]",300,-10,10); }//isDig if(isClus){ @@ -497,20 +560,21 @@ void AliRICH::ControlPlots() pPureFeeQdcH1 =new TH1F("QdcPureFee" ,"Feedback only Cluster Charge all chambers;q [QDC]",P()->MaxQdc(),0,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,P()->PcSizeX(),1000,0,P()->PcSizeY()); + }//isClus - +//end of hists booking for(Int_t iEvtN=0;iEvtN < GetLoader()->GetRunLoader()->GetAliRun()->GetEventsPerRun();iEvtN++){//events loop - GetLoader()->GetRunLoader()->GetEvent(iEvtN); //gets current event + GetLoader()->GetRunLoader()->GetEvent(iEvtN); //get current event - if(!GetLoader()->TreeH()) GetLoader()->LoadHits(); for(Int_t iPrimN=0;iPrimN < GetLoader()->TreeH()->GetEntries();iPrimN++){//hit tree loop GetLoader()->TreeH()->GetEntry(iPrimN); - for(Int_t j=0;jGetEntries();j++){//hits loop for a given primary - AliRICHhit *pHit = (AliRICHhit*)Hits()->At(j); - TParticle *pParticle = GetLoader()->GetRunLoader()->Stack()->Particle(pHit->GetTrack()); + for(Int_t j=0;jGetEntries();j++){//hits loop + AliRICHhit *pHit = (AliRICHhit*)Hits()->At(j); + TParticle *pParticle = GetLoader()->GetRunLoader()->Stack()->Particle(pHit->GetTrack());//get particle produced this hit + Double_t Radius = TMath::Sqrt(pParticle->Vx()*pParticle->Vx()+pParticle->Vy()*pParticle->Vy()+pParticle->Vz()*pParticle->Vz()); switch(pParticle->GetPdgCode()){ - case kPositron : pElecP->Fill( pParticle->P()); break; - case kElectron : pElecP->Fill(-pParticle->P()); break; + case kPositron : pElecP->Fill( pParticle->P());pelecRadius->Fill(Radius); break; + case kElectron : pElecP->Fill(-pParticle->P());pelecRadius->Fill(Radius); break; case kMuonPlus : pMuonP->Fill( pParticle->P()); break; case kMuonMinus: pMuonP->Fill(-pParticle->P()); break; @@ -521,18 +585,39 @@ void AliRICH::ControlPlots() case kKPlus : pKaonP->Fill( pParticle->P()); break; case kKMinus : pKaonP->Fill(-pParticle->P()); break; - case kProton : pProtP->Fill( pParticle->P()); break; - case kProtonBar: pProtP->Fill(-pParticle->P()); break; + case kProton : pProtP->Fill( pParticle->P()); pprotRadius->Fill(Radius); break; + case kProtonBar: pProtP->Fill(-pParticle->P()); pprotbarRadius->Fill(Radius); break; }//switch PdgCode }//hits loop }//hit tree loop - if(isClus) GetLoader()->TreeR()->GetEntry(0); + if(isSdig){ + GetLoader()->TreeS()->GetEntry(0); + for(Int_t iSdigN=0;iSdigNGetEntries();iSdigN++){//sdigits loop + AliRICHdigit *pSdig=(AliRICHdigit*)SDigits()->At(iSdigN); //get current sdigit pointer + AliRICHhit *pHit=Hit(pSdig->GetTrack(0)); //get hit of this sdigit (always one) + TVector2 hit2 =C(pHit->C())->Mrs2Pc(pHit->OutX3()); //this hit position in local system + TVector2 sdig2=P()->Pad2Loc(pSdig->Pad()); //center of pad for this sdigit + pHxSd->Fill(hit2.X()-sdig2.X()); + pHySd->Fill(hit2.Y()-sdig2.Y()); + }//sdigits loop + }//if(isSdig) + if(isDig) GetLoader()->TreeD()->GetEntry(0); + if(isClus) GetLoader()->TreeR()->GetEntry(0); for(Int_t iChamN=1;iChamN<=7;iChamN++){//chambers loop + if(isDig){ + for(Int_t iDigN=0;iDigNGetEntries();iDigN++){//digits loop + AliRICHdigit *pDig=(AliRICHdigit*)Digits(iChamN)->At(iDigN); + AliRICHhit *pHit=Hit(pDig->GetTrack(0)); //get first hit of this digit + TVector2 hitV2=C(iChamN)->Mrs2Pc(pHit->OutX3()); + TVector2 digV2=P()->Pad2Loc(pDig->Pad()); //center of pad for this digit + pHxD->Fill(hitV2.X()-digV2.X()); pHyD->Fill(hitV2.Y()-digV2.Y()); + }//digits loop + }//isDig if(isClus){ Int_t iNclusCham=Clusters(iChamN)->GetEntries(); if(iNclusCham) pNumClusH1->Fill(iNclusCham);//number of clusters per event for(Int_t iClusN=0;iClusNIsPureMip()) pPhotQdcH1->Fill(pClus->Q()); //not MIP }//clusters loop }//isClus - if(isDig){ - for(Int_t iDigN=0;iDigNGetEntries();iDigN++){//digits loop - AliRICHdigit *pDig=(AliRICHdigit*)Digits(iChamN)->At(iDigN); - AliRICHhit *pHit=Hit(pDig->GetTrack(0));//get first hit of this digit - TVector2 hitV2=C(iChamN)->Mrs2Pc(pHit->OutX3()); TVector2 digV2=P()->Pad2Loc(pDig->Pad());//center of pad for digit - 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 - + GetLoader()->UnloadHits(); + if(isSdig) GetLoader()->UnloadSDigits(); if(isDig) GetLoader()->UnloadDigits(); if(isClus) GetLoader()->UnloadRecPoints(); @@ -576,10 +654,11 @@ void AliRICH::ControlPlots() GetLoader()->GetRunLoader()->UnloadKinematics(); pFile->Write(); delete pFile; - sw.Print();time.Print(); + + gBenchmark->Show("ControlPlots"); }//ControlPlots() //__________________________________________________________________________________________________ -AliRICHhit* AliRICH::Hit(Int_t tid) +AliRICHhit* AliRICH::Hit(Int_t tid)const { //defines which hit provided by given tid for the currently loaded event GetLoader()->LoadHits(); @@ -597,8 +676,8 @@ AliRICHhit* AliRICH::Hit(Int_t tid) void AliRICH::PrintHits(Int_t iEvtN) { //Prints a list of RICH hits for a given event. Default is event number 0. + if(GetLoader()->GetRunLoader()->GetEvent(iEvtN)) return; AliInfo(Form("List of RICH hits for event %i",iEvtN)); - GetLoader()->GetRunLoader()->GetEvent(iEvtN); if(GetLoader()->LoadHits()) return; Int_t iTotalHits=0; @@ -608,14 +687,15 @@ void AliRICH::PrintHits(Int_t iEvtN) iTotalHits+=Hits()->GetEntries(); } GetLoader()->UnloadHits(); + ResetHits(); AliInfo(Form("totally %i hits",iTotalHits)); } //__________________________________________________________________________________________________ void AliRICH::PrintSDigits(Int_t iEvtN) { //prints a list of RICH sdigits for a given event + if(GetLoader()->GetRunLoader()->GetEvent(iEvtN)) return; Info("PrintSDigits","List of RICH sdigits for event %i",iEvtN); - GetLoader()->GetRunLoader()->GetEvent(iEvtN); if(GetLoader()->LoadSDigits()) return; GetLoader()->TreeS()->GetEntry(0); @@ -627,8 +707,8 @@ void AliRICH::PrintSDigits(Int_t iEvtN) void AliRICH::PrintDigits(Int_t iEvtN) { //prints a list of RICH digits for a given event + if(GetLoader()->GetRunLoader()->GetEvent(iEvtN)) return; Info("PrintDigits","List of RICH digits for event %i",iEvtN); - GetLoader()->GetRunLoader()->GetEvent(iEvtN); if(GetLoader()->LoadDigits()) return; Int_t iTotalDigits=0; @@ -644,7 +724,7 @@ void AliRICH::PrintDigits(Int_t iEvtN) 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); + AliInfo(Form("List of RICH clusters for event %i",iEvtN)); GetLoader()->GetRunLoader()->GetEvent(iEvtN); if(GetLoader()->LoadRecPoints()) return; @@ -655,7 +735,7 @@ void AliRICH::PrintClusters(Int_t iEvtN) iTotalClusters+=Clusters(iChamber)->GetEntries(); } GetLoader()->UnloadRecPoints(); - Info("PrintClusters","totally %i clusters",iTotalClusters); + AliInfo(Form("totally %i clusters for event %i",iTotalClusters,iEvtN)); } //__________________________________________________________________________________________________ void AliRICH::PrintTracks(Int_t iEvtN) @@ -667,35 +747,14 @@ void AliRICH::PrintTracks(Int_t iEvtN) if(GetLoader()->GetRunLoader()->LoadKinematics()) return; AliStack *pStack=GetLoader()->GetRunLoader()->Stack(); - for(Int_t i=0;iGetNtrack();i++){ - pStack->Particle(i)->Print(); - } + for(Int_t i=0;iGetNtrack();i++) pStack->Particle(i)->Print(); - AliInfo(Form("totally %i tracks including %i primaries",pStack->GetNtrack(),pStack->GetNprimary())); + AliInfo(Form("totally %i tracks including %i primaries for event %i",pStack->GetNtrack(),pStack->GetNprimary(),iEvtN)); GetLoader()->GetRunLoader()->UnloadHeader(); GetLoader()->GetRunLoader()->UnloadKinematics(); } //__________________________________________________________________________________________________ -Int_t AliRICH::Nparticles(Int_t iPartID,Int_t iEvtN,AliRunLoader *pRL) -{ -//counts total number of particles of given type (including secondary) for a given event - pRL->GetEvent(iEvtN); - if(pRL->LoadHeader()) return 0; - if(pRL->LoadKinematics()) return 0; - AliStack *pStack=pRL->Stack(); - - Int_t iCounter=0; - for(Int_t i=0;iGetNtrack();i++){ - if(pStack->Particle(i)->GetPdgCode()==iPartID) iCounter++; - } - - pRL->UnloadHeader(); - pRL->UnloadKinematics(); - return iCounter; -} - -//__________________________________________________________________________________________________ -void AliRICH::GeomPadPanelFrame() +void AliRICH::GeomPadPanelFrame()const { //Pad Panel frame 6 sectors Double_t cm=1,mm=0.1*cm;//default is cm @@ -731,7 +790,7 @@ void AliRICH::GeomPadPanelFrame() gMC->Gspos("PPFL",8,"RPPF", +224.5*mm, +151.875*mm, 0.85*mm, 0,"ONLY"); }//GeomPadPanelFrame() //__________________________________________________________________________________________________ -void AliRICH::GeomAmpGap() +void AliRICH::GeomAmpGap()const { //Gap - anod wires 6 copies to RICH Double_t cm=1,mm=0.1*cm,mkm=0.001*mm;//default is cm @@ -752,7 +811,7 @@ void AliRICH::GeomAmpGap() gMC->Gspos("RANO",i,"RGAP", 0*mm, -411/2*mm+i*4*mm, 0.185*mm, matrixIdReturn,"ONLY"); //WP 2099P1 }//GeomAmpGap() //__________________________________________________________________________________________________ -void AliRICH::GeomRadiators() +void AliRICH::GeomRadiators()const { //Defines radiators geometry Double_t mm=0.1;//default is cm @@ -778,7 +837,7 @@ void AliRICH::GeomRadiators() gMC->Gspos("RRSP",10*i+j,"RRAD",-1330*mm/2+116*mm+j*122*mm,(i-1)*105*mm,-0.5*mm,0,"ONLY");//spacers }//GeomRadiators() //__________________________________________________________________________________________________ -void AliRICH::GeomSandBox() +void AliRICH::GeomSandBox()const { //Defines SandBox geometry Double_t mm=0.1;//default is cm @@ -793,7 +852,7 @@ void AliRICH::GeomSandBox() gMC->Gspos("RSCO",2,"RSNB", 0*mm, 0*mm, -25*mm, 0,"ONLY"); //cover to sandbox }//GeomSandBox() //__________________________________________________________________________________________________ -void AliRICH::GeomRadioSrc() +void AliRICH::GeomRadioSrc()const { // Defines geometry for radioactive source Double_t cm=1,mm=0.1*cm,mkm=0.001*cm; @@ -820,7 +879,7 @@ void AliRICH::GeomRadioSrc() gMC->Gspos("RSSR",1,"RSSS", 0, 0, -4.5*mm, 0,"ONLY");//Sr90 in support steel screw }//GeomSr90() //__________________________________________________________________________________________________ -void AliRICH::GeomAerogel() +void AliRICH::GeomAerogel()const { //Creates detailed geometry for aerogel study. AliDebug(1,"Start."); @@ -869,3 +928,22 @@ void AliRICH::CreateGeometry() AliDebug(1,"Stop main."); }//CreateGeometry() //__________________________________________________________________________________________________ +void AliRICH::CheckPR()const +{ +//Pattern recognition with stack particles + TFile *pFile = new TFile("$(HOME)/RPR.root","RECREATE","RICH Pattern Recognition"); + TNtupleD *hn = new TNtupleD("hn","ntuple","Pmod:Charge:TrackTheta:TrackPhi:TrackX:TrackY:MinX:MinY:ChargeMIP:ThetaCerenkov:NPhotons:MipIndex:Chamber:Particle"); + printf("\n\n"); + printf("Pattern Recognition done for event %5i",0); + AliMagF * magf = gAlice->Field(); + AliTracker::SetFieldMap(magf); + for(Int_t iEvtN=0;iEvtNGetRunLoader()->GetNumberOfEvents();iEvtN++) { + GetLoader()->GetRunLoader()->GetEvent(iEvtN); + AliRICHTracker *tr = new AliRICHTracker(); + tr->RecWithStack(hn); +// Info("CheckPR","Pattern Recognition done for event %i \b",iEvtN); + printf("\b\b\b\b\b%5i",iEvtN+1); + } + printf("\n\n"); + pFile->Write();pFile->Close(); +} diff --git a/RICH/AliRICH.h b/RICH/AliRICH.h index bd42271e423..3a6bb71b480 100644 --- a/RICH/AliRICH.h +++ b/RICH/AliRICH.h @@ -1,10 +1,9 @@ #ifndef AliRICH_h #define AliRICH_h - /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * See cxx source for full Copyright notice */ -#include +//#include #include #include #include @@ -178,17 +177,18 @@ public: virtual Int_t IsVersion() const =0; //interface from virtual void StepManager() =0; //interface from AliMC virtual void Hits2SDigits(); //interface from AliSimulation + virtual void Digits2Raw(); //interface from AliSimulation virtual AliDigitizer* CreateDigitizer(AliRunDigitizer* man) const {return new AliRICHDigitizer(man);} //interface from AliSimulation 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 - void GeomPadPanelFrame(); //defines PPF geometry - void GeomAmpGap(); //defines gap geometry + anod wires - void GeomRadiators(); //defines radiators geometry - void GeomSandBox(); //defines sandbox geometry - void GeomRadioSrc(); //defines radio source geometry - void GeomAerogel(); //defines aerogel geometry + void GeomPadPanelFrame()const; //defines PPF geometry + void GeomAmpGap() const; //defines gap geometry + anod wires + void GeomRadiators() const; //defines radiators geometry + void GeomSandBox() const; //defines sandbox geometry + void GeomRadioSrc() const; //defines radio source geometry + void GeomAerogel() const; //defines aerogel geometry virtual void BuildGeometry(); //interface //private part static Float_t Fresnel(Float_t ene,Float_t pdoti, Bool_t pola); //deals with Fresnel absorption @@ -196,7 +196,6 @@ public: 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 -// 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;}} @@ -214,18 +213,19 @@ public: void PrintDigits (Int_t iEvent=0); //prints a list of digits for a given event void PrintClusters(Int_t iEvent=0); //prints a list of clusters for a given event void PrintTracks (Int_t iEvent=0); //prints a list of tracks for a given event - static Int_t Nparticles (Int_t iPID,Int_t iEvent=0,AliRunLoader *p=0); //returns a number of particles - AliRICHhit* Hit(Int_t tid); //returns pointer of the first RICH hit created by a given particle - + void CheckPR ()const; //utility-> run staff for stack ?????? + AliRICHhit* Hit(Int_t tid)const; //returns pointer of the first RICH hit created by a given particle inline void AddHit(Int_t chamber,Int_t tid,TVector3 in3,TVector3 out3,Double_t eloss=0); //add new hit inline void AddSDigit(Int_t c,TVector pad,Double_t q,Int_t pid,Int_t tid); //add new sdigit inline void AddDigit(int c,TVector pad,int q,int cfm,int *tid); //add new digit // 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);} inline void AddCluster(AliRICHcluster &cl); //add new cluster + using AliModule::AddHit; //to get rid of virtual hidden warning + using AliModule::AddDigit; //to get rid of virtual hidden warning + using AliModule::Digits; //to get rid of virtual hidden warning protected: - enum EMedia {kAir=1,kRoha,kSiO2,kC6F14,kCH4,kCsI,kGridCu,kOpSiO2,kGap,kAl,kGlass,kCu,kW - ,kSteel,kPerpex,kSr90,kMylar,kGel,kReflector}; + enum EMedia {kAir=1,kRoha,kSiO2,kC6F14,kCH4,kCsI,kGridCu,kOpSiO2,kGap,kAl,kGlass,kCu,kW,kSteel,kPerpex,kSr90,kMylar,kGel,kReflector}; enum ECounters {kStepManager=0,kCerProdTot,kCerProdRad,kCerKillTot,kCerKillRad,kCerKillRef,kEleProdTot}; AliRICHParam *fpParam; //main RICH parametrization //fHits and fDigits belong to AliDetector diff --git a/RICH/AliRICHChamber.cxx b/RICH/AliRICHChamber.cxx index 4221e04c9d0..adba02d2d18 100644 --- a/RICH/AliRICHChamber.cxx +++ b/RICH/AliRICHChamber.cxx @@ -32,10 +32,11 @@ AliRICHChamber::AliRICHChamber(Int_t iChamber):TNamed() const Double_t kAngVer=20; // vertical angle between chambers 20 grad const Double_t kAngCom=30; // common RICH rotation with respect to x axis 20 grad - RotY(90);//rotate around y since initila position is in XY plane - fCenterX3.SetXYZ(490,0,0); //shift center along x by 490 cm - fPcX3.SetXYZ(490+8-0.2,0,0);//position of the center of apmlification gap (anod wires plane) - fRadX3.SetXYZ(490-2,0,0); //position of the entrance to freon + RotY(90);//rotate around y since initial position is in XY plane + fRad .SetXYZ(490-2 ,0,0); //position of the entrance to freon 1.5 cm of freon+0.5 cm of quartz window + fCenter.SetXYZ(490 ,0,0); //shift center along x by 490 cm + fAnod .SetXYZ(490+8-0.2,0,0); //position of the center of apmlification gap (anod wires plane) + fPc .SetXYZ(490+8 ,0,0); //position of the center of PC switch(iChamber){ case 1: RotY(kAngHor); RotZ(-kAngVer); //right and down @@ -72,6 +73,9 @@ void AliRICHChamber::Print(Option_t *opt) const { // Debug printout TNamed::Print(opt); - fCenterX3.Print(opt); + fRad.Print(opt); + fCenter.Print(opt); + fAnod.Print(opt); + fPc.Print(opt); }//Print() //__________________________________________________________________________________________________ diff --git a/RICH/AliRICHChamber.h b/RICH/AliRICHChamber.h index 2f9979cd2c9..006bb8bd98c 100644 --- a/RICH/AliRICHChamber.h +++ b/RICH/AliRICHChamber.h @@ -21,42 +21,47 @@ public: virtual ~AliRICHChamber() {;} AliRICHChamber& operator=(const AliRICHChamber&) {return *this;} - static Double_t AlphaFeedback(Int_t ) {return 0.020;} //determines number of feedback photons updated to 9/11/04 by Di Mauro + static Double_t AlphaFeedback(Int_t ) {return 0.02;} //determines number of feedback photons updated to 9/11/04 by Di Mauro TRotMatrix* RotMatrix() const{return fpRotMatrix;} TString RotMatrixName() const{return "rot"+fName;} TRotation Rot() const{return fRot;} - Double_t Rho() const{return fCenterX3.Mag();} //gives distance to chamber center in MRS - Double_t ThetaD() const{return fCenterX3.Theta()*TMath::RadToDeg();} //gives polar angle of chamber center in MRS - Double_t PhiD() const{return fCenterX3.Phi() *TMath::RadToDeg();} //gives azimuthal angle of chamber center in MRS - Double_t ThetaXd() const{return fRot.ThetaX() *TMath::RadToDeg();} - Double_t PhiXd() const{return fRot.PhiX() *TMath::RadToDeg();} - Double_t ThetaYd() const{return fRot.ThetaY() *TMath::RadToDeg();} - Double_t PhiYd() const{return fRot.PhiY() *TMath::RadToDeg();} - Double_t ThetaZd() const{return fRot.ThetaZ() *TMath::RadToDeg();} - Double_t PhiZd() const{return fRot.PhiZ() *TMath::RadToDeg();} - void RotX(Double_t a) {a*=TMath::DegToRad();fRot.RotateX(a);fCenterX3.RotateX(a);fRadX3.RotateX(a);fPcX3.RotateX(a);}//degrees around X - void RotY(Double_t a) {a*=TMath::DegToRad();fRot.RotateY(a);fCenterX3.RotateY(a);fRadX3.RotateY(a);fPcX3.RotateY(a);}//degrees around Y - void RotZ(Double_t a) {a*=TMath::DegToRad();fRot.RotateZ(a);fCenterX3.RotateZ(a);fRadX3.RotateZ(a);fPcX3.RotateZ(a);}//degrees around Z - TVector3 Rad() const{return fRadX3;} //provides center of radiator position in MRS, cm - TVector3 Pc() const{return fPcX3;} //provides center of photocathond position in MRS, cm - TVector3 Center() const{return fCenterX3;} //provides center of chamber position in MRS, cm - void Print(Option_t *sOption)const; //virtual interface from TObject + Double_t Rho() const{return fCenter.Mag();} //gives distance to chamber center in MRS + Double_t ThetaD() const{return fCenter.Theta()*TMath::RadToDeg();} //gives polar angle of chamber center in MRS + Double_t PhiD() const{return fCenter.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 RotX(Double_t a) {a*=TMath::DegToRad();fRot.RotateX(a);fRad.RotateX(a);fCenter.RotateX(a);fAnod.RotateX(a);fPc.RotateX(a);}//degrees around X + void RotY(Double_t a) {a*=TMath::DegToRad();fRot.RotateY(a);fRad.RotateY(a);fCenter.RotateY(a);fAnod.RotateY(a);fPc.RotateY(a);}//degrees around Y + void RotZ(Double_t a) {a*=TMath::DegToRad();fRot.RotateZ(a);fRad.RotateZ(a);fCenter.RotateZ(a);fAnod.RotateZ(a);fPc.RotateZ(a);}//degrees around Z + TVector3 Rad() const{return fRad;} //provides center of radiator position in MRS, cm + TVector3 Anod() const{return fAnod;} //provides center of anod wires plane in MRS, cm + TVector3 Pc() const{return fPc;} //provides center of photocathode position in MRS, cm + TVector3 Center() const{return fCenter;} //provides center of chamber position (exit from quartz window) in MRS, cm + void Print(Option_t *sOption)const; + TVector3 PMrs2Loc(TVector3 p3)const{TVector3 ploc=Rot().Invert()*p3;ploc.SetXYZ(-ploc.Px(),ploc.Py(),ploc.Pz()); return ploc;} //momentum MARS-local +//Transformations for radiator plane + TVector2 Mrs2Rad(TVector3 x3)const{x3-=fRad;x3.Transform(fRot.Inverse());return TVector2(-x3.X()+0.5*AliRICHParam::PcSizeX(),x3.Y()+0.5*AliRICHParam::PcSizeY());} + TVector3 Rad2Mrs(TVector2 x2)const{TVector3 x3(-x2.X()+0.5*AliRICHParam::PcSizeX(),x2.Y()-0.5*AliRICHParam::PcSizeY(),0);x3.Transform(fRot); x3+=fRad;return x3;} +//Transformations for anod wires plane + TVector2 Mrs2Anod(TVector3 x3)const{x3-=fAnod;x3.Transform(fRot.Inverse());return TVector2(-x3.X()+0.5*AliRICHParam::PcSizeX(),x3.Y()+0.5*AliRICHParam::PcSizeY());} + TVector3 Anod2Mrs(TVector2 x2)const{TVector3 x3(-x2.X()+0.5*AliRICHParam::PcSizeX(),x2.Y()-0.5*AliRICHParam::PcSizeY(),0);x3.Transform(fRot); x3+=fAnod;return x3;} //Transformations for photcathode plane - TVector2 Mrs2Pc(TVector3 x3)const{x3-=fPcX3;x3.Transform(fRot.Inverse());return TVector2(-x3.X()+0.5*AliRICHParam::PcSizeX(),x3.Y()+0.5*AliRICHParam::PcSizeY());} - TVector3 Pc2Mrs(TVector2 x2)const{TVector3 x3(-x2.X()+0.5*AliRICHParam::PcSizeX(),x2.Y()-0.5*AliRICHParam::PcSizeY(),0);x3.Transform(fRot); x3+=fPcX3;return x3;} + TVector2 Mrs2Pc(TVector3 x3)const{x3-=fPc;x3.Transform(fRot.Inverse());return TVector2(-x3.X()+0.5*AliRICHParam::PcSizeX(),x3.Y()+0.5*AliRICHParam::PcSizeY());} + TVector3 Pc2Mrs(TVector2 x2)const{TVector3 x3(-x2.X()+0.5*AliRICHParam::PcSizeX(),x2.Y()-0.5*AliRICHParam::PcSizeY(),0);x3.Transform(fRot); x3+=fPc;return x3;} TVector2 Mrs2Pc(TLorentzVector x4) const{return Mrs2Pc(x4.Vect());} -//Transformations for radiator plane - TVector2 Mrs2Rad(TVector3 x3)const{x3-=fRadX3;x3.Transform(fRot.Inverse());return TVector2(-x3.X()+0.5*AliRICHParam::PcSizeX(),x3.Y()+0.5*AliRICHParam::PcSizeY());} - TVector3 Rad2Mrs(TVector2 x2)const{TVector3 x3(-x2.X()+0.5*AliRICHParam::PcSizeX(),x2.Y()-0.5*AliRICHParam::PcSizeY(),0);x3.Transform(fRot); x3+=fRadX3;return x3;} - TVector3 PMrs2Loc(TVector3 p3)const{TVector3 ploc=Rot().Invert()*p3;ploc.SetXYZ(-ploc.Px(),ploc.Py(),ploc.Pz()); return ploc;} protected: - TVector3 fCenterX3; //chamber center position in MRS (cm) - TVector3 fRadX3; //radiator entrance center position in MRS (cm) - TVector3 fPcX3; //PC center position in MRS (cm) + TVector3 fRad; //radiator entrance center position in MRS (cm) + TVector3 fCenter; //chamber center position (quartz window exit) in MRS (cm) + TVector3 fAnod; //anod wires plane center position in MRS (cm) + TVector3 fPc; //PC center position in MRS (cm) TRotation fRot; //chamber rotation in MRS TRotMatrix *fpRotMatrix; //rotation matrix of the chamber with respect to MRS - ClassDef(AliRICHChamber,7) //single RICH chamber description + ClassDef(AliRICHChamber,8) //single RICH chamber description };//class AliRICHChamber #endif //AliRICHChamber_h diff --git a/RICH/AliRICHClusterFinder.cxx b/RICH/AliRICHClusterFinder.cxx index 536462f02a1..5b70a4a2a26 100644 --- a/RICH/AliRICHClusterFinder.cxx +++ b/RICH/AliRICHClusterFinder.cxx @@ -39,14 +39,14 @@ AliRICHClusterFinder::AliRICHClusterFinder(AliRICH *pRICH) AliDebug(1,"main ctor Stop."); }//main ctor //__________________________________________________________________________________________________ -void AliRICHClusterFinder::Exec() +void AliRICHClusterFinder::Exec(const Option_t *) { //Main method of cluster finder. Loops on events and chambers, everything else is done in FindClusters() AliDebug(1,"Exec Start."); R()->GetLoader() ->LoadDigits(); // R()->GetLoader()->GetRunLoader()->LoadHeader(); - R()->GetLoader()->GetRunLoader()->LoadKinematics(); //header is already loaded + if(!R()->GetLoader()->GetRunLoader()->TreeK()) R()->GetLoader()->GetRunLoader()->LoadKinematics(); for(Int_t iEventN=0;iEventNGetEventsPerRun();iEventN++){//events loop AliDebug(1,Form("Processing event %i...",iEventN)); diff --git a/RICH/AliRICHClusterFinder.h b/RICH/AliRICHClusterFinder.h index e9a5598f554..5c03e195e6f 100644 --- a/RICH/AliRICHClusterFinder.h +++ b/RICH/AliRICHClusterFinder.h @@ -16,7 +16,7 @@ public: virtual ~AliRICHClusterFinder() {;} AliRICH *R() {return fRICH;} //returns pointer to RICH - void Exec(); //loop on events and chambers + void Exec(const Option_t *option=""); //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 diff --git a/RICH/AliRICHDisplFast.cxx b/RICH/AliRICHDisplFast.cxx index 7137e0c6cd4..4161408ee54 100644 --- a/RICH/AliRICHDisplFast.cxx +++ b/RICH/AliRICHDisplFast.cxx @@ -17,12 +17,13 @@ #include "AliRICHChamber.h" #include "AliRICHParam.h" #include +#include +#include #include -#include -#include #include -#include #include +#include +#include ClassImp(AliRICHDisplFast) //__________________________________________________________________________________________________ @@ -36,7 +37,6 @@ void AliRICHDisplFast::ShowEvent(Int_t iEvtNmin,Int_t iEvtNmax) if(!isDigits){Error("ShoEvent","No digits. Nothing to display.");return;} TCanvas *canvas = new TCanvas("RICHDisplay","RICH Display",0,0,1226,900); -// TCanvas *canvas = (TCanvas*)gROOT->FindObjectAny("RICHDisplay"); gStyle->SetPalette(1); @@ -54,8 +54,8 @@ void AliRICHDisplFast::ShowEvent(Int_t iEvtNmin,Int_t iEvtNmax) if(iEvtNmax>gAlice->GetEventsPerRun()) iEvtNmax=gAlice->GetEventsPerRun(); - TText *t = new TText(); - t->SetTextSize(0.10); + TLatex t; + t.SetTextSize(0.10); TText *tit = new TText(0.1,0.6,"RICH Display"); tit->SetTextSize(0.10); @@ -64,8 +64,8 @@ void AliRICHDisplFast::ShowEvent(Int_t iEvtNmin,Int_t iEvtNmax) canvas->cd(1); sprintf(titdisp,"Event Number %i",iEventN); - t->SetText(0.2,0.4,titdisp); - t->Draw(); + t.SetText(0.2,0.4,titdisp); + t.Draw(); tit->Draw(); pRich->GetLoader()->GetRunLoader()->GetEvent(iEventN); pRich->GetLoader()->TreeD()->GetEntry(0); @@ -75,7 +75,7 @@ void AliRICHDisplFast::ShowEvent(Int_t iEvtNmin,Int_t iEvtNmax) AliRICHdigit *pDig = (AliRICHdigit*)pRich->Digits(iChamber)->At(j); TVector2 x2=AliRICHParam::Pad2Loc(pDig->Pad()); pDigitsH2[iChamber]->Fill(x2.X(),x2.Y()); - } + }//digits loop if(iChamber==1) canvas->cd(7); if(iChamber==2) canvas->cd(8); if(iChamber==3) canvas->cd(4); @@ -89,10 +89,8 @@ void AliRICHDisplFast::ShowEvent(Int_t iEvtNmin,Int_t iEvtNmax) canvas->Modified(); if(iEvtNminWaitPrimitive(); -// canvas->cd(3); -// gPad->Clear(); } -} +}//ShowEvent() //__________________________________________________________________________________________________ void AliRICHDisplFast::Exec(Option_t *) { @@ -151,7 +149,7 @@ void AliRICHDisplFast::Exec(Option_t *) pHitsH2->SetTitle(Form("event %i chamber %2i",iEventN,iChamber)); pHitsH2->SetMarkerColor(kRed); pHitsH2->SetMarkerStyle(29); pHitsH2->SetMarkerSize(0.4); pHitsH2->Draw(); - DrawSectors(); + AliRICHParam::DrawSectors(); TLatex l; l.SetNDC(); l.SetTextSize(0.02); if(!isHits) {l.SetTextColor(kRed) ;l.DrawLatex(0.1,0.01,"No Hits" );} if(!isDigits) {l.SetTextColor(kGreen);l.DrawLatex(0.4,0.01,"No DIGITS" );} @@ -194,21 +192,20 @@ void AliRICHDisplFast::Exec(Option_t *) if(isClusters) pRich->GetLoader()->UnloadRecPoints(); }//Exec() //__________________________________________________________________________________________________ -void AliRICHDisplFast::DrawSectors() -{ - 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()}; +Int_t AliRICHDisplFast::Nparticles(Int_t iPartID,Int_t iEvtN,AliRunLoader *pRL) +{ +//counts total number of particles of given type (including secondary) for a given event + pRL->GetEvent(iEvtN); + if(pRL->LoadHeader()) return 0; + if(pRL->LoadKinematics()) return 0; + AliStack *pStack=pRL->Stack(); - 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()}; + Int_t iCounter=0; + for(Int_t i=0;iGetNtrack();i++){ + if(pStack->Particle(i)->GetPdgCode()==iPartID) iCounter++; + } - 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() + pRL->UnloadHeader(); + pRL->UnloadKinematics(); + return iCounter; +} diff --git a/RICH/AliRICHDisplFast.h b/RICH/AliRICHDisplFast.h index 78feb0014b7..2f270cf8cf2 100644 --- a/RICH/AliRICHDisplFast.h +++ b/RICH/AliRICHDisplFast.h @@ -5,23 +5,21 @@ * See cxx source for full Copyright notice */ #include -#include "AliRICH.h" -#include -#include "TH2.h" -class AliRICH; + +class AliRunLoader; class AliRICHDisplFast : public TTask { public : AliRICHDisplFast() {;} virtual ~AliRICHDisplFast() {;} - static void DrawSectors(); //Draw sectors in plot - void ShowEvent(Int_t iEvtNmin,Int_t iEvtNmax); // Display looping on events - void ShowEvent(Int_t iEvent) {ShowEvent(iEvent,iEvent);} // Display only one event - virtual void Exec(Option_t *opt=0); //virtual do the main job + void ShowEvent(Int_t iEvtNmin,Int_t iEvtNmax); //Display looping on events + void ShowEvent(Int_t iEvent) {ShowEvent(iEvent,iEvent);} //Display only one event + static Int_t Nparticles (Int_t iPID,Int_t iEvent=0,AliRunLoader *p=0); //returns a number of particles of given type + virtual void Exec(Option_t *opt=0); //virtual do the main job protected: ClassDef(AliRICHDisplFast,0) //Utility class to draw the current event topology }; -#endif // #ifdef AliRICHDisplFast_cxx +#endif //AliRICHDisplFast_cxx diff --git a/RICH/AliRICHHelix.cxx b/RICH/AliRICHHelix.cxx index c5ad6d6501a..92b67f626d7 100644 --- a/RICH/AliRICHHelix.cxx +++ b/RICH/AliRICHHelix.cxx @@ -3,10 +3,14 @@ ClassImp(AliRICHHelix) -void AliRICHHelix::Print(Option_t *opt) const +void AliRICHHelix::Print(Option_t *) const { -// Debug printout - fX0.Print(opt);fP0.Print(opt); - AliInfo("Point of interest:"); - fX.Print(opt); fP.Print(opt); -}//Print() +//Debug printout + AliInfo(Form("Q=%i, in B(0,0,%5.2f) x0=(%5.2f,%5.2f,%5.2f) p0=(%5.2f,%5.2f,%5.2f)",fQ,fBz, + fX0.X(),fX0.Y(),fX0.Z(), + fP0.X(),fP0.Y(),fP0.Z() )); + AliInfo(Form("At length %7.2f gives x=(%5.2f,%5.2f,%5.2f) p=(%5.2f,%5.2f,%5.2f)",fLen, + fX.X(),fX.Y(),fX.Z(), + fP.X(),fP.Y(),fP.Z() )); + +} diff --git a/RICH/AliRICHHelix.h b/RICH/AliRICHHelix.h index 09737df2528..18f30bc34d4 100644 --- a/RICH/AliRICHHelix.h +++ b/RICH/AliRICHHelix.h @@ -9,37 +9,36 @@ class AliRICHHelix: public TObject { public: - AliRICHHelix() {;} - AliRICHHelix(TVector3 x0,TVector3 p0,Int_t q,Double_t b=0.4) {fX0=x0;fP0=p0;fQ=q;fBz=b;} //takes position and momentum at parametrised point - virtual ~AliRICHHelix() {;} - inline void Propagate(Double_t lenght); - Bool_t Intersection(TVector3 plane) {return Intersection(plane,plane.Unit());} // special plane given by point only - inline Bool_t Intersection(TVector3 planePoint,TVector3 planeNorm); //intersection with plane given by point and normal vector - inline Int_t RichIntersect(AliRICHParam *pParam); - TVector3 X()const {return fX;} - TVector3 P()const {return fP;} - TVector3 X0()const {return fX0;} - TVector3 P0()const {return fP0;} - TVector2 PosPc() const {return fPosPc;} //returns position of intersection with PC - TVector2 PosRad()const {return fPosRad;} //returns position of intersection with radiator - TVector3 Ploc()const {return fPloc;} //returns momentum at the position of intersection with radiator - - void Print(Option_t *sOption)const; //virtual interface from TObject + AliRICHHelix():TObject() {;} + AliRICHHelix(const TVector3 &x0,const TVector3 &p0,Int_t q=1,Double_t b=0.2):TObject(),fX0(x0),fP0(p0),fX(x0),fP(p0), + fLen(0),fQ(q),fBz(b) {;} + virtual ~AliRICHHelix() {;} + + inline void Propagate(Double_t len); //propogate helix by length len along it + inline Int_t RichIntersect(AliRICHParam *pParam); //search intersection with any RICH chamber + inline Bool_t Intersection(TVector3 planePnt,TVector3 planeNorm); //intersection with plane given by point and normal vector + Bool_t Intersection(TVector3 pl) {return Intersection(pl,pl.Unit());} // special plane given by point only + TVector2 PosRad() const{return fPosRad;} //returns position of intersection with radiator (local system) + TVector2 PosAnod() const{return fPosAnod;} //returns position of intersection with anod wires plane (local system) + TVector2 PosPc() const{return fPosPc;} //returns position of intersection with PC (local system) + TVector3 Ploc() const{return fPloc;} //returns momentum at the position of intersection with radiator + void Print(Option_t *sOption)const; //virtual interface from TObject protected: TVector3 fX0; //helix position in parametrised point, cm in MRS TVector3 fP0; //helix momentum in parametrised point, GeV/c in MRS TVector3 fX; //helix position in point of interest, cm in MRS TVector3 fP; //helix momentum in point of interest, GeV/c in MRS - Double_t fLength; //helix length in point of interest + Double_t fLen; //helix length in point of interest Int_t fQ; //sign of track charge (value not provided by current ESD) Double_t fBz; //magnetic field along z value in Tesla under assumption of uniformity - TVector2 fPosPc; //position on PC in local system - TVector2 fPosRad; //position on radiator in local system + TVector2 fPosRad; //track intersection with radiator (local system) + TVector2 fPosAnod; //track intersection with anod wires plane (local system) + TVector2 fPosPc; //track intersection with PC (local system) TVector3 fPloc; //momentum in local system ClassDef(AliRICHHelix,0) //General helix };//class AliRICHHelix //__________________________________________________________________________________________________ -void AliRICHHelix::Propagate(Double_t length) +void AliRICHHelix::Propagate(Double_t len) { // Propogates the helix to the position of interest defined by helix length s // Assumes uniform magnetic field along z direction. @@ -47,13 +46,13 @@ void AliRICHHelix::Propagate(Double_t length) Double_t a = -c*fBz*fQ; Double_t rho = a/fP0.Mag(); - fX.SetX( fX0.X()+fP0.X()*TMath::Sin(rho*length)/a-fP0.Y()*(1-TMath::Cos(rho*length))/a ); - fX.SetY( fX0.Y()+fP0.Y()*TMath::Sin(rho*length)/a+fP0.X()*(1-TMath::Cos(rho*length))/a ); - fX.SetZ( fX0.Z()+fP0.Z()*length/fP0.Mag() ); - fP.SetX( fP0.X()*TMath::Cos(rho*length)-fP0.Y()*TMath::Sin(rho*length) ); - fP.SetY( fP0.Y()*TMath::Cos(rho*length)+fP0.X()*TMath::Sin(rho*length) ); - fP.SetZ( fP0.Z() ); - fLength=length; + fX.SetX( fX0.X()+fP0.X()*TMath::Sin(rho*len)/a-fP0.Y()*(1-TMath::Cos(rho*len))/a ); + fX.SetY( fX0.Y()+fP0.Y()*TMath::Sin(rho*len)/a+fP0.X()*(1-TMath::Cos(rho*len))/a ); + fX.SetZ( fX0.Z()+fP0.Z()*len/fP0.Mag() ); + fP.SetX( fP0.X()*TMath::Cos(rho*len)-fP0.Y()*TMath::Sin(rho*len) ); + fP.SetY( fP0.Y()*TMath::Cos(rho*len)+fP0.X()*TMath::Sin(rho*len) ); + fP.SetZ( fP0.Z() ); + fLen=len; } //__________________________________________________________________________________________________ Bool_t AliRICHHelix::Intersection(TVector3 planePoint,TVector3 planeNorm) @@ -62,9 +61,9 @@ Bool_t AliRICHHelix::Intersection(TVector3 planePoint,TVector3 planeNorm) // Returns kTrue if helix intersects the plane, kFALSE otherwise. // Stores result in current helix fields fX and fP. - Double_t s=(planePoint-fX0)*planeNorm,dist=99999,distPrev=dist; + Double_t s=(planePoint-fX0)*planeNorm,dist=99999,distPrev=dist;//estimates initial distance to plane - while(TMath::Abs(dist)>0.0001){ + while(TMath::Abs(dist)>0.00001){ Propagate(s); //calculates helix at the distance s from x0 ALONG the helix dist=(fX-planePoint)*planeNorm; //distance between current helix position and plane if(TMath::Abs(dist) >= TMath::Abs(distPrev)) { return kFALSE;} @@ -83,13 +82,19 @@ Int_t AliRICHHelix::RichIntersect(AliRICHParam *pParam) if(Intersection(pParam->C(iChamberN)->Rad())){//there is intersection with radiator plane fPosRad=pParam->C(iChamberN)->Mrs2Rad(fX);//position on radiator plane if(pParam->IsAccepted(fPosRad)){//intersection within radiator (even if in dead zone) + if(Intersection(pParam->C(iChamberN)->Pc())){//there is intersection with photocathode fPosPc=pParam->C(iChamberN)->Mrs2Pc(fX);//position on photcathode plane if(pParam->IsAccepted(fPosPc)){//intersection within pc (even if in dead zone) - fPloc=pParam->C(iChamberN)->PMrs2Loc(fP); + + Intersection(pParam->C(iChamberN)->Anod()); //search for anod intersection position + fPosAnod=pParam->C(iChamberN)->Mrs2Anod(fX); + + fPloc=pParam->C(iChamberN)->PMrs2Loc(fP);//trasform p to local system return iChamberN; }//if inside PC - }//if for PC + }//if intersects PC + }//if inside radiator }//if for radiator }//chamber loop diff --git a/RICH/AliRICHMap.h b/RICH/AliRICHMap.h index cca94dd8164..d57edb4f7b2 100644 --- a/RICH/AliRICHMap.h +++ b/RICH/AliRICHMap.h @@ -1,7 +1,5 @@ #ifndef AliRICHMap_h #define AliRICHMap_h - - /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * See cxx source for full Copyright notice */ @@ -24,7 +22,7 @@ public: void FlagHit(Int_t ix,Int_t iy) {(*fMap)(ix, iy)=-TMath::Abs((*fMap)(ix,iy));} //virtual Bool_t ValidateHit(Int_t,Int_t) {return 1;} //virtual inline FlagType TestHit(Int_t ix,Int_t iy); //virtual - void Print() const{fMap->Print();} + void Print(const Option_t *) const{fMap->Print();} protected: TClonesArray *fDigits; //List of digits Int_t fNdigits; //Number of digits diff --git a/RICH/AliRICHParam.cxx b/RICH/AliRICHParam.cxx index c13575596b4..51ba87e41a3 100644 --- a/RICH/AliRICHParam.cxx +++ b/RICH/AliRICHParam.cxx @@ -22,10 +22,12 @@ #include #include #include +#include ClassImp(AliRICHParam) Bool_t AliRICHParam::fgIsWireSag =kTRUE; //take ware sagita into account? Bool_t AliRICHParam::fgIsResolveClusters =kTRUE; //do cluster resolving? +Bool_t AliRICHParam::fgIsFeedback =kTRUE; //generate feedback photons? Bool_t AliRICHParam::fgIsRadioSrc =kFALSE; //put radioactive source instead of radiators? Bool_t AliRICHParam::fgIsAerogel =kFALSE; //special aerogel configuration Bool_t AliRICHParam::fgIsTestBeam =kFALSE; //special test beam configuration @@ -38,7 +40,7 @@ Float_t AliRICHParam::fgSigmaThSpread =0.035; // //__________________________________________________________________________________________________ void AliRICHParam::Print(Option_t*) const { - AliInfo(Form("Pads in chamber (%3i,%3i) in sector (%2i,%2i)",NpadsX(),NpadsY(),NpadsXsec(),NpadsYsec())); + AliInfo(Form("Pads in chamber (%3i,%3i) in sector (%2i,%2i) pad size (%4.2f,%4.2f)",NpadsX(),NpadsY(),NpadsXsec(),NpadsYsec(),PadSizeX(),PadSizeY())); AliInfo(Form("Resolve clusters %i sagita %i Radio source %i Aerogel %i TestBeam %i", IsResolveClusters(),IsWireSag(),IsRadioSrc(),IsAerogel(),IsTestBeam())); fpChambers->Print(); @@ -82,100 +84,59 @@ Float_t AliRICHParam::AbsCH4(Float_t eV) }//AbsoCH4() //__________________________________________________________________________________________________ void AliRICHParam::TestSeg() -{ - - new TCanvas("name","PC segmentation"); - gPad->Range(-20,-20,200,150); - AliRICHDisplFast::DrawSectors(); - - TLatex t; t.SetTextSize(0.02); - t.DrawText(0,140,"View from interaction point"); - t.DrawLatex(PcSizeX()+10,120,Form("Pc %6.2fx%6.2fcm %3ix%3ipads",PcSizeX() ,PcSizeY(), NpadsX() ,NpadsY())); - t.DrawLatex(PcSizeX()+10,115,Form("Sec %6.2fx%5.2fcm %3ix%2ipads",SectorSizeX(),SectorSizeY(),NpadsXsec(),NpadsYsec())); - t.DrawLatex(PcSizeX()+10,110,Form("Pad %6.2fx%4.2fcm DeadZone %6.2fcm",PadSizeX() ,PadSizeY(),DeadZone())); - - TVector2 v2; - t.SetTextAlign(12); - v2=AliRICHParam::Pad2Loc( 1,24); t.DrawText(v2.X(),v2.Y(),"sec 1"); - v2=AliRICHParam::Pad2Loc(81,24); t.DrawText(v2.X(),v2.Y(),"sec 2"); - v2=AliRICHParam::Pad2Loc( 1,70); t.DrawText(v2.X(),v2.Y(),"sec 3"); - v2=AliRICHParam::Pad2Loc(81,70); t.DrawText(v2.X(),v2.Y(),"sec 4"); - v2=AliRICHParam::Pad2Loc( 1,120); t.DrawText(v2.X(),v2.Y(),"sec 5"); - v2=AliRICHParam::Pad2Loc(81,120); t.DrawText(v2.X(),v2.Y(),"sec 6"); +{ + new TCanvas("pads","PC segmentation - pads display",700,600); + gPad->Range(-5,-5,PcSizeX()+5,PcSizeY()+15); + TVector p(2); TVector2 c; TVector2 b; //current: pad, pad center, pad boundary +// list of corners: + Double_t x0=0,x1=SectorSizeX(),x2=SectorSizeX()+DeadZone(), x3=PcSizeX(); + Double_t y0=0,y1=SectorSizeY(),y2=SectorSizeY()+DeadZone(),y3=2*SectorSizeY()+DeadZone(),y4=PcSizeY()-SectorSizeY(),y5=PcSizeY(); + DrawSectors(); +//header + TLatex t; + t.SetTextSize(0.02); t.SetTextColor(kBlack); t.SetTextAlign(11); + t.DrawLatex(0,PcSizeY()+10,Form("IP in front of this page. pad size %.2fx%.2fcm dead zone %.2fcm",PadSizeX(),PadSizeY(),DeadZone())); + t.DrawLatex(0,PcSizeY()+ 5,Form("Pc %.2fx%.2f cm %ix%i pads Sec %.2fx%.2f cm %ix%i pads", + PcSizeX() , PcSizeY() , NpadsX() , NpadsY() , + SectorSizeX() , SectorSizeY() , NpadsXsec() , NpadsYsec() )); +//sectors + t.SetTextSize(0.015); t.SetTextColor(kRed); t.SetTextAlign(22); + c=Pad2Loc( 40, 24); t.DrawText(c.X(),c.Y(),Form("sec 1 (%.2f,%.2f)",c.X(),c.Y() )); + c=Pad2Loc( 40, 75); t.DrawText(c.X(),c.Y(),Form("sec 3 (%.2f,%.2f)",c.X(),c.Y() )); + c=Pad2Loc( 40,121); t.DrawText(c.X(),c.Y(),Form("sec 5 (%.2f,%.2f)",c.X(),c.Y() )); + c=Pad2Loc(120, 24); t.DrawText(c.X(),c.Y(),Form("sec 2 (%.2f,%.2f)",c.X(),c.Y() )); + c=Pad2Loc(120, 75); t.DrawText(c.X(),c.Y(),Form("sec 4 (%.2f,%.2f)",c.X(),c.Y() )); + c=Pad2Loc(120,121); t.DrawText(c.X(),c.Y(),Form("sec 6 (%.2f,%.2f)",c.X(),c.Y() )); +//coners + t.SetTextSize(0.015); t.SetTextColor(kBlue); + + b.Set(x0,y0);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(11);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x0,y1);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(13);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x0,y2);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(11);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x0,y3);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(13);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x0,y4);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(11);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x0,y5);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(13);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); -// TGaxis *pAx=new TGaxis(0,0,140, 0,0,140,510,"-="); pAx->SetTitle("x, cm"); pAx->SetTextSize(0.05); pAx->Draw(); -// TGaxis *pAy=new TGaxis(0,0, 0,140,0,140,510,"-="); pAy->SetTitle("y, cm"); pAy->SetTextSize(0.05); pAy->Draw(); - + b.Set(x1,y0);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(31);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x1,y1);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(33);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x1,y2);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(31);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x1,y3);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(33);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x1,y4);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(31);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x1,y5);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(33);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); - t.SetTextColor(kBlue); - - Double_t margin=5; - 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(); - -//write pad numbers along x - t.SetTextAlign(13); t.DrawText(x0,y0,"1"); - t.SetTextAlign(33); t.DrawText(x1,y0,"80"); - t.SetTextAlign(13); t.DrawText(x2,y0,"81"); - t.SetTextAlign(33); t.DrawText(x3,y0,"160"); -//write pad numbers along y - t.SetTextAlign(31); t.DrawText(x0,y0,"1"); - t.SetTextAlign(33); t.DrawText(x0,y1,"48"); - t.SetTextAlign(31); t.DrawText(x0,y2,"49"); - t.SetTextAlign(33); t.DrawText(x0,y3,"96"); - t.SetTextAlign(31); t.DrawText(x0,y4,"97"); - t.SetTextAlign(33); t.DrawText(x0,y5,"144"); + b.Set(x2,y0);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(11);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x2,y1);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(13);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x2,y2);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(11);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x2,y3);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(13);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x2,y4);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(11);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x2,y5);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(13);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); - -//positions along x - t.SetTextColor(kGreen); - t.SetTextAlign(11);t.DrawText(x0,y0-margin,Form("%5.2f",x0)); - t.SetTextAlign(31);t.DrawText(x1,y0-margin,Form("%5.2f",x1)); - t.SetTextAlign(11);t.DrawText(x2,y0-margin,Form("%5.2f",x2)); - t.SetTextAlign(31);t.DrawText(x3,y0-margin,Form("%5.2f",x3)); -//positions along y - t.SetTextAlign(31);t.DrawText(x0-margin,y0,Form("%5.2f",y0)); - t.SetTextAlign(33);t.DrawText(x0-margin,y1,Form("%5.2f",y1)); - t.SetTextAlign(31);t.DrawText(x0-margin,y2,Form("%5.2f",y2)); - t.SetTextAlign(33);t.DrawText(x0-margin,y3,Form("%5.2f",y3)); - t.SetTextAlign(31);t.DrawText(x0-margin,y4,Form("%5.2f",y4)); - t.SetTextAlign(33);t.DrawText(x0-margin,y5,Form("%5.2f",y5)); -//coners - t.SetTextColor(kRed); - t.SetTextSize(0.01); - TVector pad(2); -//sector 1 - v2=Pad2Loc(Loc2Pad(TVector2(x0,y0)));t.SetTextAlign(11);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - v2=Pad2Loc(Loc2Pad(TVector2(x0,y1)));t.SetTextAlign(13);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - v2=Pad2Loc(Loc2Pad(TVector2(x1,y1)));t.SetTextAlign(33);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - v2=Pad2Loc(Loc2Pad(TVector2(x1,y0)));t.SetTextAlign(31);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); -//secr 3 - v2=Pad2Loc(Loc2Pad(TVector2(x0,y2)));t.SetTextAlign(11);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - v2=Pad2Loc(Loc2Pad(TVector2(x0,y3)));t.SetTextAlign(13);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - v2=Pad2Loc(Loc2Pad(TVector2(x1,y3)));t.SetTextAlign(33);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - v2=Pad2Loc(Loc2Pad(TVector2(x1,y2)));t.SetTextAlign(31);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); -//secr 5 - v2=Pad2Loc(Loc2Pad(TVector2(x0,y4)));t.SetTextAlign(11);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - v2=Pad2Loc(Loc2Pad(TVector2(x0,y5)));t.SetTextAlign(13);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - v2=Pad2Loc(Loc2Pad(TVector2(x1,y5)));t.SetTextAlign(33);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - v2=Pad2Loc(Loc2Pad(TVector2(x1,y4)));t.SetTextAlign(31);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - - v2=Pad2Loc(Loc2Pad(TVector2(x2,y4)));t.SetTextAlign(11);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - v2=Pad2Loc(Loc2Pad(TVector2(x2,y5)));t.SetTextAlign(13);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - v2=Pad2Loc(Loc2Pad(TVector2(x3,y5)));t.SetTextAlign(33);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - v2=Pad2Loc(Loc2Pad(TVector2(x3,y4)));t.SetTextAlign(31);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - - v2=Pad2Loc(Loc2Pad(TVector2(x2,y2)));t.SetTextAlign(11);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - v2=Pad2Loc(Loc2Pad(TVector2(x2,y3)));t.SetTextAlign(13);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - v2=Pad2Loc(Loc2Pad(TVector2(x3,y3)));t.SetTextAlign(33);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - v2=Pad2Loc(Loc2Pad(TVector2(x3,y2)));t.SetTextAlign(31);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - - v2=Pad2Loc(Loc2Pad(TVector2(x2,y0)));t.SetTextAlign(11);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - v2=Pad2Loc(Loc2Pad(TVector2(x2,y1)));t.SetTextAlign(13);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - v2=Pad2Loc(Loc2Pad(TVector2(x3,y1)));t.SetTextAlign(33);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); - v2=Pad2Loc(Loc2Pad(TVector2(x3,y0)));t.SetTextAlign(31);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y())); + b.Set(x3,y0);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(31);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x3,y1);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(33);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x3,y2);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(31);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x3,y3);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(33);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x3,y4);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(31);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); + b.Set(x3,y5);p=Loc2Pad(b);c=Pad2Loc(p);t.SetTextAlign(33);t.DrawText(c.X(),c.Y(),Form("(%.2f,%.2f)-(%.0f,%.0f)-(%.2f,%.2f)",b.X(),b.Y(),p(0),p(1),c.X(),c.Y())); }//TestSeg() //__________________________________________________________________________________________________ void AliRICHParam::TestResp() @@ -250,3 +211,21 @@ void AliRICHParam::DrawAxis() TPolyLine3D *pYaxis=new TPolyLine3D(2,Y);pYaxis->SetLineColor(kGreen); pYaxis->Draw(); TPolyLine3D *pZaxis=new TPolyLine3D(2,Z);pZaxis->SetLineColor(kBlue); pZaxis->Draw(); } +//__________________________________________________________________________________________________ +void AliRICHParam::DrawSectors() +{ + Double_t xLeft[5] = {0,0,SectorSizeX(),SectorSizeX(),0}; + Double_t xRight[5] = {SectorSizeX()+DeadZone(),SectorSizeX()+DeadZone(),PcSizeX(),PcSizeX(),SectorSizeX()+DeadZone()}; + + Double_t yDown[5] = {0,SectorSizeY(),SectorSizeY(),0,0}; + Double_t yCenter[5] = { SectorSizeY()+DeadZone(),2*SectorSizeY()+DeadZone(),2*SectorSizeY()+DeadZone(), + SectorSizeY()+DeadZone(),SectorSizeY()+DeadZone()}; + Double_t yUp[5] = {2*SectorSizeY()+2*DeadZone(),PcSizeY(),PcSizeY(),2*SectorSizeY()+2*DeadZone(),2*SectorSizeY()+2*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/AliRICHParam.h b/RICH/AliRICHParam.h index a8418171824..8053d39d3bc 100644 --- a/RICH/AliRICHParam.h +++ b/RICH/AliRICHParam.h @@ -40,7 +40,8 @@ public: void TestResp(); //test the response group of methodes void TestSeg(); //test the segmentation group of methodes void TestTrans(); //test the transform group of methodes - void DrawAxis(); + static void DrawAxis(); + static void DrawSectors(); //flags staff static void SetAerogel(Bool_t a) {fgIsAerogel=a;} static Bool_t IsAerogel() {return fgIsAerogel;} @@ -105,8 +106,10 @@ public: static Float_t DenGel() {return (RefIdxGel(0)-1)/0.21;} //aerogel density gr/cm^3 parametrization by E.Nappi //trasformation methodes inline static TVector Loc2Area(const TVector2 &x2); //return area affected by hit x2 - inline static Int_t Loc2Sec(const TVector2 &x2); //return sector for given position + inline static Int_t Loc2Sec(const TVector2 &x2); //return sector containing given position + static Int_t Loc2Sec(Double_t x,Double_t y) {return Loc2Sec(TVector2(x,y));} //return sector containing given position inline static TVector Loc2Pad(const TVector2 &x2); //return pad containing given position + static TVector Loc2Pad(Double_t x,Double_t y) {return Loc2Pad(TVector2(x,y));} //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);}//return center of the pad (x,y) inline static Int_t Pad2Sec(const TVector &pad); //return sector of given pad @@ -136,6 +139,7 @@ protected: static Bool_t fgIsTestBeam; //test beam geometry instead of normal RICH flag static Bool_t fgIsWireSag; //wire sagitta ON/OFF flag static Bool_t fgIsResolveClusters; //declustering ON/OFF flag + static Bool_t fgIsFeedback; //generate feedback photon? TObjArray *fpChambers; //list of chambers static Int_t fgHV[6]; //HV applied to anod wires @@ -200,9 +204,9 @@ TVector AliRICHParam::Loc2Pad(const TVector2 &loc) if(sec==1||sec==3||sec==5) pad[0]= Int_t( loc.X() / PadSizeX() )+1; //sector 1 or 3 or 5 else pad[0]=NpadsX() - Int_t( (PcSizeX()-loc.X()) / PadSizeX() ) ; //sector 2 or 4 or 6 //second deal with y - if(sec==1||sec==2) pad[1]=Int_t( loc.Y() / PadSizeY() )+1; //sector 1 or 2 - else if(sec==3||sec==4) pad[1]=Int_t( (loc.Y()-SectorSizeY()-DeadZone()) / PadSizeY() )+NpadsYsec(); //sector 3 or 4 - else pad[1]=NpadsY() - Int_t( (PcSizeY()-loc.Y()) / PadSizeY() ); //sector 5 or 6 + if(sec==1||sec==2) pad[1]=Int_t( loc.Y() / PadSizeY())+1; //sector 1 or 2 + else if(sec==3||sec==4) pad[1]=Int_t( (loc.Y()-SectorSizeY()-DeadZone()) / PadSizeY())+NpadsYsec()+1; //sector 3 or 4 + else pad[1]=NpadsY() - Int_t( (PcSizeY()-loc.Y()) / PadSizeY()); //sector 5 or 6 return pad; } //__________________________________________________________________________________________________ diff --git a/RICH/AliRICHReconstructor.h b/RICH/AliRICHReconstructor.h index 6bd045d8905..4d5c8b71345 100644 --- a/RICH/AliRICHReconstructor.h +++ b/RICH/AliRICHReconstructor.h @@ -17,6 +17,8 @@ public: virtual AliTracker* CreateTracker(AliRunLoader*)const {AliDebug(1,"Start.");return new AliRICHTracker;} //virtual from AliReconstructor virtual void Reconstruct(AliRunLoader* pAL) const; //virtual from AliReconstructor virtual void FillESD(AliRunLoader* pAL, AliESD* pESD) const; //virtual from AliReconstructor + using AliReconstructor::Reconstruct; //to get rid of virtual hidden warning + using AliReconstructor::FillESD; //to get rid of virtual hidden warning private: AliRICH* GetRICH(AliRunLoader* pAL) const; diff --git a/RICH/AliRICHTracker.cxx b/RICH/AliRICHTracker.cxx index 9982e2e446d..d4ea8b39667 100644 --- a/RICH/AliRICHTracker.cxx +++ b/RICH/AliRICHTracker.cxx @@ -23,44 +23,42 @@ Int_t AliRICHTracker::PropagateBack(AliESD *pESD) else RecWithStack(0); AliDebug(1,"Stop pattern recognition"); - return 0; // error code: 0=no error; }//PropagateBack() //__________________________________________________________________________________________________ void AliRICHTracker::RecWithESD(AliESD *pESD) { - //recontruction from ESD - // +//recontruction from ESD- primary way to reconstruct particle ID signal from tracks provided by core detectors + Int_t iNtracks=pESD->GetNumberOfTracks(); Double_t b=GetFieldMap()->SolenoidField()/10;// magnetic field in Tesla AliDebug(1,Form("Start with %i tracks in %f Tesla field",iNtracks,b)); - Double_t xb[3],pb[3];//tmp storage for track parameters - TVector3 x0(0,0,0); TVector3 p0(0,0,0);//tmp storage for AliRICHHelix AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH")); for(Int_t iTrackN=0;iTrackNGetTrack(iTrackN);// get next reconstructed track // if((pTrack->GetStatus()&AliESDtrack::kTOFout)==0) continue; //ignore tracks not recontructed by TOF - pTrack->GetXYZ(xb); - pTrack->GetPxPyPz(pb); +// pTrack->GetXYZ(xb); +// pTrack->GetPxPyPz(pb); Int_t status=pTrack->GetStatus()&AliESDtrack::kTOFout;//get running track parameters Int_t charge = (Int_t)(-TMath::Sign(1.,pTrack->GetSign()*b)); AliDebug(1,Form("Track %i pmod=%f charge=%i stat=%i",iTrackN,pTrack->GetP(),charge,status)); - x0.SetXYZ(xb[0],xb[1],xb[2]); p0.SetXYZ(pb[0],pb[1],pb[2]); - AliRICHHelix helix(x0,p0,charge,b); + AliRICHHelix helix(pTrack->X3(),pTrack->P3(),charge,b); Int_t iChamber=helix.RichIntersect(pRich->P()); AliDebug(1,Form("intersection with %i chamber found",iChamber)); if(!iChamber) continue;//intersection with no chamber found - - Double_t distMip=9999; //min distance between clusters and track position on PC +//find MIP cluster candidate (closest to track intersection point cluster) + Double_t distMip=9999,distX=0,distY=0; //min distance between clusters and track position on PC Int_t iMipId=0; //index of that min distance cluster for(Int_t iClusN=0;iClusNClusters(iChamber)->GetEntries();iClusN++){//clusters loop for intersected chamber AliRICHcluster *pClus=(AliRICHcluster*)pRich->Clusters(iChamber)->UncheckedAt(iClusN);//get pointer to current cluster - Double_t distCurrent=pClus->DistTo(helix.PosPc());//ditance between current cluster and helix intersection with PC + Double_t distCurrent=pClus->DistTo(helix.PosPc());//distance between current cluster and helix intersection with PC if(distCurrentDistX(helix.PosPc()); + distY=pClus->DistY(helix.PosPc()); }//find cluster nearest to the track AliDebug(1,Form("Ploc (%f,%f,%f) dist= %f",helix.Ploc().Mag(),helix.Ploc().Theta()*TMath::RadToDeg(), helix.Ploc().Phi()*TMath::RadToDeg(),pClus->DistTo(helix.PosPc()))); @@ -68,10 +66,15 @@ void AliRICHTracker::RecWithESD(AliESD *pESD) AliDebug(1,Form("Min distance cluster: %i dist is %f",iMipId,distMip)); - AliRICHRecon recon(&helix,pRich->Clusters(iChamber),iMipId); //create reconstruction object for helix, list of clusters - Double_t thetaCerenkov=recon.ThetaCerenkov(); //search for mean Cerenkov angle for this track - AliDebug(1,Form("FINAL Theta Cerenkov=%f",thetaCerenkov)); - pTrack->SetRICHsignal(thetaCerenkov); + AliRICHRecon recon(&helix,pRich->Clusters(iChamber),iMipId); //actual job is done there + + pTrack->SetRICHcluster(iMipId); + pTrack->SetRICHdxdy(distX,distY); + pTrack->SetRICHthetaPhi(helix.Ploc().Theta(),helix.Ploc().Phi()); + pTrack->SetRICHsignal(recon.ThetaCerenkov()); + pTrack->SetRICHnclusters(recon.GetHoughPhotons()); + + AliDebug(1,Form("FINAL Theta Cerenkov=%f",pTrack->GetRICHsignal())); // Double_t richPID[5]={0.2,0.2,0.2,0.2,0.2}; //start with equal probs for (e,mu,pi,k,p) // CalcProb(thetaCerenkov,p0.Mag(),richPID); @@ -83,11 +86,13 @@ void AliRICHTracker::RecWithESD(AliESD *pESD) //__________________________________________________________________________________________________ void AliRICHTracker::RecWithStack(TNtupleD *hn) { -// Reconstruction for particles from STACK +//Reconstruction for particles from STACK. This methode is to be used for RICH standalone when no other detectors are switched on, +//so normal tracking is not available + AliDebug(1,"Start."); AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH")); // pRich->GetLoader()->GetRunLoader()->LoadHeader(); - pRich->GetLoader()->GetRunLoader()->LoadKinematics(); + if(!pRich->GetLoader()->GetRunLoader()->TreeK()) pRich->GetLoader()->GetRunLoader()->LoadKinematics(); AliStack *pStack = pRich->GetLoader()->GetRunLoader()->Stack(); if(!pStack) {AliDebug(1,Form("No STACK found in AliRoot"));return;} Int_t iNtracks=pStack->GetNtrack(); diff --git a/RICH/RichMake b/RICH/RichMake index adbc0ffbe46..6c353453cb8 100644 --- a/RICH/RichMake +++ b/RICH/RichMake @@ -10,7 +10,7 @@ include lib$(Module)rec.pkg SrcRec :=$(SRCS) RootTarget :=$(shell root-config --arch) -DirOut :=$(LIB_MY)/$(Module) +DirOut :=/tmp/$(Module) LibBase :=$(LIB_MY)/lib$(Module)base.so LibSim :=$(LIB_MY)/lib$(Module)sim.so LibRec :=$(LIB_MY)/lib$(Module)rec.so @@ -44,7 +44,7 @@ ifeq ($(RootTarget),linuxicc) LibOpt :=-O -g -shared -Wl else Compiler :=g++ - CompilerOpt :=-g -W -Wall -Werror -fPIC -pipe -fmessage-length=0 -Wno-long-long -pedantic-errors -ansi -I$(shell root-config --incdir) -I$(ALICE_ROOT)/include + CompilerOpt :=-g -W -Wall -Woverloaded-virtual -fPIC -pipe -fmessage-length=0 -Wno-long-long -pedantic-errors -ansi -I$(shell root-config --incdir) -I$(ALICE_ROOT)/include LibOpt :=-O -g -shared -Wl endif @@ -71,6 +71,7 @@ $(LibRec): $(ObjRec) $(DepFile): $(HdrBase) $(HdrSim) $(HdrRec) @[ -d $(DirOut) ] || mkdir -p $(DirOut) + @[ -d $(LIB_MY) ] || mkdir -p $(LIB_MY) @touch $(DepFile) @echo "Generating dependency $@" $(Mute)rmkdepend -f$(DepFile) -p$(DirOut)/ -- $(CompilerOpt) -- $(SrcBase) $(SrcSim) $(SrcRec) 2>/dev/null @@ -117,8 +118,9 @@ show: @echo -e "Rec Library: $(LibRec)\n" spec: $(SrcBase) - @echo $^ - @echo $@ + @echo "^ $^" + @echo "@ $@" + @echo "< $<" clean: @echo "Cleaning..." diff --git a/RICH/RichMenu.C b/RICH/RichMenu.C index ebc8827197d..4a4a8744a75 100644 --- a/RICH/RichMenu.C +++ b/RICH/RichMenu.C @@ -10,14 +10,17 @@ void ps(Int_t event=0) {r->PrintSDigits(event);} //utility print sdigits void pd(Int_t event=0) {r->PrintDigits(event);} //utility print digits void pc(Int_t event=0) {r->PrintClusters(event);}//utility print clusters void pt(Int_t event=0) {r->PrintTracks(event);} //utility print tracks -Int_t ne(Int_t event=0) {AliRICH::Nparticles(kElectron,event,al);} //utility number of electrons -Int_t npi0(Int_t event=0) {AliRICH::Nparticles(kPi0,event,al);} //utility number of electrons -Int_t npip(Int_t event=0) {AliRICH::Nparticles(kPiPlus,event,al);} //utility number of electrons -Int_t npim(Int_t event=0) {AliRICH::Nparticles(kPiMinus,event,al);} //utility number of electrons -Int_t nk0(Int_t event=0) {AliRICH::Nparticles(kK0,event,al);} //utility number of electrons -Int_t nkp(Int_t event=0) {AliRICH::Nparticles(kKPlus,event,al);} //utility number of electrons -Int_t nkm(Int_t event=0) {AliRICH::Nparticles(kKMinus,event,al);} //utility number of electrons -Int_t npp(Int_t event=0) {AliRICH::Nparticles(kProton,event,al);} //utility number of protons +Int_t nem(Int_t event=0) {AliRICHDisplFast::Nparticles(kElectron ,event,al);} //utility number of electrons +Int_t nep(Int_t event=0) {AliRICHDisplFast::Nparticles(kPositron ,event,al);} //utility number of positrons +Int_t nmup(Int_t event=0) {AliRICHDisplFast::Nparticles(kMuonPlus ,event,al);} //utility number of positive muons +Int_t nmum(Int_t event=0) {AliRICHDisplFast::Nparticles(kMuonMinus ,event,al);} //utility number of negative muons +Int_t npi0(Int_t event=0) {AliRICHDisplFast::Nparticles(kPi0 ,event,al);} //utility number of electrons +Int_t npip(Int_t event=0) {AliRICHDisplFast::Nparticles(kPiPlus ,event,al);} //utility number of electrons +Int_t npim(Int_t event=0) {AliRICHDisplFast::Nparticles(kPiMinus ,event,al);} //utility number of electrons +Int_t nk0(Int_t event=0) {AliRICHDisplFast::Nparticles(kK0 ,event,al);} //utility number of electrons +Int_t nkp(Int_t event=0) {AliRICHDisplFast::Nparticles(kKPlus ,event,al);} //utility number of electrons +Int_t nkm(Int_t event=0) {AliRICHDisplFast::Nparticles(kKMinus ,event,al);} //utility number of electrons +Int_t npp(Int_t event=0) {AliRICHDisplFast::Nparticles(kProton ,event,al);} //utility number of protons //__________________________________________________________________________________________________ void pp(int tid) { @@ -73,30 +76,7 @@ void Show() Bool_t isClusters=!rl->LoadRecPoints(); for(Int_t iEventN=0;iEventNGetEventsPerRun();iEventN++){//events loop - Int_t iNparticles=al->Stack()->GetNtrack(); - Int_t iNprims=al->Stack()->GetNprimary(); - Int_t iElectronCounter=0,iPositronCounter=0; - Int_t iPiPlusCounter=0,iPiMinusCounter=0; - Int_t iKPlusCounter=0,iKMinusCounter=0; - Int_t iProtonCounter=0,iProtonBarCounter=0; - Int_t nP=0; - Info("Show-STA"," %i particles to be read",iNparticles); - for(Int_t iParticleN=0;iParticleNStack()->Particle(iParticleN); - nP++; - if(nP%10000==0) Info("Show-STA"," %i particles read",nP); - switch(pPart->GetPdgCode()){ - case kKPlus: iKPlusCounter++; break; - case kKMinus:iKMinusCounter++; break; - case kProton: iProtonCounter++; break; - case kProtonBar:iProtonBarCounter++; break; - case kElectron: iElectronCounter++; break; - case kPositron:iPositronCounter++; break; - case kPiPlus: iPiPlusCounter++; break; - case kPiMinus:iPiMinusCounter++; break; - } - }//stack loop Info("Show-STA","Evt %i-> %i particles %i primaries %i e- %i e+", iEventN, iNparticles, iNprims, iElectronCounter, iPositronCounter); @@ -108,19 +88,6 @@ void Show() iEventN, iNparticles, iNprims, iProtonCounter, iProtonBarCounter); Int_t iHitsCounter=0; - Int_t iNentries=rl->TreeH()->GetEntries(); - for(Int_t iEntryN=0;iEntryNTreeH()->GetEntry(iEntryN);//get current entry (prim) - - for(Int_t iHitN=0;iHitNHits()->GetEntries();iHitN++){//hits loop - iHitsCounter++; - AliRICHhit *pHit = (AliRICHhit*)r->Hits()->At(iHitN);//get current hit - TParticle *pPart=al->Stack()->Particle(pHit->GetTrack());//get stack particle which produced the current hit - }//hits loop - - if(iEntryN<7) Info("Show","Evt %i-> prim %4i has %4i hits from %s (,%7.2f,%7.2f)", - iEventN,iEntryN, r->Hits()->GetEntries(), pPart->GetName(), pPart->Theta()*TMath::RadToDeg(),pPart->Phi()*TMath::RadToDeg()); - }//TreeH loop Info("Show-HIT","Evt %i-> %i particles %i primaries %i entries in TreeH %i hits", iEventN, iNparticles, iNprims, iNentries, iHitsCounter); @@ -202,7 +169,7 @@ void RichMenu() pMenu->AddButton("Show", "Show()", "Shows the structure of events in files"); pMenu->AddButton("Display Fast", "AliRICHDisplFast *d = new AliRICHDisplFast(); d->Exec();", "Display Fast"); pMenu->AddButton("Control Plots", "ControlPlots()", "Create some control histograms"); - + pMenu->AddButton("Recon with stack","r->CheckPR()", "Create RSR.root with ntuple hn"); }else{//it's aliroot, simulate pMenu->AddButton("Debug ON", "DebugON();", "Switch debug on-off"); pMenu->AddButton("Debug OFF", "DebugOFF();", "Switch debug on-off"); @@ -218,60 +185,6 @@ void RichMenu() //__________________________________________________________________________________________________ void DebugOFF(){ Info("DebugOFF",""); AliLog::SetGlobalDebugLevel(0);} void DebugON() { Info("DebugON",""); AliLog::SetGlobalDebugLevel(AliLog::kDebug);} -//__________________________________________________________________________________________________ -TObjArray * CreateHists(Double_t pcut=0.9) -{ - TH2F *pPosH2 =new TH2F("pos" ,"Pos mixture",5,0,5, 9,0,9); pPosH2->SetStats(0); - TH2F *pNegH2 =new TH2F("neg" ,"Neg mixture",5,0,5, 9,0,9); pNegH2->SetStats(0); - TH2F *pPosCutH2 =new TH2F("poscut",Form("Pos mixture with P>%4.2f GeV",pcut),5,0,5, 9,0,9); pPosCutH2->SetStats(0); - TH2F *pNegCutH2 =new TH2F("negcut",Form("Neg mixture with P>%4.2f GeV",pcut),5,0,5, 9,0,9); pNegCutH2->SetStats(0); - pPosH2->GetXaxis()->SetBinLabel(1,"e^{+}"); pNegH2->GetXaxis()->SetBinLabel(1,"e^{-}"); - pPosH2->GetXaxis()->SetBinLabel(2,"#mu^{+}"); pNegH2->GetXaxis()->SetBinLabel(2,"#mu^{-}"); - pPosH2->GetXaxis()->SetBinLabel(3,"#pi^{+}"); pNegH2->GetXaxis()->SetBinLabel(3,"#pi^{-}"); - pPosH2->GetXaxis()->SetBinLabel(4,"K^{+}"); pNegH2->GetXaxis()->SetBinLabel(4,"K^{-}"); - pPosH2->GetXaxis()->SetBinLabel(5,"p^{+}"); pNegH2->GetXaxis()->SetBinLabel(5,"p^{-}"); - - pPosCutH2->GetXaxis()->SetBinLabel(1,"e^{+}"); pNegCutH2->GetXaxis()->SetBinLabel(1,"e^{-}"); - pPosCutH2->GetXaxis()->SetBinLabel(2,"#mu^{+}"); pNegCutH2->GetXaxis()->SetBinLabel(2,"#mu^{-}"); - pPosCutH2->GetXaxis()->SetBinLabel(3,"#pi^{+}"); pNegCutH2->GetXaxis()->SetBinLabel(3,"#pi^{-}"); - pPosCutH2->GetXaxis()->SetBinLabel(4,"K^{+}"); pNegCutH2->GetXaxis()->SetBinLabel(4,"K^{-}"); - pPosCutH2->GetXaxis()->SetBinLabel(5,"p^{+}"); pNegCutH2->GetXaxis()->SetBinLabel(5,"p^{-}"); - - pPosH2->GetYaxis()->SetBinLabel(1,"ch1"); pNegH2->GetYaxis()->SetBinLabel(1,"ch1"); - pPosH2->GetYaxis()->SetBinLabel(2,"ch2"); pNegH2->GetYaxis()->SetBinLabel(2,"ch2"); - pPosH2->GetYaxis()->SetBinLabel(3,"ch3"); pNegH2->GetYaxis()->SetBinLabel(3,"ch3"); - pPosH2->GetYaxis()->SetBinLabel(4,"ch4"); pNegH2->GetYaxis()->SetBinLabel(4,"ch4"); - pPosH2->GetYaxis()->SetBinLabel(5,"ch5"); pNegH2->GetYaxis()->SetBinLabel(5,"ch5"); - pPosH2->GetYaxis()->SetBinLabel(6,"ch6"); pNegH2->GetYaxis()->SetBinLabel(6,"ch6"); - pPosH2->GetYaxis()->SetBinLabel(7,"ch7"); pNegH2->GetYaxis()->SetBinLabel(7,"ch7"); - pPosH2->GetYaxis()->SetBinLabel(8,"prim"); pNegH2->GetYaxis()->SetBinLabel(8,"prim"); - pPosH2->GetYaxis()->SetBinLabel(9,"tot"); pNegH2->GetYaxis()->SetBinLabel(9,"tot"); - - pPosCutH2->GetYaxis()->SetBinLabel(1,"ch1"); pNegCutH2->GetYaxis()->SetBinLabel(1,"ch1"); - pPosCutH2->GetYaxis()->SetBinLabel(2,"ch2"); pNegCutH2->GetYaxis()->SetBinLabel(2,"ch2"); - pPosCutH2->GetYaxis()->SetBinLabel(3,"ch3"); pNegCutH2->GetYaxis()->SetBinLabel(3,"ch3"); - pPosCutH2->GetYaxis()->SetBinLabel(4,"ch4"); pNegCutH2->GetYaxis()->SetBinLabel(4,"ch4"); - pPosCutH2->GetYaxis()->SetBinLabel(5,"ch5"); pNegCutH2->GetYaxis()->SetBinLabel(5,"ch5"); - pPosCutH2->GetYaxis()->SetBinLabel(6,"ch6"); pNegCutH2->GetYaxis()->SetBinLabel(6,"ch6"); - pPosCutH2->GetYaxis()->SetBinLabel(7,"ch7"); pNegCutH2->GetYaxis()->SetBinLabel(7,"ch7"); - pPosCutH2->GetYaxis()->SetBinLabel(8,"prim"); pNegCutH2->GetYaxis()->SetBinLabel(8,"prim"); - pPosCutH2->GetYaxis()->SetBinLabel(9,"tot"); pNegCutH2->GetYaxis()->SetBinLabel(9,"tot"); - TObjArray *pOA=new TObjArray; - pOA->Add(pPosH2);pOA->Add(pNegH2);pOA->Add(pPosCutH2);pOA->Add(pNegCutH2); - return pOA; -}//ParticleContribs() -//__________________________________________________________________________________________________ - -void ShowHists() -{ - pC=new TCanvas("c1","Particle composition"); - pC->Divide(2,2); - pC->cd(1); pos->Draw("text"); gPad->SetGrid(); - pC->cd(2); neg->Draw("text"); gPad->SetGrid(); - pC->cd(3); poscut->Draw("text"); gPad->SetGrid(); - pC->cd(4); negcut->Draw("text"); gPad->SetGrid(); -} - //__________________________________________________________________________________________________ void GeomGui() { @@ -287,3 +200,130 @@ void ControlPlots() r->ControlPlots(); new TBrowser; } +//__________________________________________________________________________________________________ +void ParticleComposition() +{ + TH2F *pFluxH2 =new TH2F("flux","Charged flux for central Hijing event with Vertex<470 P>5MeV for electrons and positrons, P>1GeV for others",10,-5,5, 10,0,10); pFluxH2->SetStats(0); + pFluxH2->GetXaxis()->SetBinLabel(1,Form("p^{-}>%dGeV" ,cutPproton)); + pFluxH2->GetXaxis()->SetBinLabel(2,Form("K^{-}>%dGeV" ,cutPkaonminus)); + pFluxH2->GetXaxis()->SetBinLabel(3,Form("#pi^{-}>%dGeV" , )); + pFluxH2->GetXaxis()->SetBinLabel(4,Form("#mu^{-}>%dGeV" , )); + pFluxH2->GetXaxis()->SetBinLabel(5,Form("e^{+}>%dGeV" ,)); + + pFluxH2->GetXaxis()->SetBinLabel(6,Form("e^{-}>%dGeV")); + pFluxH2->GetXaxis()->SetBinLabel(7,Form("#mu^{+}>%dGeV")); + pFluxH2->GetXaxis()->SetBinLabel(8,Form("#pi^{+}>%dGeV")); + pFluxH2->GetXaxis()->SetBinLabel(9,Form("K^{+}>%dGeV")); + pFluxH2->GetXaxis()->SetBinLabel(10,Form("p^{+}>%dGeV")); + + pFluxH2->GetYaxis()->SetBinLabel(1,"sum"); + pFluxH2->GetYaxis()->SetBinLabel(2,"ch1"); + pFluxH2->GetYaxis()->SetBinLabel(3,"ch2"); + pFluxH2->GetYaxis()->SetBinLabel(4,"ch3"); + pFluxH2->GetYaxis()->SetBinLabel(5,"ch4"); + pFluxH2->GetYaxis()->SetBinLabel(6,"ch5"); + pFluxH2->GetYaxis()->SetBinLabel(7,"ch6"); + pFluxH2->GetYaxis()->SetBinLabel(8,"ch7"); + pFluxH2->GetYaxis()->SetBinLabel(9,"prim"); + pFluxH2->GetYaxis()->SetBinLabel(10,"tot"); + + TH1F *pElecP=new TH1F("Pelec","Electrons made hit in RICH;p [GeV]",1000,-10,10); + TH1F *pMuonP=new TH1F("Pmuon","Muons made hit in RICH;p [GeV]" ,1000,-10,10); + TH1F *pPionP=new TH1F("Ppion","Pions made hit in RICH;p [GeV]" ,1000,-10,10); + TH1F *pKaonP=new TH1F("Pkaon","Kaon made hit in RICH;p [GeV]" ,1000,-10,10); + TH1F *pProtP=new TH1F("Pprot","Protons made hit in RICH;p [GeV]" ,1000,-10,10); + + al->LoadHeader(); + al->LoadKinematics(); + Int_t iNparticles=al->Stack()->GetNtrack(); + Int_t iNprims=al->Stack()->GetNprimary(); + + for(Int_t iParticleN=0;iParticleNStack()->Particle(iParticleN); + + if(iParticleN%10000==0) Info("Show-STA"," %i particles read",iParticleN); + + switch(pPart->GetPdgCode()){ + case kPositron: pFluxH2->Fill(-1,9); break; + case kElectron: pFluxH2->Fill( 0,9); break; + + case kMuonMinus: pFluxH2->Fill(-2,9); break; + case kMuonPlus: pFluxH2->Fill( 1,9); break; + + case kPiMinus: pFluxH2->Fill(-3,9); break; + case kPiPlus: pFluxH2->Fill( 2,9); break; + + case kKMinus: pFluxH2->Fill(-4,9); break; + case kKPlus: pFluxH2->Fill( 3,9); break; + + case kProtonBar: pFluxH2->Fill(-5,9); break; + case kProton: pFluxH2->Fill( 4,9); break; + }//switch + }//stack loop + + + + + rl->LoadHits(); + + + for(Int_t iEntryN=0;iEntryN < rl->TreeH()->GetEntries();iEntryN++){//TreeH loop + rl->TreeH()->GetEntry(iEntryN);//get current entry (prim) + for(Int_t iHitN=0;iHitNHits()->GetEntries();iHitN++){//hits loop + AliRICHhit *pHit = (AliRICHhit*)r->Hits()->At(iHitN);//get current hit + TParticle *pPart=al->Stack()->Particle(pHit->GetTrack());//get stack particle which produced the current hit + + if(TMath::Sqrt(pPart->Vx()*pPart->Vx()+pPart->Vy()*pPart->Vy()+pPart->Vz()*pPart->Vz()) > 470) continue; //cut on vertex position + if(pPart->GetPdgCode()==kElectron && pPart->P()<0.005) continue; //cut on electron momentum 5MeV + if(pPart->GetPdgCode()==kPositron && pPart->P()<0.005) continue; //cut on electron momentum 5MeV + if(pPart->GetPdgCode()!=kElectron && pPart->P()<1) continue; //cut on others momentum 1GeV + + switch(pPart->GetPdgCode()){ + case kPositron : pElecP->Fill(-pPart->P()); pFluxH2->Fill(-1,pHit->C());break; + case kElectron : pElecP->Fill( pPart->P()); pFluxH2->Fill( 0,pHit->C());break; + + case kMuonPlus : pMuonP->Fill( pPart->P()); pFluxH2->Fill(-2,pHit->C());break; + case kMuonMinus: pMuonP->Fill(-pPart->P()); pFluxH2->Fill( 1,pHit->C());break; + + case kPiMinus : pPionP->Fill(-pPart->P()); pFluxH2->Fill(-3,pHit->C());break; + case kPiPlus : pPionP->Fill( pPart->P()); pFluxH2->Fill( 2,pHit->C());break; + + case kKPlus : pKaonP->Fill( pPart->P()); pFluxH2->Fill(-4,pHit->C());break; + case kKMinus : pKaonP->Fill(-pPart->P()); pFluxH2->Fill( 3,pHit->C());break; + + case kProton : pProtP->Fill( pPart->P()); pFluxH2->Fill(-5,pHit->C());break; + case kProtonBar: pProtP->Fill(-pPart->P()); pFluxH2->Fill( 4,pHit->C());break; + } + }//hits loop + }//TreeH loop + + rl->UnloadHits(); + al->UnloadHeader(); + al->UnloadKinematics(); + + + for(Int_t i=1;i<=pFluxH2->GetNbinsX();i++){ + Stat_t sum=0; + for(Int_t j=2;j<=8;j++) sum+=pFluxH2->GetBinContent(i,j); + pFluxH2->SetBinContent(i,1,sum); + } + + TCanvas *pC1=new TCanvas("canvas1",Form("Event Nprims=%i",iNprims),1000,900); + pFluxH2->Draw("text"); gPad->SetGrid(); + + new TCanvas("celec","",200,100); pElecP->Draw(); + new TCanvas("cmuon","",200,100); pMuonP->Draw(); + new TCanvas("cpion","",200,100); pPionP->Draw(); + new TCanvas("ckaon","",200,100); pKaonP->Draw(); + new TCanvas("cprot","",200,100); pProtP->Draw(); + +} + + + + + + + + + diff --git a/RICH/api.txt b/RICH/api.txt index fd84278ec10..05a0e24f063 100644 --- a/RICH/api.txt +++ b/RICH/api.txt @@ -1,8 +1,8 @@ How to open session: - use static methode AliRunLoader::Open("galice.root","AlicE","update") -How to retrive pointer to alice run loader: - use pRICH->GetLoader()->GetRunLoader() (all detector classes inherit from AliDetector wich has GetLoader()) - use methode AliRun::GetRunLoader for gAlice (depricated) + use static method AliRunLoader::Open("galice.root","AlicE","update") +How to retrieve pointer to alice run loader: + use pRICH->GetLoader()->GetRunLoader() (all detector classes inherit from AliDetector which has GetLoader()) + use method AliRun::GetRunLoader for gAlice (deprecated) How to get pointers to different root trees: TreeE belongs to AliRunLoader, available after AliRunLoader::LoadHeader() TreeK belongs to AliRunLoader, available after AliRunLoader::LoadKinematics() @@ -10,37 +10,48 @@ How to get pointers to different root trees: TreeS belongs to AliLoader , available after AliLoader::LoadSDigits() TreeD belongs to AliLoader , available after AliLoader::LoadDigits() TreeR belongs to AliLoader , available after AliLoader::LoadRecPoints() - all methodes return 0 on success. + all methods return 0 on success. +How to get event of interest: + AliRunLoader::GetEvent(event_number) returns 0 on success How to deal with the stack of particles? - first of all, the stack includes primary as well as secondary particles - pointer to the stack is taken: - AliRun::Stack() (global gAlice of type AliRun - depricated way to do) + AliRun::Stack() (global gAlice of type AliRun - deprecated way to do) AliRunLoader::Stack() but before one needs to load event header by AliRunLoader::LoadHeader() otherwise both methods return 0. Moreover loading header gives the information about number of particles only. - To retrive the list of particle one also needs to load kinematics by AliRunLoader::LoadKinematics() + To retrieve the list of particle one also needs to load kinematics by AliRunLoader::LoadKinematics() - total amount of particles in stack for a given event: AliStack::GetNtrack() AliRun::GetEvent() (after LoadHeader()) - total amount of primary particles in stack for a given event (after LoadHeader()): AliStack::GetNprimary() -How to retrive hits: - -Hits are stored on primary by primary basis. Hits for the given primary is TClonesArray. - To retrieve all hits one needs to do: - -initialise the root tree and containers: AliLoader::LoadHits() - -read number of entries in TreeH: TreeH()->GetEntries() - -loop on the list of entires: -How to retrive sdigits? +How to retrieve hits: + Hits are stored on primary by primary basis. Hits for the given primary is TClonesArray. + To retrieve all hits one needs to do: + -initialize the root tree and containers: pRich->GetLoader()->LoadHits(); (AliLoader::LoadHits() returns 0 on success) + -read number of entries in TreeH: pRich->GetLoader()->TreeH()->GetEntries() + -then for each entry: pRich->GetLoader()->TreeH()->GetEntry(i) +How to retrieve sdigits? Sdigits stored in tree S with the branch of TClonesArray, all sdigits in a single TClonesArray So the tree has only one entry. One needs to say: pRich->GetLoader()->LoadSDigits(); this one open file, get the tree and invoke AliRICH::SetTreeAddress() -What are the debug methodes avail: +How to retrieve digits? + Digits stored in tree D with the 7 branches of TClonesArray, one per chamber, all digits of a given chamber in a single TClonesArray + So the tree has only one entry. + -One needs to say: + pRich->GetLoader()->LoadDigits(); this one opens file, gets the tree and invoke AliRICH::SetTreeAddress() which in turn corresponds + branches of the tree to the digits containers in memory. There are 7 containers, one per chamber, all of them belong to AliRICH. + -Then one needs to take the tree entry (only one) to the memory: + pRich->GetLoader()->TreeD()->GetEntry(0) + -Finally pRich->Digits(chamber_number) returns the pointer to TClonesArray of AliRICHdigit +What are the debug methods avail: AliModule::GetDebug() now depricated, use AliDebug printout instead AliModule::SetDebug() AliRun::GetDebug() AliRun::SetDebug() How to get info for a given particle number: - Header and Kinematics trees must be loaded, then possible to retrive pointer to Stack of particles + Header and Kinematics trees must be loaded, then possible to retrieve pointer to Stack of particles Int_t AliRunLoader::LoadHeader(); Int_t AliRunLoader::LoadKinematics() AliStack *AliRunLoader::Stack() TParticle *AliStack::Particle(tid) @@ -53,15 +64,16 @@ What are hte meanings of different VMC flags: gMC->IsTrackAlive() gMC->IsTrackStop() gMC->IsTrackDisappeared() -How is sdigit prodiced from hit: - Responsible methode is AliRICH:: +How is sdigit produced from hit: + Responsible method is AliRICH::Hits2SDigits One hit may affect one or more pads. - Hit position is taken on the anod wires plane as the most of avalanche is developed there. - This poistion is not directly avalable, track intersections with entrance and exit of amplification gap are only store. - So the position in the middle of the gap is calculated as average out of In and Out position. - Then, total charge collected for this hit is calculated by AliRICHParam::Hit2Qdc. + Hit position is taken on the anode wires plane as the most of avalanche is developed there. + This position is not directly available, track intersections with entrance and exit of amplification gap are only stored. + So the position in the middle of the gap is calculated as average out of pHit->In() and pHit->Out() positions. + Then, total charge collected for this hit is calculated by AliRICHParam::Hit2Qdc. + Area of disintegration is a list of pads affected by current hit. This is a parameter of Mathienson How to get pad number for a local position: - use static TVector AliRICHParam::Loc2Pad(TVector2 position) + use static TVector AliRICHParam::Loc2Pad(TVector2 position); Why list of chambers belongs to AliRICHParam: diff --git a/RICH/libRICHsim.pkg b/RICH/libRICHsim.pkg index 74a0084f375..a185647f42c 100644 --- a/RICH/libRICHsim.pkg +++ b/RICH/libRICHsim.pkg @@ -1,5 +1,5 @@ SRCS:= AliRICH.cxx AliRICHv0.cxx AliRICHv1.cxx AliRICHDigitizer.cxx \ - AliGenRadioactive.cxx AliRICHDisplFast.cxx + AliGenRadioactive.cxx AliRICHDisplFast.cxx HDRS:= $(SRCS:.cxx=.h) DHDR:= RICHsimLinkDef.h -- 2.43.0