From 1d46d1f6ab7562f2a85c85fce9268a5eb5bebd3b Mon Sep 17 00:00:00 2001 From: pavlinov Date: Tue, 21 Nov 2006 22:36:04 +0000 Subject: [PATCH] New methods in geometry, update geometry for 3x3 case --- EMCAL/AliEMCALClusterizerv1.cxx | 151 ++++--- EMCAL/AliEMCALClusterizerv1.h | 10 +- EMCAL/AliEMCALGeometry.cxx | 540 +++++++++++++++++++------ EMCAL/AliEMCALGeometry.h | 46 ++- EMCAL/AliEMCALRecPoint.cxx | 120 ++++-- EMCAL/AliEMCALRecPoint.h | 1 + EMCAL/AliEMCALShishKebabTrd1Module.cxx | 84 +++- EMCAL/AliEMCALShishKebabTrd1Module.h | 23 +- EMCAL/AliEMCALv0.cxx | 63 +-- EMCAL/AliEMCALv2.cxx | 19 +- 10 files changed, 770 insertions(+), 287 deletions(-) diff --git a/EMCAL/AliEMCALClusterizerv1.cxx b/EMCAL/AliEMCALClusterizerv1.cxx index 700aacc0857..75919792e9f 100644 --- a/EMCAL/AliEMCALClusterizerv1.cxx +++ b/EMCAL/AliEMCALClusterizerv1.cxx @@ -193,8 +193,11 @@ Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t AbsId) Int_t ieta = -1; Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nTower, nIphi, nIeta) ; - if(!bCell) - Error("DigitizeEnergy","Wrong cell id number") ; + if(!bCell) { + fGeom->PrintGeometry(); + Error("Calibrate()"," Wrong cell id number : %i", AbsId); + assert(0); + } fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nTower,nIphi, nIeta,iphi,ieta); fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi); @@ -248,7 +251,7 @@ void AliEMCALClusterizerv1::Exec(Option_t * option) WriteRecPoints() ; - if(strstr(option,"deb")) + if(strstr(option,"deb") || strstr(option,"all")) PrintRecPoints(option) ; //increment the total number of recpoints per run @@ -412,10 +415,10 @@ void AliEMCALClusterizerv1::InitParameters() fNTowerInGroup = 36; //Produces maximum of 80 pseudoclusters per event - fECAClusteringThreshold = 0.5; // value obtained from Alexei + fECAClusteringThreshold = 0.1; // value obtained from Aleksei fECALocMaxCut = 0.03; // ?? - fECAW0 = 4.5 ; + fECAW0 = 4.5; fTimeGate = 1.e-8 ; fToUnfold = kFALSE ; fRecPointsInRun = 0 ; @@ -427,9 +430,9 @@ void AliEMCALClusterizerv1::InitParameters() //____________________________________________________________________________ Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2) const { - // Gives the neighbourness of two digits = 0 are not neighbour but continue searching + // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching // = 1 are neighbour - // = 2 is in different SM, quit from searching + // = 2 is in different SM; continue searching // neighbours are defined as digits having at least a common vertex // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster // which is compared to a digit (d2) not yet in a cluster @@ -450,7 +453,8 @@ Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d rowdiff = TMath::Abs(iphi1 - iphi2); coldiff = TMath::Abs(ieta1 - ieta2) ; - if (( coldiff <= 1 ) && ( rowdiff <= 1 )) rv = 1; // neighbours with at least commom vertex + // if (( coldiff <= 1 ) && ( rowdiff <= 1 )) rv = 1; // neighbours with at least commom vertex + if(coldiff + rowdiff <= 1) rv = 1; // neighbours with at least commom side; Nov 6,2006 if (gDebug == 2 && rv==1) printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n", @@ -579,7 +583,7 @@ void AliEMCALClusterizerv1::MakeClusters(char* option) aECARecPoints->Clear(); //Start with pseudoclusters, if option - if(strstr(option,"pseudo")) { + if(strstr(option,"pseudo")) { // ?? Nov 3, 2006 TClonesArray * digits = emcalLoader->Digits() ; TClonesArray * digitsC = dynamic_cast(digits->Clone()); @@ -590,7 +594,7 @@ void AliEMCALClusterizerv1::MakeClusters(char* option) int ineb=0; nextdigit.Reset(); - // PseudoClusterization starts + // PseudoClusterization starts while ( (digit = dynamic_cast(nextdigit())) ) { // scan over the list of digitsC recPoint = 0 ; TArrayI clusterECAdigitslist(1000); // what is this @@ -648,17 +652,17 @@ void AliEMCALClusterizerv1::MakeClusters(char* option) //Now do real clusters TClonesArray * digits = emcalLoader->Digits() ; - TClonesArray * digitsC = dynamic_cast(digits->Clone()); + TClonesArray * digitsC = dynamic_cast(digits->Clone()); // will work with thic copy - TIter nextdigit(digitsC) ; + TIter nextdigitC(digitsC) ; AliEMCALDigit * digit; AliEMCALRecPoint * recPoint = 0 ; int ineb=0; - nextdigit.Reset(); + nextdigitC.Reset(); - double e=0.0, ehs = 0.0; - while ( (digit = dynamic_cast(nextdigit())) ) { // clean up digits + double e = 0.0, ehs = 0.0; + while ( (digit = dynamic_cast(nextdigitC())) ) { // clean up digits e = Calibrate(digit->GetAmp(), digit->GetId()); AliEMCALHistoUtilities::FillH1(fHists, 10, digit->GetAmp()); AliEMCALHistoUtilities::FillH1(fHists, 11, e); @@ -673,9 +677,9 @@ void AliEMCALClusterizerv1::MakeClusters(char* option) // Clusterization starts // cout << "Outer Loop" << endl; ineb=0; - nextdigit.Reset(); + nextdigitC.Reset(); recPoint = 0 ; - while ( (digit = dynamic_cast(nextdigit())) ) { // scan over the list of digitsC + while ( (digit = dynamic_cast(nextdigitC())) ) { // scan over the list of digitsC recPoint = 0 ; TArrayI clusterECAdigitslist(1000); // what is this @@ -695,32 +699,33 @@ void AliEMCALClusterizerv1::MakeClusters(char* option) digitsC->Remove(digit) ; AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , thre = %f \n", digit->GetId(), Calibrate(digit->GetAmp(),digit->GetId()), fECAClusteringThreshold)); - nextdigit.Reset(); // will start from beggining + nextdigitC.Reset(); // will start from beggining AliEMCALDigit * digitN = 0; // digi neighbor Int_t index = 0 ; // Find the neighbours + NEIGHBOURS: while (index < iDigitInECACluster){ // scan over digits already in cluster digit = (AliEMCALDigit*)digits->At(clusterECAdigitslist[index]); index++ ; - while ( (digitN = (AliEMCALDigit *)nextdigit())) { // scan over the reduced list of digits + while ( (digitN = (AliEMCALDigit *)nextdigitC())) { // scan over the reduced list of digits + ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!! - switch (ineb ) { - case 0 : // not a neighbour - break ; - case 1 : // are neighbours - recPoint->AddDigit(*digitN, Calibrate(digitN->GetAmp(),digitN->GetId()) ) ; - clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ; - iDigitInECACluster++ ; - digitsC->Remove(digitN) ; - break ; - case 2 : // different SM - break ; - } // switch + if(ineb==1) { + recPoint->AddDigit(*digitN, Calibrate(digitN->GetAmp(),digitN->GetId()) ) ; + clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ; + iDigitInECACluster++ ; + digitsC->Remove(digitN) ; + // Have to start from begining for clusters digits too ; Nov 3,2006 + index = 0; + nextdigitC.Reset() ; + goto NEIGHBOURS; + } // if(ineb==1) + } // scan over the reduced list of digits + nextdigitC.Reset() ; // will start from beginning } // scan over digits already in cluster - nextdigit.Reset() ; // will start from beggining } } // while digit if(recPoint) @@ -815,24 +820,30 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); TObjArray * aECARecPoints = emcalLoader->RecPoints() ; - printf("PrintRecPoints: Clusterization result:") ; + if(strstr(option,"deb")) { + printf("PrintRecPoints: Clusterization result:") ; - printf("event # %d\n", emcalLoader->GetRunLoader()->GetEventNumber() ) ; - printf(" Found %d ECA Rec Points\n ", + printf("event # %d\n", emcalLoader->GetRunLoader()->GetEventNumber() ) ; + printf(" Found %d ECA Rec Points\n ", aECARecPoints->GetEntriesFast()) ; + } fRecPointsInRun += aECARecPoints->GetEntriesFast() ; if(strstr(option,"all")) { - Int_t index =0; - printf("\n-----------------------------------------------------------------------\n") ; - printf("Clusters in ECAL section\n") ; - printf("Index Ene(GeV) Multi Module phi r theta X Y Z Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ; + if(strstr(option,"deb")) { + printf("\n-----------------------------------------------------------------------\n") ; + printf("Clusters in ECAL section\n") ; + printf("Index Ene(GeV) Multi Module GX GY GZ lX lY lZ Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ; + } + Int_t index =0; Float_t maxE=0; Float_t maxL1=0; Float_t maxL2=0; Float_t maxDis=0; + AliEMCALHistoUtilities::FillH1(fHists, 12, double(aECARecPoints->GetEntries())); + for (index = 0 ; index < aECARecPoints->GetEntries() ; index++) { AliEMCALRecPoint * rp = dynamic_cast(aECARecPoints->At(index)) ; TVector3 globalpos; @@ -844,6 +855,7 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) Int_t * primaries; Int_t nprimaries; primaries = rp->GetPrimaries(nprimaries); + if(strstr(option,"deb")) printf("\n%6d %8.4f %3d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ", rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), @@ -861,9 +873,11 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) fPointDis->Fill(rp->GetDispersion()); fPointMult->Fill(rp->GetMultiplicity()); ///////////// - for (Int_t iprimary=0; iprimaryFill(maxE); @@ -871,7 +885,7 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) fMaxL2->Fill(maxL2); fMaxDis->Fill(maxDis); - + if(strstr(option,"deb")) printf("\n-----------------------------------------------------------------------\n"); } } @@ -881,19 +895,20 @@ TList* AliEMCALClusterizerv1::BookHists() gROOT->cd(); - fPointE = new TH1F("pointE","point energy", 2000, 0.0, 150.); - fPointL1 = new TH1F("pointL1","point L1", 1000, 0.0, 3.); - fPointL2 = new TH1F("pointL2","point L2", 1000, 0.0, 3.); - fPointDis = new TH1F("pointDis","point Dis", 1000, 0.0, 3.); - fPointMult = new TH1F("pointMult","point Mult", 100, 0.0, 100.); - fDigitAmp = new TH1F("digitAmp","Digit Amplitude", 2000, 0.0, 5000.); - fMaxE = new TH1F("maxE","Max point energy", 2000, 0.0, 150.); - fMaxL1 = new TH1F("maxL1","Max point L1", 1000, 0.0, 3.); - fMaxL2 = new TH1F("maxL2","Max point L2", 1000, 0.0, 3.); - fMaxDis = new TH1F("maxDis","Max point Dis", 1000, 0.0, 3.); // 9 + fPointE = new TH1F("00_pointE","point energy", 2000, 0.0, 150.); + fPointL1 = new TH1F("01_pointL1","point L1", 1000, 0.0, 3.); + fPointL2 = new TH1F("02_pointL2","point L2", 1000, 0.0, 3.); + fPointDis = new TH1F("03_pointDisp","point dispersion", 1000, 0.0, 10.); + fPointMult = new TH1F("04_pointMult","#cell in point(cluster)", 101, -0.5, 100.5); + fDigitAmp = new TH1F("05_digitAmp","Digit Amplitude", 2000, 0.0, 5000.); + fMaxE = new TH1F("06_maxE","Max point energy", 2000, 0.0, 150.); + fMaxL1 = new TH1F("07_maxL1","Largest (first) of eigenvalue of covariance matrix", 1000, 0.0, 3.); + fMaxL2 = new TH1F("08_maxL2","Smalest (second) of eigenvalue of covariace matrix", 1000, 0.0, 3.); + fMaxDis = new TH1F("09_maxDis","Point dispersion", 1000, 0.0, 10.); // 9 // - new TH1F("adcOfDigits","adc of digits(threshold control)", 1001, -0.5, 1000.5); // 10 - new TH1F("energyOfDigits","energy of digits(threshold control)", 1000, 0.0, 1.); // 11 + new TH1F("10_adcOfDigits","adc of digits(threshold control)", 1001, -0.5, 1000.5); // 10 + new TH1F("11_energyOfDigits","energy of digits(threshold control)", 1000, 0.0, 1.); // 11 + new TH1F("12_numberOfPoints","number of points(clusters)", 101, -0.5, 100.5); // 12 return AliEMCALHistoUtilities::MoveHistsToList("EmcalClusterizerv1ControlHists", kFALSE); } @@ -903,8 +918,34 @@ void AliEMCALClusterizerv1::SaveHists(const char *fn) AliEMCALHistoUtilities::SaveListOfHists(fHists, fn, kTRUE); } +void AliEMCALClusterizerv1::PrintRecoInfo() +{ + printf(" AliEMCALClusterizerv1::PrintRecoInfo() : version %s \n", Version() ); + TH1F *h = (TH1F*)fHists->At(12); + if(h) { + printf(" ## Multiplicity of RecPoints ## \n"); + for(int i=1; i<=h->GetNbinsX(); i++) { + int nbin = int((*h)[i]); + int mult = int(h->GetBinCenter(i)); + if(nbin > 0) printf(" %i : %5.5i %6.3f %% \n", mult, nbin, 100.*nbin/h->GetEntries()); + } + } +} + +void AliEMCALClusterizerv1::DrawLambdasHists() +{ + if(fMaxL1) { + fMaxL1->Draw(); + if(fMaxL2) fMaxL2->Draw("same"); + if(fMaxDis) { + fMaxDis->Draw("same"); + } + } +} + void AliEMCALClusterizerv1::Browse(TBrowser* b) { if(fHists) b->Add(fHists); + if(fGeom) b->Add(fGeom); TTask::Browse(b); } diff --git a/EMCAL/AliEMCALClusterizerv1.h b/EMCAL/AliEMCALClusterizerv1.h index bf33d2abc99..ee2f4f2d9e0 100644 --- a/EMCAL/AliEMCALClusterizerv1.h +++ b/EMCAL/AliEMCALClusterizerv1.h @@ -85,6 +85,8 @@ public: TList* BookHists(); void SaveHists(const char *fn="reco.root"); //*MENU* + void PrintRecoInfo(); //*MENU* + void DrawLambdasHists(); //*MENU* protected: void WriteRecPoints() ; @@ -100,9 +102,9 @@ protected: TH1F* fPointMult; //histogram of point multiplicity TH1F* fDigitAmp; //histogram of digit ADC Amplitude TH1F* fMaxE; //histogram of maximum point energy - TH1F* fMaxL1; //histogram of maximum point L1 - TH1F* fMaxL2; //histogram of maximum point L2 - TH1F* fMaxDis; //histogram of maximum point dispersion + TH1F* fMaxL1; //histogram of largest (first) of eigenvalue of covariance matrix + TH1F* fMaxL2; //histogram of smalest (second) of eigenvalue of covariace matrix + TH1F* fMaxDis; //histogram of point dispersion /////////////////////// @@ -140,7 +142,7 @@ private: Float_t fECAClusteringThreshold ; // minimum energy to seed a EC digit in a cluster Float_t fECALocMaxCut ; // minimum energy difference to distinguish local maxima in a cluster Float_t fECAW0 ; // logarithmic weight for the cluster center of gravity calculation - Int_t fRecPointsInRun ; //! Total number of recpoints in one run + Int_t fRecPointsInRun ; //! Total number of recpoints in one run Float_t fTimeGate ; // Maximum time difference between the digits in ont EMC cluster Float_t fMinECut; // Minimum energy for a digit to be a member of a cluster diff --git a/EMCAL/AliEMCALGeometry.cxx b/EMCAL/AliEMCALGeometry.cxx index c1f9a03a090..9c16ac2c548 100644 --- a/EMCAL/AliEMCALGeometry.cxx +++ b/EMCAL/AliEMCALGeometry.cxx @@ -23,19 +23,22 @@ // -0.7 to 0.7 in eta // Number of Modules and Layers may be controlled by // the name of the instance defined +// EMCAL geometry tree: +// EMCAL -> superModule -> module -> cell +// Indexes +// absId -> nSupMod -> nTower -> (nIphi,nIeta) +// //*-- Author: Sahal Yacoob (LBL / UCT) // and : Yves Schutz (SUBATECH) // and : Jennifer Klay (LBL) // SHASHLYK : Aleksei Pavlinov (WSU) -// SuperModules -> module(or tower) -> cell - +// // --- AliRoot header files --- #include -#include "Riostream.h" +#include #include #include - //#include #include #include #include @@ -43,6 +46,7 @@ #include #include #include +#include // -- ALICE Headers. //#include "AliConst.h" @@ -70,8 +74,10 @@ AliEMCALGeometry::AliEMCALGeometry() fSteelFrontThick(0.),fFrontSteelStrip(0.),fLateralSteelStrip(0.),fPassiveScintThick(0.),fPhiModuleSize(0.), fEtaModuleSize(0.),fPhiTileSize(0.),fEtaTileSize(0.),fLongModuleSize(0.),fNPhiSuperModule(0),fNPHIdiv(0),fNETAdiv(0), fNCells(0),fNCellsInSupMod(0),fNCellsInTower(0),fNTRU(0),fNTRUEta(0),fNTRUPhi(0),fTrd1Angle(0.),f2Trd1Dx2(0.), - fPhiGapForSM(0.),fKey110DEG(0),fTrd2AngleY(0.),f2Trd2Dy2(0.),fEmptySpace(0.),fTubsR(0.),fTubsTurnAngle(0.),fEtaCentersOfCells(0), - fXCentersOfCells(0),fPhiCentersOfCells(0),fShishKebabTrd1Modules(),fNAdditionalOpts(0) + fPhiGapForSM(0.),fKey110DEG(0),fPhiBoundariesOfSM(0), fPhiCentersOfSM(0),fEtaMaxOfTRD1(0), + fTrd2AngleY(0.),f2Trd2Dy2(0.),fEmptySpace(0.),fTubsR(0.),fTubsTurnAngle(0.),fCentersOfCellsEtaDir(0), + fCentersOfCellsXDir(0),fCentersOfCellsPhiDir(0),fEtaCentersOfCells(0),fPhiCentersOfCells(0), + fShishKebabTrd1Modules(0),fNAdditionalOpts(0) { // default ctor only for internal usage (singleton) // must be kept public for root persistency purposes, but should never be called by the outside world @@ -87,13 +93,22 @@ AliEMCALGeometry::AliEMCALGeometry(const Text_t* name, const Text_t* title) fSteelFrontThick(0.),fFrontSteelStrip(0.),fLateralSteelStrip(0.),fPassiveScintThick(0.),fPhiModuleSize(0.), fEtaModuleSize(0.),fPhiTileSize(0.),fEtaTileSize(0.),fLongModuleSize(0.),fNPhiSuperModule(0),fNPHIdiv(0),fNETAdiv(0), fNCells(0),fNCellsInSupMod(0),fNCellsInTower(0),fNTRU(0),fNTRUEta(0),fNTRUPhi(0),fTrd1Angle(0.),f2Trd1Dx2(0.), - fPhiGapForSM(0.),fKey110DEG(0),fTrd2AngleY(0.),f2Trd2Dy2(0.),fEmptySpace(0.),fTubsR(0.),fTubsTurnAngle(0.),fEtaCentersOfCells(0), - fXCentersOfCells(0),fPhiCentersOfCells(0),fShishKebabTrd1Modules(),fNAdditionalOpts(0) + fPhiGapForSM(0.),fKey110DEG(0),fPhiBoundariesOfSM(0), fPhiCentersOfSM(0), fEtaMaxOfTRD1(0), + fTrd2AngleY(0.),f2Trd2Dy2(0.),fEmptySpace(0.),fTubsR(0.),fTubsTurnAngle(0.),fCentersOfCellsEtaDir(0), + fCentersOfCellsXDir(0),fCentersOfCellsPhiDir(0),fEtaCentersOfCells(0),fPhiCentersOfCells(0), + fShishKebabTrd1Modules(0),fNAdditionalOpts(0) { // ctor only for internal usage (singleton) AliDebug(2, Form("AliEMCALGeometry(%s,%s) ", name,title)); + Init(); + CreateListOfTrd1Modules(); + + if (AliDebugLevel()>=2) { + PrintGeometry(); + } + } //______________________________________________________________________ AliEMCALGeometry::AliEMCALGeometry(const AliEMCALGeometry& geom) @@ -138,13 +153,18 @@ AliEMCALGeometry::AliEMCALGeometry(const AliEMCALGeometry& geom) f2Trd1Dx2(geom.f2Trd1Dx2), fPhiGapForSM(geom.fPhiGapForSM), fKey110DEG(geom.fKey110DEG), + fPhiBoundariesOfSM(geom.fPhiBoundariesOfSM), + fPhiCentersOfSM(geom.fPhiCentersOfSM), + fEtaMaxOfTRD1(geom.fEtaMaxOfTRD1), fTrd2AngleY(geom.fTrd2AngleY), f2Trd2Dy2(geom.f2Trd2Dy2), fEmptySpace(geom.fEmptySpace), fTubsR(geom.fTubsR), fTubsTurnAngle(geom.fTubsTurnAngle), + fCentersOfCellsEtaDir(geom.fCentersOfCellsEtaDir), + fCentersOfCellsXDir(geom.fCentersOfCellsXDir), + fCentersOfCellsPhiDir(geom.fCentersOfCellsPhiDir), fEtaCentersOfCells(geom.fEtaCentersOfCells), - fXCentersOfCells(geom.fXCentersOfCells), fPhiCentersOfCells(geom.fPhiCentersOfCells), fShishKebabTrd1Modules(geom.fShishKebabTrd1Modules), fNAdditionalOpts(geom.fNAdditionalOpts) @@ -167,6 +187,9 @@ void AliEMCALGeometry::Init(void){ // SHISH_25 or SHISH_62 // 11-oct-05 - correction for pre final design // Feb 06,2006 - decrease the weight of EMCAL + // + // Oct 30,2006 - SHISH_TRD1_CURRENT_1X1, SHISH_TRD1_CURRENT_2X2 or SHISH_TRD1_CURRENT_3X3; + // fAdditionalOpts[0] = "nl="; // number of sampling layers (fNECLayers) fAdditionalOpts[1] = "pbTh="; // cm, Thickness of the Pb (fECPbRadThick) @@ -179,14 +202,14 @@ void AliEMCALGeometry::Init(void){ fGeoName = GetName(); fGeoName.ToUpper(); fKey110DEG = 0; - if(fGeoName.Contains("110DEG")) fKey110DEG = 1; // for GetAbsCellId + if(fGeoName.Contains("110DEG") || fGeoName.Contains("CURRENT")) fKey110DEG = 1; // for GetAbsCellId fShishKebabTrd1Modules = 0; fTrd2AngleY = f2Trd2Dy2 = fEmptySpace = fTubsR = fTubsTurnAngle = 0; fNZ = 114; // granularity along Z (eta) fNPhi = 168; // granularity in phi (azimuth) - fArm1PhiMin = 60.0; // degrees, Starting EMCAL Phi position - fArm1PhiMax = 180.0; // degrees, Ending EMCAL Phi position + fArm1PhiMin = 80.0; // degrees, Starting EMCAL Phi position + fArm1PhiMax = 190.0; // degrees, Ending EMCAL Phi position fArm1EtaMin = -0.7; // pseudorapidity, Starting EMCAL Eta position fArm1EtaMax = +0.7; // pseudorapidity, Ending EMCAL Eta position fIPDistance = 454.0; // cm, Radial distance to inner surface of EMCAL @@ -229,7 +252,7 @@ void AliEMCALGeometry::Init(void){ if(fGeoName.Contains("TRD1")) { // 30-jan-05 // for final design fPhiGapForSM = 2.; // cm, only for final TRD1 geometry - if(fGeoName.Contains("MAY05") || fGeoName.Contains("WSUC") || fGeoName.Contains("FINAL")){ + if(fGeoName.Contains("MAY05") || fGeoName.Contains("WSUC") || fGeoName.Contains("FINAL") || fGeoName.Contains("CURRENT")){ fNumberOfSuperModules = 12; // 20-may-05 if(fGeoName.Contains("WSUC")) fNumberOfSuperModules = 1; // 27-may-05 fNECLayers = 77; // (13-may-05 from V.Petrov) @@ -242,13 +265,21 @@ void AliEMCALGeometry::Init(void){ fNZ = 24; fTrd1Angle = 1.5; // 1.3 or 1.5 - if(fGeoName.Contains("FINAL")) { // 9-sep-05 + if(fGeoName.Contains("FINAL") || fGeoName.Contains("CURRENT")) { // 9-sep-05 fNumberOfSuperModules = 10; - if(fGeoName.Contains("110DEG")) { + if(GetKey110DEG()) { fNumberOfSuperModules = 12;// last two modules have size 10 degree in phi (180 indexes fNCellsInTower = fNPHIdiv*fNETAdiv; fNCellsInSupMod = fNCellsInTower*fNPhi*fNZ; fNCells = fNCellsInSupMod*fNumberOfSuperModules; - if(fGeoName.Contains("110DEG")) fNCells -= fNCellsInSupMod; + if(GetKey110DEG()) fNCells -= fNCellsInSupMod; fLongModuleSize = fNECLayers*(fECScintThick + fECPbRadThickness); if(fGeoName.Contains("MAY05")) fLongModuleSize += (fFrontSteelStrip + fPassiveScintThick); @@ -308,7 +339,7 @@ void AliEMCALGeometry::Init(void){ fNPhiSuperModule = fNumberOfSuperModules/2; if(fNPhiSuperModule<1) fNPhiSuperModule = 1; - //There is always one more scintillator than radiator layer because of the first block of aluminium + fShellThickness = fAlFrontThick + fGap2Active + fNECLayers*GetECScintThick()+(fNECLayers-1)*GetECPbRadThick(); if(fGeoName.Contains("SHISH")) { fShellThickness = fSteelFrontThick + fLongModuleSize; @@ -331,62 +362,129 @@ void AliEMCALGeometry::Init(void){ fEnvelop[2] = 1.00001*fZLength; // add some padding for mother volume. fNumberOfSuperModules = 12; - + + // SM phi boundaries - (0,1),(2,3) .. (10,11) - has the same boundaries; Nov 7, 2006 + fPhiBoundariesOfSM.Set(fNumberOfSuperModules); + fPhiCentersOfSM.Set(fNumberOfSuperModules/2); + fPhiBoundariesOfSM[0] = TMath::PiOver2() - TMath::ATan2(fParSM[1] , fIPDistance); // 1th and 2th modules) + fPhiBoundariesOfSM[1] = TMath::PiOver2() + TMath::ATan2(fParSM[1] , fIPDistance); + fPhiCentersOfSM[0] = TMath::PiOver2(); + for(int i=1; i<=4; i++) { // from 2th ro 9th + fPhiBoundariesOfSM[2*i] = fPhiBoundariesOfSM[0] + 20.*TMath::DegToRad()*i; + fPhiBoundariesOfSM[2*i+1] = fPhiBoundariesOfSM[1] + 20.*TMath::DegToRad()*i; + fPhiCentersOfSM[i] = fPhiCentersOfSM[0] + 20.*TMath::DegToRad()*i; + } + fPhiBoundariesOfSM[11] = 190.*TMath::DegToRad(); + fPhiBoundariesOfSM[10] = fPhiBoundariesOfSM[11] - TMath::ATan2((fParSM[1]) , fIPDistance); + fPhiCentersOfSM[5] = (fPhiBoundariesOfSM[10]+fPhiBoundariesOfSM[11])/2.; + fgInit = kTRUE; - if (AliDebugLevel()>=2) { - printf("Init: geometry of EMCAL named %s is as follows:\n", fGeoName.Data()); - printf( " ECAL : %d x (%f cm Pb, %f cm Sc) \n", - GetNECLayers(), GetECPbRadThick(), GetECScintThick() ) ; - printf(" fSampling %5.2f \n", fSampling ); - if(fGeoName.Contains("SHISH")){ - printf(" fIPDistance %6.3f cm \n", fIPDistance); - if(fSteelFrontThick>0.) - printf(" fSteelFrontThick %6.3f cm \n", fSteelFrontThick); - printf(" fNPhi %i | fNZ %i \n", fNPhi, fNZ); - printf(" fNCellsInTower %i : fNCellsInSupMod %i : fNCells %i\n",fNCellsInTower, fNCellsInSupMod, fNCells); - if(fGeoName.Contains("MAY05")){ - printf(" fFrontSteelStrip %6.4f cm (thickness of front steel strip)\n", - fFrontSteelStrip); - printf(" fLateralSteelStrip %6.4f cm (thickness of lateral steel strip)\n", - fLateralSteelStrip); - printf(" fPassiveScintThick %6.4f cm (thickness of front passive Sc tile)\n", - fPassiveScintThick); - } - printf(" X:Y module size %6.3f , %6.3f cm \n", fPhiModuleSize, fEtaModuleSize); - printf(" X:Y tile size %6.3f , %6.3f cm \n", fPhiTileSize, fEtaTileSize); - printf(" #of sampling layers %i(fNECLayers) \n", fNECLayers); - printf(" fLongModuleSize %6.3f cm \n", fLongModuleSize); - printf(" #supermodule in phi direction %i \n", fNPhiSuperModule ); - } - if(fGeoName.Contains("TRD")) { - printf(" fTrd1Angle %7.4f\n", fTrd1Angle); - printf(" f2Trd1Dx2 %7.4f\n", f2Trd1Dx2); - if(fGeoName.Contains("TRD2")) { - printf(" fTrd2AngleY %7.4f\n", fTrd2AngleY); - printf(" f2Trd2Dy2 %7.4f\n", f2Trd2Dy2); - printf(" fTubsR %7.2f cm\n", fTubsR); - printf(" fTubsTurnAngle %7.4f\n", fTubsTurnAngle); - printf(" fEmptySpace %7.4f cm\n", fEmptySpace); - } else if(fGeoName.Contains("TRD1") && fGeoName.Contains("FINAL")){ - printf("SM dimensions(TRD1) : dx %7.2f dy %7.2f dz %7.2f (SMOD, BOX)\n", - fParSM[0],fParSM[1],fParSM[2]); - printf(" fPhiGapForSM %7.4f cm \n", fPhiGapForSM); - if(fGeoName.Contains("110DEG"))printf(" Last two modules have size 10 degree in phi (180 for EMCAL envelope only\n", + GetArm1PhiMin(), GetArm1PhiMax(),GetArm1EtaMin(), GetArm1EtaMax(), GetIPDistance() ); + printf( " ECAL : %d x (%f cm Pb, %f cm Sc) \n", + GetNECLayers(), GetECPbRadThick(), GetECScintThick() ) ; + printf(" fSampling %5.2f \n", fSampling ); + if(fGeoName.Contains("SHISH")){ + printf(" fIPDistance %6.3f cm \n", fIPDistance); + if(fSteelFrontThick>0.) + printf(" fSteelFrontThick %6.3f cm \n", fSteelFrontThick); + printf(" fNPhi %i | fNZ %i \n", fNPhi, fNZ); + printf(" fNCellsInTower %i : fNCellsInSupMod %i : fNCells %i\n",fNCellsInTower, fNCellsInSupMod, fNCells); + if(fGeoName.Contains("MAY05")){ + printf(" fFrontSteelStrip %6.4f cm (thickness of front steel strip)\n", + fFrontSteelStrip); + printf(" fLateralSteelStrip %6.4f cm (thickness of lateral steel strip)\n", + fLateralSteelStrip); + printf(" fPassiveScintThick %6.4f cm (thickness of front passive Sc tile)\n", + fPassiveScintThick); + } + printf(" X:Y module size %6.3f , %6.3f cm \n", fPhiModuleSize, fEtaModuleSize); + printf(" X:Y tile size %6.3f , %6.3f cm \n", fPhiTileSize, fEtaTileSize); + printf(" #of sampling layers %i(fNECLayers) \n", fNECLayers); + printf(" fLongModuleSize %6.3f cm \n", fLongModuleSize); + printf(" #supermodule in phi direction %i \n", fNPhiSuperModule ); + } + if(fGeoName.Contains("TRD")) { + printf(" fTrd1Angle %7.4f\n", fTrd1Angle); + printf(" f2Trd1Dx2 %7.4f\n", f2Trd1Dx2); + if(fGeoName.Contains("TRD2")) { + printf(" fTrd2AngleY %7.4f\n", fTrd2AngleY); + printf(" f2Trd2Dy2 %7.4f\n", f2Trd2Dy2); + printf(" fTubsR %7.2f cm\n", fTubsR); + printf(" fTubsTurnAngle %7.4f\n", fTubsTurnAngle); + printf(" fEmptySpace %7.4f cm\n", fEmptySpace); + } else if(fGeoName.Contains("TRD1")){ + printf("SM dimensions(TRD1) : dx %7.2f dy %7.2f dz %7.2f (SMOD, BOX)\n", + fParSM[0],fParSM[1],fParSM[2]); + printf(" fPhiGapForSM %7.4f cm (%7.4f <- phi size in degree)\n", + fPhiGapForSM, TMath::ATan2(fPhiGapForSM,fIPDistance)*TMath::RadToDeg()); + if(GetKey110DEG()) printf(" Last two modules have size 10 degree in phi (180 %7.5f(%7.2f) : center %7.5f(%7.2f) \n", i, + fPhiBoundariesOfSM[2*i], fPhiBoundariesOfSM[2*i]*TMath::RadToDeg(), + fPhiBoundariesOfSM[2*i+1], fPhiBoundariesOfSM[2*i+1]*TMath::RadToDeg(), + fPhiCentersOfSM[i], fPhiCentersOfSM[i]*TMath::RadToDeg()); + } + printf(" fShishKebabTrd1Modules has %i modules : max eta %5.4f \n", + fShishKebabTrd1Modules->GetSize(),fEtaMaxOfTRD1); + + printf("\n Cells grid in eta directions : size %i\n", fCentersOfCellsEtaDir.GetSize()); + for(Int_t i=0; i nSupMod %i nTower %i nIphi %i nIeta %i \n", tit, absId, nSupMod, nTower, nIphi, nIeta); + if(pri>0) { + GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta); + printf(" local SM index : iphi %i : ieta %i \n", iphi,ieta); + GetGlobal(absId, vg); + printf(" vglob : mag %7.2f : perp %7.2f : z %7.2f : eta %6.4f : phi %6.4f(%6.2f) \n", + vg.Mag(), vg.Perp(), vg.Z(), vg.Eta(), vg.Phi(), vg.Phi()*TMath::RadToDeg()); + } +} + +//______________________________________________________________________ void AliEMCALGeometry::CheckAdditionalOptions() { // Feb 06,2006 @@ -704,42 +802,80 @@ Bool_t AliEMCALGeometry::GetCellIndex(Int_t absId,Int_t &nSupMod,Int_t &nTower,I void AliEMCALGeometry::GetModulePhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower, int &iphim, int &ietam) const { - // added nSupMod; have to check - 19-oct-05 ! + // added nSupMod; - 19-oct-05 ! // Alice numbering scheme - Jun 01,2006 + // ietam, iphi - indexes of module in two dimensional grid of SM + // ietam - have to change from 0 to fNZ-1 + // iphim - have to change from 0 to nphi-1 (fNPhi-1 or fNPhi/2-1) static Int_t nphi; if(fKey110DEG == 1 && nSupMod>=10) nphi = fNPhi/2; else nphi = fNPhi; - ietam = nTower/nphi; // have to change from 0 to fNZ-1 - iphim = nTower%nphi; // have to change from 0 to fNPhi-1 + ietam = nTower/nphi; + iphim = nTower%nphi; } void AliEMCALGeometry::GetCellPhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower, Int_t nIphi, Int_t nIeta, int &iphi, int &ieta) const { - // added nSupMod; Nov 25, 05 - // Alice numbering scheme - Jun 01,2006 + // + // Added nSupMod; Nov 25, 05 + // Alice numbering scheme - Jun 01,2006 + // ieta, iphi - indexes of cell(tower) in two dimensional grid of SM + // ieta - have to change from 0 to (fNZ*fNETAdiv-1) + // iphi - have to change from 0 to (fNPhi*fNPHIdiv-1 or fNPhi*fNPHIdiv/2-1) + // static Int_t iphim, ietam; GetModulePhiEtaIndexInSModule(nSupMod,nTower, iphim, ietam); - // have to change from 0 to (fNZ*fNETAdiv-1) - ieta = ietam*fNETAdiv + (1-nIeta); // x(module) = -z(SM) - // iphi - have to change from 0 to (fNPhi*fNPHIdiv-1) + // ieta = ietam*fNETAdiv + (1-nIeta); // x(module) = -z(SM) + ieta = ietam*fNETAdiv + (fNETAdiv - 1 - nIeta); // x(module) = -z(SM) iphi = iphim*fNPHIdiv + nIphi; // y(module) = y(SM) + + if(iphi<0 || ieta<0) + AliDebug(1,Form(" nSupMod %i nTower %i nIphi %i nIeta %i => ieta %i iphi %i\n", + nSupMod, nTower, nIphi, nIeta, ieta, iphi)); } Int_t AliEMCALGeometry::GetSuperModuleNumber(Int_t absId) const { - //return the number of the - //supermodule given the absolute - //ALICE numbering + // Return the number of the supermodule given the absolute + // ALICE numbering id static Int_t nSupMod, nTower, nIphi, nIeta; GetCellIndex(absId, nSupMod, nTower, nIphi, nIeta); return nSupMod; } +void AliEMCALGeometry::GetModuleIndexesFromCellIndexesInSModule(Int_t nSupMod, Int_t iphi, Int_t ieta, + Int_t &iphim, Int_t &ietam, Int_t &nTower) const +{ + // Transition from cell indexes (ieta,iphi) to module indexes (ietam,iphim, nTower) + static Int_t nphi; + nphi = GetNumberOfModuleInPhiDirection(nSupMod); + + ietam = ieta/fNETAdiv; + iphim = iphi/fNPHIdiv; + nTower = ietam * nphi + iphim; +} + +Int_t AliEMCALGeometry::GetAbsCellIdFromCellIndexes(Int_t nSupMod, Int_t iphi, Int_t ieta) const +{ + // Transition from super module number(nSupMod) and cell indexes (ieta,iphi) to absId + static Int_t ietam, iphim, nTower; + static Int_t nIeta, nIphi; // cell indexes in module + + GetModuleIndexesFromCellIndexesInSModule(nSupMod, iphi, ieta, ietam, iphim, nTower); + + nIeta = ieta%fNETAdiv; + nIeta = fNETAdiv - 1 - nIeta; + nIphi = iphi%fNPHIdiv; + + return GetAbsCellId(nSupMod, nTower, nIphi, nIeta); +} + + // Methods for AliEMCALRecPoint - Feb 19, 2006 Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t &xr, Double_t &yr, Double_t &zr) const { @@ -755,16 +891,16 @@ Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t &xr, Double_t GetCellIndex(absId, nSupMod, nTower, nIphi, nIeta); GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi, ieta); - xr = fXCentersOfCells.At(ieta); - zr = fEtaCentersOfCells.At(ieta); + xr = fCentersOfCellsXDir.At(ieta); + zr = fCentersOfCellsEtaDir.At(ieta); if(nSupMod<10) { - yr = fPhiCentersOfCells.At(iphi); + yr = fCentersOfCellsPhiDir.At(iphi); } else { - yr = fPhiCentersOfCells.At(iphi + phiIndexShift); - // cout<<" absId "<SetName("ListOfTRD1"); for(int iz=0; iz< GetNZ(); iz++) { if(iz==0) { mod = new AliEMCALShishKebabTrd1Module(TMath::Pi()/2.,this); @@ -815,47 +952,89 @@ void AliEMCALGeometry::CreateListOfTrd1Modules() } else { AliDebug(2,Form(" Already exits : ")); } - AliDebug(2,Form(" fShishKebabTrd1Modules has %i modules \n", - fShishKebabTrd1Modules->GetSize())); + mod = (AliEMCALShishKebabTrd1Module*)fShishKebabTrd1Modules->At(fShishKebabTrd1Modules->GetSize()-1); + fEtaMaxOfTRD1 = mod->GetMaxEtaOfModule(0); + + AliDebug(2,Form(" fShishKebabTrd1Modules has %i modules : max eta %5.4f \n", + fShishKebabTrd1Modules->GetSize(),fEtaMaxOfTRD1)); // Feb 20,2006; // Jun 01, 2006 - ALICE numbering scheme // define grid for cells in eta(z) and x directions in local coordinates system of SM - // fEtaCentersOfCells = new TArrayD(fNZ *fNETAdiv); - // fXCentersOfCells = new TArrayD(fNZ *fNETAdiv); - fEtaCentersOfCells.Set(fNZ *fNETAdiv); - fXCentersOfCells.Set(fNZ *fNETAdiv); - AliDebug(2,Form(" Cells grid in eta directions : size %i\n", fEtaCentersOfCells.GetSize())); + // Works just for 2x2 case only -- ?? start here + // + // + // Define grid for cells in phi(y) direction in local coordinates system of SM + // as for 2X2 as for 3X3 - Nov 8,2006 + // + AliDebug(2,Form(" Cells grid in phi directions : size %i\n", fCentersOfCellsPhiDir.GetSize())); + Int_t ind=0; // this is phi index Int_t iphi=0, ieta=0, nTower=0; - Double_t xr, zr; - for(Int_t it=0; itGetCenterOfCellInLocalCoordinateofSM(ic, xr, zr); - GetCellPhiEtaIndexInSModule(0, nTower, 0, ic, iphi, ieta); // don't depend from phi - ieta in action - fXCentersOfCells.AddAt(float(xr) - fParSM[0],ieta); - fEtaCentersOfCells.AddAt(float(zr) - fParSM[2],ieta); + for(Int_t ic=0; icGetCenterOfCellInLocalCoordinateofSM(ic, xr, zr); // case of 2X2 + GetCellPhiEtaIndexInSModule(0, nTower, 0, ic, iphi, ieta); // don't depend from phi - ieta in action + fCentersOfCellsXDir.AddAt(float(xr) - fParSM[0],ieta); + fCentersOfCellsEtaDir.AddAt(float(zr) - fParSM[2],ieta); + } if(fNPHIdiv==3) { + trd1->GetCenterOfCellInLocalCoordinateofSM_3X3(ic, xr, zr); // case of 3X3 + GetCellPhiEtaIndexInSModule(0, nTower, 0, ic, iphi, ieta); // don't depend from phi - ieta in action + fCentersOfCellsXDir.AddAt(float(xr) - fParSM[0],ieta); + fCentersOfCellsEtaDir.AddAt(float(zr) - fParSM[2],ieta); + } + // Define grid on eta direction for each bin in phi + for(int iphi=0; iphiGetRadius(); + y = fCentersOfCellsPhiDir[iphi]; + r = TMath::Sqrt(x*x + y*y + zr*zr); + theta = TMath::ACos(zr/r); + eta = AliEMCALShishKebabTrd1Module::ThetaToEta(theta); + // ind = ieta*fCentersOfCellsPhiDir.GetSize() + iphi; + ind = iphi*fCentersOfCellsEtaDir.GetSize() + ieta; + fEtaCentersOfCells.AddAt(eta, ind); + } + //printf(" ieta %i : xr + trd1->GetRadius() %f : zr %f : eta %f \n", ieta, xr + trd1->GetRadius(), zr, eta); } } - for(Int_t i=0; i11) return kFALSE; + i = nSupMod/2; + phiMin = fPhiBoundariesOfSM[2*i]; + phiMax = fPhiBoundariesOfSM[2*i+1]; + return kTRUE; +} + +Bool_t AliEMCALGeometry::GetPhiBoundariesOfSMGap(Int_t nPhiSec, Double_t &phiMin, Double_t &phiMax) const +{ + // 0<= nPhiSec <=4; phi in rad + // 0; gap boundaries between 0th&2th | 1th&3th SM + // 1; gap boundaries between 2th&4th | 3th&5th SM + // 2; gap boundaries between 4th&6th | 5th&7th SM + // 3; gap boundaries between 6th&8th | 7th&9th SM + // 4; gap boundaries between 8th&10th | 9th&11th SM + if(nPhiSec<0 || nPhiSec >4) return kFALSE; + phiMin = fPhiBoundariesOfSM[2*nPhiSec+1]; + phiMax = fPhiBoundariesOfSM[2*nPhiSec+2]; + return kTRUE; +} + +Bool_t AliEMCALGeometry::SuperModuleNumberFromEtaPhi(Double_t eta, Double_t phi, Int_t &nSupMod) const +{ + // Return false if phi belongs a phi cracks between SM + + static Int_t i; + + if(TMath::Abs(eta) > fEtaMaxOfTRD1) return kFALSE; + + phi = TVector2::Phi_0_2pi(phi); // move phi to (0,2pi) boundaries + for(i=0; i<6; i++) { + if(phi>=fPhiBoundariesOfSM[2*i] && phi<=fPhiBoundariesOfSM[2*i+1]) { + nSupMod = 2*i; + if(eta < 0.0) nSupMod++; + // Info("SuperModuleNumberFromEtaPhi", Form(" nSupMod %i : eta %5.3f : phi %5.3f(%5.1f) ", + // nSupMod, eta, phi, phi*TMath::RadToDeg())); + return kTRUE; + } + } + AliDebug(1,Form("eta %f phi %f(%5.2f) : nSupMod %i : #bound %i", eta,phi,phi*TMath::RadToDeg(), nSupMod,i)); + return kFALSE; +} + +Bool_t AliEMCALGeometry::GetAbsCellIdFromEtaPhi(Double_t eta, Double_t phi, Int_t &absId) const +{ + // Nov 17,2006 + // stay here - phi problem as usual + static Int_t nSupMod, i, ieta, iphi, etaShift, nphi; + static Double_t absEta=0.0, d=0.0, dmin=0.0, phiLoc; + absId = nSupMod = - 1; + if(SuperModuleNumberFromEtaPhi(eta, phi, nSupMod)) { + // phi index first + phi = TVector2::Phi_0_2pi(phi); + phiLoc = phi - fPhiCentersOfSM[nSupMod/2]; + nphi = fPhiCentersOfCells.GetSize(); + if(nSupMod>=10) { + phiLoc = phi - 190.*TMath::DegToRad(); + nphi /= 2; + } + + dmin = TMath::Abs(fPhiCentersOfCells[0]-phiLoc); + iphi = 0; + for(i=1; iAdd(fShishKebabTrd1Modules); +} + +Bool_t AliEMCALGeometry::IsFolder() const +{ + if(fShishKebabTrd1Modules) return kTRUE; + else return kFALSE; +} diff --git a/EMCAL/AliEMCALGeometry.h b/EMCAL/AliEMCALGeometry.h index c86f9220696..db79b9e9814 100644 --- a/EMCAL/AliEMCALGeometry.h +++ b/EMCAL/AliEMCALGeometry.h @@ -41,7 +41,11 @@ public: Fatal("operator =", "not implemented"); return *this; }; - + void PrintGeometry(); //*MENU* + void PrintCellIndexes(Int_t absId=0, int pri=0, char *tit=""); //*MENU* + virtual void Browse(TBrowser* b); + virtual Bool_t IsFolder() const; + void FillTRU(const TClonesArray * digits, TClonesArray * amptru, TClonesArray * timeRtru) ; //Fills Trigger Unit matrices with digit amplitudes and time void GetCellPhiEtaIndexInSModuleFromTRUIndex(Int_t itru, Int_t iphitru, Int_t ietatru, Int_t &ietaSM, Int_t &iphiSM) const ; // Tranforms Eta-Phi Cell index in TRU into Eta-Phi index in Super Module @@ -50,8 +54,15 @@ public: void GetGlobal(const TVector3 &vloc, TVector3 &vglob, int ind) const; void GetGlobal(Int_t absId, Double_t glob[3]) const; void GetGlobal(Int_t absId, TVector3 &vglob) const; - // for a given tower index it returns eta and phi of center of that tower. - void EtaPhiFromIndex(Int_t absId,Float_t &eta,Float_t &phi) const; + // for a given tower index absId returns eta and phi of gravity center of tower. + void EtaPhiFromIndex(Int_t absId, Double_t &eta, Double_t &phi) const; + void EtaPhiFromIndex(Int_t absId, Float_t &eta, Float_t &phi) const; + // + Bool_t GetPhiBoundariesOfSM (Int_t nSupMod, Double_t &phiMin, Double_t &phiMax) const; + Bool_t GetPhiBoundariesOfSMGap(Int_t nPhiSec, Double_t &phiMin, Double_t &phiMax) const; + Bool_t SuperModuleNumberFromEtaPhi(Double_t eta, Double_t phi, Int_t &nSupMod) const; + + Bool_t GetAbsCellIdFromEtaPhi(Double_t eta,Double_t phi, Int_t &absId) const; // virtual void GetGlobal(const AliEMCALRecPoint *rp, TVector3 &vglob) const; @@ -143,7 +154,17 @@ public: void GetModulePhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower, Int_t &iphim, Int_t &ietam) const; void GetCellPhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower, Int_t nIphi, Int_t nIeta, Int_t &iphi, Int_t &ieta) const ; - Int_t GetSuperModuleNumber(Int_t absId) const; + Int_t GetSuperModuleNumber(Int_t absId) const; + Int_t GetNumberOfModuleInPhiDirection(Int_t nSupMod) const + { + // inline function + if(fKey110DEG == 1 && nSupMod>=10) return fNPhi/2; + else return fNPhi; + } + // From cell indexes to abs cell id + void GetModuleIndexesFromCellIndexesInSModule(Int_t nSupMod, Int_t iphi, Int_t ieta, + Int_t &iphim, Int_t &ietam, Int_t &nTower) const; + Int_t GetAbsCellIdFromCellIndexes(Int_t nSupMod, Int_t iphi, Int_t ieta) const; // Methods for AliEMCALRecPoint - Feb 19, 2006 Bool_t RelPosCellInSModule(Int_t absId, Double_t &xr, Double_t &yr, Double_t &zr) const; Bool_t RelPosCellInSModule(Int_t absId, Double_t loc[3]) const; @@ -166,7 +187,8 @@ public: void SetSampling(Float_t samp) { fSampling = samp; printf("SetSampling: Sampling factor set to %f", fSampling) ; } Int_t GetNCellsInSupMod() const {return fNCellsInSupMod;} - Int_t GetNCellsInTower() const {return fNCellsInTower; } + Int_t GetNCellsInTower() const {return fNCellsInTower; } + Int_t GetKey110DEG() const {return fKey110DEG;} AliEMCALGeometry(); // default ctor only for internal usage (singleton) @@ -232,6 +254,9 @@ private: Float_t f2Trd1Dx2; // 2*dx2 for TRD1 Float_t fPhiGapForSM; // Gap betweeen supermodules in phi direction Int_t fKey110DEG; // for calculation abs cell id; 19-oct-05 + TArrayD fPhiBoundariesOfSM; // phi boundaries of SM in rad; size is fNumberOfSuperModules; + TArrayD fPhiCentersOfSM; // phi of centers of SMl size is fNumberOfSuperModules/2 + Float_t fEtaMaxOfTRD1; // max eta in case of TRD1 geometry (see AliEMCALShishKebabTrd1Module) // TRD2 options - 27-jan-07 Float_t fTrd2AngleY; // angle in y-z plane (in degree) Float_t f2Trd2Dy2; // 2*dy2 for TRD2 @@ -240,9 +265,12 @@ private: Float_t fTubsR; // radius of tubs Float_t fTubsTurnAngle; // turn angle of tubs in degree // Local Coordinates of SM - TArrayD fEtaCentersOfCells; // size fNEta*fNETAdiv (for TRD1 only) (eta or z in SM) - TArrayD fXCentersOfCells; // size fNEta*fNETAdiv (for TRD1 only) ( x in SM) - TArrayD fPhiCentersOfCells; // size fNPhi*fNPHIdiv (for TRD1 only) (phi or y in SM) + TArrayD fCentersOfCellsEtaDir; // size fNEta*fNETAdiv (for TRD1 only) (eta or z in SM, in cm) + TArrayD fCentersOfCellsXDir; // size fNEta*fNETAdiv (for TRD1 only) ( x in SM, in cm) + TArrayD fCentersOfCellsPhiDir; // size fNPhi*fNPHIdiv (for TRD1 only) (phi or y in SM, in cm) + // + TArrayD fEtaCentersOfCells; // [fNEta*fNETAdiv*fNPhi*fNPHIdiv], positive direction (eta>0); eta depend from phi position; + TArrayD fPhiCentersOfCells; // [fNPhi*fNPHIdiv] from center of SM (-10. < phi < +10.) // Move from AliEMCALv0 - Feb 19, 2006 TList *fShishKebabTrd1Modules; //! list of modules // Local coordinates of SM for TRD1 @@ -252,7 +280,7 @@ private: char *fAdditionalOpts[4]; //! some additional options for the geometry type and name int fNAdditionalOpts; //! size of additional options parameter - ClassDef(AliEMCALGeometry, 10) // EMCAL geometry class + ClassDef(AliEMCALGeometry, 11) // EMCAL geometry class }; #endif // AliEMCALGEOMETRY_H diff --git a/EMCAL/AliEMCALRecPoint.cxx b/EMCAL/AliEMCALRecPoint.cxx index cb399124ff9..2ceb877ea4b 100644 --- a/EMCAL/AliEMCALRecPoint.cxx +++ b/EMCAL/AliEMCALRecPoint.cxx @@ -406,28 +406,53 @@ void AliEMCALRecPoint::EvalAll(Float_t logWeight,TClonesArray * digits) void AliEMCALRecPoint::EvalDispersion(Float_t logWeight, TClonesArray * digits) { // Calculates the dispersion of the shower at the origin of the RecPoint + // in cell units - Nov 16,2006 - Double_t d = 0., wtot = 0., w = 0., xyzi[3], diff=0.; + Double_t d = 0., wtot = 0., w = 0.; Int_t iDigit=0, nstat=0, i=0; AliEMCALDigit * digit ; - - // Calculates the centre of gravity in the local EMCAL-module coordinates - if (!fLocPos.Mag()) - EvalLocalPosition(logWeight, digits) ; - // Calculates the dispersion in coordinates + // Calculates the dispersion in cell units + Double_t etai, phii, etaMean=0.0, phiMean=0.0; + int nSupMod=0, nTower=0, nIphi=0, nIeta=0; + int iphi=0, ieta=0; + // Calculate mean values for(iDigit=0; iDigit < fMulDigit; iDigit++) { digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit]) ; - fGeomPtr->RelPosCellInSModule(digit->GetId(), xyzi[0], xyzi[1], xyzi[2]); - w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; + if (fAmp>0 && fEnergyList[iDigit]>0) { + fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nTower,nIphi,nIeta); + fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta); + etai=(Double_t)ieta; + phii=(Double_t)iphi; + w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; + + if(w>0.0) { + phiMean += phii*w; + etaMean += etai*w; + wtot += w; + } + } + } + if (wtot>0) { + phiMean /= wtot ; + etaMean /= wtot ; + } else AliError(Form("Wrong weight %f\n", wtot)); - if(w>0.0) { - wtot += w; - nstat++; - for(i=0; i<3; i++ ) { - diff = xyzi[i] - double(fLocPos[i]); - d += w * diff*diff; + // Calculate dispersion + for(iDigit=0; iDigit < fMulDigit; iDigit++) { + digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit]) ; + + if (fAmp>0 && fEnergyList[iDigit]>0) { + fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nTower,nIphi,nIeta); + fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta); + etai=(Double_t)ieta; + phii=(Double_t)iphi; + w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; + + if(w>0.0) { + nstat++; + d += w*((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean)); } } } @@ -505,9 +530,10 @@ void AliEMCALRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits) // i.e. within a radius rad = fCoreEnergy around the center. Beyond this radius // in accordance with shower profile the energy deposition // should be less than 2% + // Unfinished - Nov 15,2006 + // Distance is calculate in (phi,eta) units AliEMCALDigit * digit ; - const Float_t kDeg2Rad = TMath::DegToRad() ; Int_t iDigit; @@ -515,16 +541,16 @@ void AliEMCALRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits) EvalLocalPosition(logWeight, digits); } + Double_t phiPoint = fLocPos.Phi(), etaPoint = fLocPos.Eta(); + Double_t eta, phi, distance; for(iDigit=0; iDigit < fMulDigit; iDigit++) { digit = (AliEMCALDigit *) ( digits->At(fDigitsList[iDigit]) ) ; - Float_t ieta = 0.0; - Float_t iphi = 0.0; - //fGeomPtr->PosInAlice(digit->GetId(), ieta, iphi); - fGeomPtr->EtaPhiFromIndex(digit->GetId(),ieta, iphi) ; - iphi = iphi * kDeg2Rad; + eta = phi = 0.0; + fGeomPtr->EtaPhiFromIndex(digit->GetId(),eta, phi) ; + phi = phi * TMath::DegToRad(); - Float_t distance = TMath::Sqrt((ieta-fLocPos.X())*(ieta-fLocPos.X())+(iphi-fLocPos.Y())*(iphi-fLocPos.Y())) ; + distance = TMath::Sqrt((eta-etaPoint)*(eta-etaPoint)+(phi-phiPoint)*(phi-phiPoint)); if(distance < fCoreRadius) fCoreEnergy += fEnergyList[iDigit] ; } @@ -534,40 +560,44 @@ void AliEMCALRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits) void AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits) { // Calculates the axis of the shower ellipsoid in eta and phi + // in cell units - Double_t wtot = 0. ; + static TString gn(fGeomPtr->GetName()); + + Double_t wtot = 0.; Double_t x = 0.; Double_t z = 0.; Double_t dxx = 0.; Double_t dzz = 0.; Double_t dxz = 0.; - const Float_t kDeg2Rad = TMath::DegToRad(); - AliEMCALDigit * digit ; + AliEMCALDigit * digit = 0; - TString gn(fGeomPtr->GetName()); - - Int_t iDigit; - - for(iDigit=0; iDigitAt(fDigitsList[iDigit]) ; - Float_t etai = 0. ; - Float_t phii = 0. ; - if(gn.Contains("SHISH")) { // have to be change - Feb 28, 2006 - //copied for shish-kebab geometry, ieta,iphi is cast as float as eta,phi conversion - // for this geometry does not exist - int nSupMod=0, nTower=0, nIphi=0, nIeta=0; - int iphi=0, ieta=0; + etai = phii = 0.; + if(gn.Contains("SHISH")) { + // Nov 15,2006 - use cell numbers as coordinates + // Copied for shish-kebab geometry, ieta,iphi is cast as double as eta,phi + // We can use the eta,phi(or coordinates) of cell + nSupMod = nTower = nIphi = nIeta = iphi = ieta = 0; + fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nTower,nIphi,nIeta); fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta); - etai=(Float_t)ieta; - phii=(Float_t)iphi; - } else { + etai=(Double_t)ieta; + phii=(Double_t)iphi; + } else { // fGeomPtr->EtaPhiFromIndex(digit->GetId(), etai, phii); - phii = phii * kDeg2Rad; + phii = phii * TMath::DegToRad(); } - Double_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; + w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; + // fAmp summed amplitude of digits, i.e. energy of recpoint + // Gives smaller value of lambda than log weight + // w = fEnergyList[iDigit] / fAmp; // Nov 16, 2006 - try just energy dxx += w * etai * etai ; x += w * etai ; @@ -967,3 +997,11 @@ void AliEMCALRecPoint::Print(Option_t *) const message += " Stored at position %d" ; Info("Print", message.Data(), fClusterType, fMulDigit, fAmp, fCoreEnergy, fCoreRadius, fMulTrack, GetIndexInList() ) ; } + +Double_t AliEMCALRecPoint::GetPointEnergy() const +{ + static double e; + e=0.0; + for(int ic=0; icGetTrd1Angle())*TMath::DegToRad(); fgtanBetta = TMath::Tan(fgangle/2.); fgr = (Double_t)fgGeometry->GetIPDistance(); + if(!sn.Contains("TRD2")) fgr += fgGeometry->GetSteelFrontThickness(); + fga2 = Double_t(fgGeometry->Get2Trd1Dx2()); //PH PrintShish(0); return kTRUE; @@ -215,13 +236,19 @@ void AliEMCALShishKebabTrd1Module::PrintShish(int pri) const printf(" %i |%s| theta %f : fOK.Phi = %f(%5.2f)\n", GetUniqueID(), GetName(), fTheta, fOK.Phi(),fOK.Phi()*TMath::RadToDeg()); printf(" A %f B %f | fThetaA %7.6f(%5.2f)\n", fA,fB, fThetaA,fThetaA*TMath::RadToDeg()); - printf(" fOK : X %9.4f: Y %9.4f \n", fOK.X(), fOK.Y()); - printf(" fOK1 : X %9.4f: Y %9.4f (local, ieta=2)\n", fOK1.X(), fOK1.Y()); - printf(" fOK2 : X %9.4f: Y %9.4f (local, ieta=1)\n\n", fOK2.X(), fOK2.Y()); + printf(" fOK : X %9.4f: Y %9.4f : eta %5.3f\n", fOK.X(), fOK.Y(), GetEtaOfCenterOfModule()); + printf(" fOK1 : X %9.4f: Y %9.4f : (local, ieta=2)\n", fOK1.X(), fOK1.Y()); + printf(" fOK2 : X %9.4f: Y %9.4f : (local, ieta=1)\n\n", fOK2.X(), fOK2.Y()); printf(" fOB : X %9.4f: Y %9.4f \n", fOB.X(), fOB.Y()); printf(" fOB1 : X %9.4f: Y %9.4f (local, ieta=2)\n", fOB1.X(), fOB1.Y()); printf(" fOB2 : X %9.4f: Y %9.4f (local, ieta=1)\n", fOB2.X(), fOB2.Y()); + // 3X3 + printf(" 3X3 \n"); + for(int ieta=0; ieta<3; ieta++) { + printf(" fOK3X3[%i] : X %9.4f: Y %9.4f (local) \n", ieta, fOK3X3[ieta].X(), fOK3X3[ieta].Y()); + } // fOK.Dump(); + GetMaxEtaOfModule(pri); } } } @@ -237,3 +264,26 @@ Double_t AliEMCALShishKebabTrd1Module::GetEtaOfCenterOfModule() const { return -TMath::Log(TMath::Tan(fOK.Phi()/2.)); } + +//_____________________________________________________________________________ +Double_t AliEMCALShishKebabTrd1Module::GetMaxEtaOfModule(int pri) const +{ + // Right bottom point of module + Double_t xBottom = (fgr - fB) / fA; + Double_t thetaBottom = TMath::ATan2(fgr, xBottom); + Double_t etaBottom = ThetaToEta(thetaBottom); + // Right top point of module + Double_t l = fgb / TMath::Cos(fgangle/2.); // length of lateral module side + Double_t xTop = xBottom + l*TMath::Cos(TMath::ATan(fA)); + Double_t yTop = fA*xTop + fB; + Double_t thetaTop = TMath::ATan2(yTop, xTop); + Double_t etaTop = ThetaToEta(thetaTop); + + if(pri) { + printf(" Right bottom point of module : eta %5.4f : theta %6.4f (%6.2f) \n", + etaBottom, thetaBottom, thetaBottom * TMath::RadToDeg()); + printf(" Right top point of module : eta %5.4f : theta %6.4f (%6.2f) \n", + etaTop, thetaTop, thetaTop * TMath::RadToDeg()); + } + return etaBottom>etaTop ? etaBottom : etaTop; +} diff --git a/EMCAL/AliEMCALShishKebabTrd1Module.h b/EMCAL/AliEMCALShishKebabTrd1Module.h index 8ee884c536e..0e89a47ffa5 100644 --- a/EMCAL/AliEMCALShishKebabTrd1Module.h +++ b/EMCAL/AliEMCALShishKebabTrd1Module.h @@ -20,6 +20,7 @@ class AliEMCALShishKebabTrd1Module : public TNamed { AliEMCALShishKebabTrd1Module(Double_t theta=0.0, AliEMCALGeometry *g=0); AliEMCALShishKebabTrd1Module(AliEMCALShishKebabTrd1Module &leftNeighbor); void Init(Double_t A, Double_t B); + void DefineAllStaff(); AliEMCALShishKebabTrd1Module(const AliEMCALShishKebabTrd1Module& mod); AliEMCALShishKebabTrd1Module & operator = (const AliEMCALShishKebabTrd1Module& /*rvalue*/) { @@ -40,17 +41,25 @@ class AliEMCALShishKebabTrd1Module : public TNamed { Double_t GetPosXfromR() const {return fOK.Y() - fgr;} Double_t GetA() const {return fA;} Double_t GetB() const {return fB;} + Double_t GetRadius() const {return fgr;} // Additional offline staff // ieta=0 or 1 - Jun 02, 2006 TVector2& GetCenterOfCellInLocalCoordinateofSM(Int_t ieta) { if(ieta<=0) return fOK2; else return fOK1;} - void GetCenterOfCellInLocalCoordinateofSM(Int_t ieta, Double_t &xr, Double_t &zr) + void GetCenterOfCellInLocalCoordinateofSM(Int_t ieta, Double_t &xr, Double_t &zr) const { if(ieta<=0) {xr = fOK2.Y(); zr = fOK2.X(); } else {xr = fOK1.Y(); zr = fOK1.X(); } } + void GetCenterOfCellInLocalCoordinateofSM_3X3(Int_t ieta, Double_t &xr, Double_t &zr) const + { // 3X3 case - Nov 9,2006 + ieta = ieta<0? ieta=0 : ieta; // check index + ieta = ieta>2? ieta=2 : ieta; + xr = fOK3X3[2-ieta].Y(); + zr = fOK3X3[2-ieta].X(); + } // 15-may-06 TVector2& GetCenterOfModuleFace() {return fOB;} TVector2& GetCenterOfModuleFace(Int_t ieta) { @@ -62,8 +71,11 @@ class AliEMCALShishKebabTrd1Module : public TNamed { Double_t Getb() const {return fgb;} // service methods void PrintShish(Int_t pri=1) const; // *MENU* - Double_t GetThetaInDegree() const; + Double_t GetThetaInDegree() const; Double_t GetEtaOfCenterOfModule() const; + Double_t GetMaxEtaOfModule(int pri=0) const; + static Double_t ThetaToEta(Double_t theta) + {return -TMath::Log(TMath::Tan(theta/2.));} protected: // geometry info @@ -79,7 +91,7 @@ class AliEMCALShishKebabTrd1Module : public TNamed { Double_t fA; // parameters of right line : y = A*z + B Double_t fB; // system where zero point is IP. Double_t fThetaA; // angle coresponding fA - for convinience - Double_t fTheta; // theta angle of perependicular to SK module + Double_t fTheta; // theta angle of perpendicular to SK module // position of towers(cells) with differents ieta (1 or 2) in local coordinate of SM // Nov 04,2004; Feb 19,2006 TVector2 fOK1; // ieta=1 @@ -88,9 +100,10 @@ class AliEMCALShishKebabTrd1Module : public TNamed { TVector2 fOB; // module TVector2 fOB1; // ieta=1 TVector2 fOB2; // ieta=0 - + // 3X3 case - Nov 9,2006 + TVector2 fOK3X3[3]; // public: - ClassDef(AliEMCALShishKebabTrd1Module,0) // TRD1 Shish-Kebab module + ClassDef(AliEMCALShishKebabTrd1Module,1) // TRD1 Shish-Kebab module }; #endif diff --git a/EMCAL/AliEMCALv0.cxx b/EMCAL/AliEMCALv0.cxx index 366d63bba59..6a9efd9b4af 100644 --- a/EMCAL/AliEMCALv0.cxx +++ b/EMCAL/AliEMCALv0.cxx @@ -403,10 +403,10 @@ void AliEMCALv0::CreateShishKebabGeometry() CreateEmod("SMOD","EMOD"); // 18-may-05 - if(gn.Contains("110DEG")) CreateEmod("SM10","EMOD"); // 12-oct-05 + if(g->GetKey110DEG()) CreateEmod("SM10","EMOD"); // Nov 1,2006 // Sensitive SC (2x2 tiles) - double parSCM0[5], *dummy = 0, parTRAP[11]; + double parSCM0[5]={0,0,0,0}, *dummy = 0, parTRAP[11]; Double_t trd1Angle = g->GetTrd1Angle()*TMath::DegToRad(), tanTmp = TMath::Tan(trd1Angle/2.); if(!gn.Contains("TRD")) { // standard module par[0] = (g->GetECPbRadThick()+g->GetECScintThick())*g->GetNECLayers()/2.; @@ -443,7 +443,8 @@ void AliEMCALv0::CreateShishKebabGeometry() gMC->Gsvolu("SCM0", "TRD1", fIdTmedArr[kIdAIR], parSCM0, 4); gMC->Gspos("SCM0", 1, "EMOD", 0., 0., dzTmp/2., 0, "ONLY") ; } else { // before MAY 2005 - double wallThickness = g->GetPhiModuleSize()/2. - g->GetPhiTileSize(); // Need check + // double wallThickness = g->GetPhiModuleSize()/2. - g->GetPhiTileSize(); // Need check + double wallThickness = g->GetPhiModuleSize()/g->GetNPHIdiv() - g->GetPhiTileSize(); for(int i=0; i<3; i++) parSCM0[i] = fParEMOD[i] - wallThickness; parSCM0[3] = fParEMOD[3]; gMC->Gsvolu("SCM0", "TRD1", fIdTmedArr[kIdAIR], parSCM0, 4); @@ -502,13 +503,16 @@ void AliEMCALv0::CreateShishKebabGeometry() AliDebug(2,Form(" Number of Pb tiles in SCMX %i \n", nr)); } } else if(g->GetNPHIdiv()==3 && g->GetNETAdiv()==3) { - Trd1Tower3X3(parSCM0); + printf(" before AliEMCALv0::Trd1Tower3X3() : parSCM0"); + for(int i=0; i<4; i++) printf(" %7.4f ", parSCM0[i]); + printf("\n"); + Trd1Tower3X3(parSCM0); // Stay here } else if(g->GetNPHIdiv()==4 && g->GetNETAdiv()==4) { Trd1Tower4X4(); } } else if(gn.Contains("TRD2")) { // TRD2 - 14-jan-05 // Scm0InTrd2(g, fParEMOD, parSCM0); // First dessin - PbmoInTrd2(g, fParEMOD, parSCM0); // Second dessin + PbmoInTrd2(g, fParEMOD, parSCM0); // Second design } } @@ -532,7 +536,7 @@ void AliEMCALv0::CreateSmod(const char* mother) int nphism = g->GetNumberOfSuperModules()/2; // 20-may-05 if(nphism>0) { dphi = (g->GetArm1PhiMax() - g->GetArm1PhiMin())/nphism; - // if(gn.Contains("110DEG")) dphi = (g->GetArm1PhiMax() - g->GetArm1PhiMin())/(nphism-1); + // if(g->GetKey110DEG()) dphi = (g->GetArm1PhiMax() - g->GetArm1PhiMin())/(nphism-1); rpos = (g->GetEnvelop(0) + g->GetEnvelop(1))/2.; AliDebug(2,Form(" rpos %8.2f : dphi %6.1f degree \n", rpos, dphi)); } @@ -580,7 +584,7 @@ void AliEMCALv0::CreateSmod(const char* mother) fIdTmedArr[kIdAIR], par[0],par[1],par[2])); fSmodPar0 = par[0]; fSmodPar2 = par[2]; - if(gn.Contains("110DEG")) { // 12-oct-05 + if(g->GetKey110DEG()) { // 12-oct-05 par1C = par[1]; par[1] /= 2.; gMC->Gsvolu("SM10", "BOX", fIdTmedArr[kIdAIR], par, 3); @@ -653,7 +657,7 @@ void AliEMCALv0::CreateSmod(const char* mother) nr++; } else { // TRD1 TString smName("SMOD"); // 12-oct-05 - if(i==5 && gn.Contains("110DEG")) { + if(i==5 && g->GetKey110DEG()) { smName = "SM10"; nrsmod = nr; nr = 0; @@ -666,7 +670,7 @@ void AliEMCALv0::CreateSmod(const char* mother) xpos = rpos * TMath::Cos(phiRad); ypos = rpos * TMath::Sin(phiRad); zpos = fSmodPar2; // 21-sep-04 - if(i==5 && gn.Contains("110DEG")) { + if(i==5 && g->GetKey110DEG()) { xpos += (par1C/2. * TMath::Sin(phiRad)); ypos -= (par1C/2. * TMath::Cos(phiRad)); } @@ -780,7 +784,7 @@ void AliEMCALv0::CreateEmod(const char* mother, const char* child) zpos = mod->GetPosZ() - fSmodPar2; int iyMax = g->GetNPhi(); - if(strcmp(mother,"SMOD") && gn.Contains("110DEG")) { + if(strcmp(mother,"SMOD") && g->GetKey110DEG()) { iyMax /= 2; } for(int iy=0; iyGetName()), scmx; @@ -880,16 +888,13 @@ void AliEMCALv0::Trd1Tower3X3(const double parSCM0[4]) double ndiv=3., xpos=0.0; // should be defined once gMC->Gsvolu("PBTI", "BOX", fIdTmedArr[kIdPB], dummy, 0); - if(gn.Contains("TEST")==0) { // one name for all trapesoid - scmx = "SCMX"; - gMC->Gsvolu(scmx.Data(), "TRAP", fIdTmedArr[kIdSC], dummy, 0); - } - - for(int ix=1; ix<=3; ix++) { // 3X3 + scmx = "SCX"; // Nov 10,2006 // ix=1 parTRAP[0] = dz; - parTRAP[1] = TMath::ATan2((dx2-dx1)/2.,2.*dz)*TMath::RadToDeg(); // theta + double xCentBot = 2.*dx1/3.; + double xCentTop = 2.*(dx2/4. + dx1/12.); + parTRAP[1] = TMath::ATan2((xCentTop-xCentBot),2.*dz)*TMath::RadToDeg(); // theta parTRAP[2] = 0.; // phi // bottom parTRAP[3] = dy/ndiv; // H1 @@ -898,10 +903,10 @@ void AliEMCALv0::Trd1Tower3X3(const double parSCM0[4]) parTRAP[6] = 0.0; // ALP1 // top parTRAP[7] = dy/ndiv; // H2 - parTRAP[8] = dx2/ndiv; // BL2 + parTRAP[8] = dx2/2 - dx1/6.;// BL2 parTRAP[9] = parTRAP[8]; // TL2 parTRAP[10]= 0.0; // ALP2 - xpos = +(dx1+dx2)/3.; // 6 or 3 + xpos = (xCentBot+xCentTop)/2.; if (ix==3) { parTRAP[1] = -parTRAP[1]; @@ -914,13 +919,11 @@ void AliEMCALv0::Trd1Tower3X3(const double parSCM0[4]) } AliDebug(2,Form(" ** TRAP ** xpos %9.3f\n", xpos)); for(int i=0; i<11; i++) AliDebug(2,Form(" par[%2.2i] %9.4f\n", i, parTRAP[i])); - if(gn.Contains("TEST")){ - scmx = "SCX"; scmx += ix; - gMC->Gsvolu(scmx.Data(), "TRAP", fIdTmedArr[kIdSC], parTRAP, 11); - gMC->Gspos(scmx.Data(), 1, "SCMY", xpos, 0.0, 0.0, 0, "ONLY") ; - } else { - gMC->Gsposp(scmx.Data(), ix, "SCMY", xpos, 0.0, 0.0, 0, "ONLY", parTRAP, 11) ; - } + + scmx += ix; + gMC->Gsvolu(scmx.Data(), "TRAP", fIdTmedArr[kIdSC], parTRAP, 11); + gMC->Gspos(scmx.Data(), 1, "SCMY", xpos, 0.0, 0.0, 0, "ONLY") ; + PbInTrap(parTRAP, scmx); } AliDebug(2,Form("Trd1Tower3X3()", "Ver. 1.0 : was tested.")); @@ -1244,9 +1247,7 @@ void AliEMCALv0::AddAlignableVolumes() const AliFatal("Unable to set alignable entry!!"); } - TString gn( ((AliEMCALGeometry*)GetGeometry())->GetName() ); - gn.ToUpper(); - if(gn.Contains("110DEG")) { + if(GetGeometry()->GetKey110DEG()) { TString vpstr2 = "ALIC_1/XEN1_1/SM10_"; TString snstr2 = "EMCAL/HalfSupermodule"; for (Int_t smodnum=0; smodnum < 2; smodnum++) { diff --git a/EMCAL/AliEMCALv2.cxx b/EMCAL/AliEMCALv2.cxx index 40c09ce5169..1e8715b10f2 100644 --- a/EMCAL/AliEMCALv2.cxx +++ b/EMCAL/AliEMCALv2.cxx @@ -77,7 +77,7 @@ AliEMCALv2::AliEMCALv2(const char *name, const char *title) // if (gDebug>0){ if (1){ TH1::AddDirectory(0); - fHDe = new TH1F("fHDe","De in EMCAL", 1000, 0., 1.); + fHDe = new TH1F("fHDe","De in EMCAL", 1000, 0., 10.); fHNhits = new TH1F("fHNhits","#hits in EMCAL", 2001, -0.5, 2000.5); fHistograms->Add(fHDe); fHistograms->Add(fHNhits); @@ -143,7 +143,6 @@ void AliEMCALv2::StepManager(void){ static int supModuleNumber, moduleNumber, yNumber, xNumber, absid; static int keyGeom=1; static char *vn = "SCMX"; // Apr 13, 2006 - only TRD1 case now - // static char *vn = "SX"; // 15-mar-05 static int nSMOP[7]={1,3,5,7,9,11}; // 30-mar-05 static int nSMON[7]={2,4,6,8,10,12}; static Float_t depositedEnergy=0.0; @@ -161,7 +160,7 @@ void AliEMCALv2::StepManager(void){ Int_t tracknumber = gAlice->GetMCApp()->GetCurrentTrackNumber(); curVolName = gMC->CurrentVolName(); - if(curVolName.Contains(vn)) { // We are in a scintillator layer + if(curVolName.Contains(vn) || curVolName.Contains("SCX")) { // We are in a scintillator layer; SCX for 3X3 // printf(" keyGeom %i : Sensetive volume %s (%s) \n", keyGeom, curVolName.Data(), vn); if( ((depositedEnergy = gMC->Edep()) > 0.) && (gMC->TrackTime() < fTimeCut)){// Track is inside a scintillator and deposits some energy @@ -221,6 +220,14 @@ void AliEMCALv2::StepManager(void){ gMC->CurrentVolOffID(1, yNumber); gMC->CurrentVolOffID(0, xNumber); // really x number now if(strcmp(gMC->CurrentVolOffName(4),"SM10")==0) supModuleNumber += 10; // 13-oct-05 + // Nov 10,2006 + xNumber == 0; + if (strcmp(gMC->CurrentVolOffName(0),"SCX1")==0) xNumber=1; + else if(strcmp(gMC->CurrentVolOffName(0),"SCX2")==0) xNumber=2; + else if(strcmp(gMC->CurrentVolOffName(0),"SCX3")==0) xNumber=3; + if(xNumber==0) { + Fatal("StepManager()", "Wrong name SCX : %s ", gMC->CurrentVolOffName(0)) ; + } } else { gMC->CurrentVolOffID(5, supModuleNumber); gMC->CurrentVolOffID(4, moduleNumber); @@ -232,7 +239,11 @@ void AliEMCALv2::StepManager(void){ } absid = fGeometry->GetAbsCellId(supModuleNumber-1, moduleNumber-1, yNumber-1, xNumber-1); - if (absid < 0) Fatal("StepManager()", "Wrong id ") ; + if (absid < 0) { + printf(" supModuleNumber %i : moduleNumber %i : yNumber %i : xNumber %i \n", + supModuleNumber, moduleNumber, yNumber, xNumber); + Fatal("StepManager()", "Wrong id : %i ", absid) ; + } Float_t lightYield = depositedEnergy ; // Apply Birk's law (copied from G3BIRK) -- 2.43.0