New methods in geometry, update geometry for 3x3 case
authorpavlinov <pavlinov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 21 Nov 2006 22:36:04 +0000 (22:36 +0000)
committerpavlinov <pavlinov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 21 Nov 2006 22:36:04 +0000 (22:36 +0000)
EMCAL/AliEMCALClusterizerv1.cxx
EMCAL/AliEMCALClusterizerv1.h
EMCAL/AliEMCALGeometry.cxx
EMCAL/AliEMCALGeometry.h
EMCAL/AliEMCALRecPoint.cxx
EMCAL/AliEMCALRecPoint.h
EMCAL/AliEMCALShishKebabTrd1Module.cxx
EMCAL/AliEMCALShishKebabTrd1Module.h
EMCAL/AliEMCALv0.cxx
EMCAL/AliEMCALv2.cxx

index 700aacc..7591979 100644 (file)
@@ -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<TClonesArray*>(digits->Clone());
 
@@ -590,7 +594,7 @@ void AliEMCALClusterizerv1::MakeClusters(char* option)
     int ineb=0;
     nextdigit.Reset();
 
-    // PseudoClusterization starts    
+    // PseudoClusterization starts  
     while ( (digit = dynamic_cast<AliEMCALDigit *>(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<TClonesArray*>(digits->Clone());
+  TClonesArray * digitsC =  dynamic_cast<TClonesArray*>(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<AliEMCALDigit *>(nextdigit())) ) { // clean up digits
+  double e = 0.0, ehs = 0.0;
+  while ( (digit = dynamic_cast<AliEMCALDigit *>(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<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
+  while ( (digit = dynamic_cast<AliEMCALDigit *>(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<AliEMCALLoader*>(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<AliEMCALRecPoint * >(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; iprimary<nprimaries; iprimary++) {
-       printf("%d ", primaries[iprimary] ) ; 
-      } 
+      if(strstr(option,"deb")){ 
+        for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
+         printf("%d ", primaries[iprimary] ) ; 
+        }
+      }
     }
 
       fMaxE->Fill(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);
 }
index bf33d2a..ee2f4f2 100644 (file)
@@ -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
 
index c1f9a03..9c16ac2 100644 (file)
 // -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 <assert.h>
-#include "Riostream.h"
+#include <Riostream.h>
 
 #include <TMath.h>
 #include <TVector3.h>
-             //#include <TArrayD.h>
 #include <TObjArray.h>
 #include <TGeoManager.h>
 #include <TGeoNode.h>
@@ -43,6 +46,7 @@
 #include <TMatrixD.h>
 #include <TObjString.h>
 #include <TClonesArray.h>
+#include <TBrowser.h>
 
 // -- 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<phi<190)
-              fArm1PhiMax = 200.0; // for XEN1 and turn angle of super modules
+              fArm1PhiMax = 200.0;       // for XEN1 and turn angle of super modules
+           }
+            if(fGeoName.Contains("FINAL")) {
+              fPhiModuleSize = 12.26 - fPhiGapForSM / Float_t(fNPhi); // first assumption
+           } else if(fGeoName.Contains("CURRENT")) {
+              fECScintThick      = 0.176; // 10% of weight reduction
+              fECPbRadThickness  = 0.144; //
+              fLateralSteelStrip = 0.015; // 0.015cm  = 0.15mm (Oct 30, from Fred)
+              fPhiModuleSize     = 12.00;
+              fPhiGapForSM       = (12.26 - fPhiModuleSize)*fNPhi; // have to check
            }
-            fPhiModuleSize = 12.26 - fPhiGapForSM / Float_t(fNPhi); // first assumption
             fEtaModuleSize = fPhiModuleSize;
             if(fGeoName.Contains("HUGE")) fNECLayers *= 3; // 28-oct-05 for analysing leakage    
           }
@@ -285,14 +316,14 @@ void AliEMCALGeometry::Init(void){
     CheckAdditionalOptions();
     DefineSamplingFraction();
 
-    fPhiTileSize = fPhiModuleSize/2. - fLateralSteelStrip; // 13-may-05 
-    fEtaTileSize = fEtaModuleSize/2. - fLateralSteelStrip; // 13-may-05 
+    fPhiTileSize = fPhiModuleSize/double(fNPHIdiv) - fLateralSteelStrip; // 13-may-05 
+    fEtaTileSize = fEtaModuleSize/double(fNETAdiv) - fLateralSteelStrip; // 13-may-05 
 
     // constant for transition absid <--> 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<phi<190)\n");
-      }
-    }
-    printf("Granularity: %d in eta and %d in phi\n", GetNZ(), GetNPhi()) ;
-    printf("Layout: phi = (%7.1f, %7.1f), eta = (%5.2f, %5.2f), IP = %7.2f\n",  
-          GetArm1PhiMin(), GetArm1PhiMax(),GetArm1EtaMin(), GetArm1EtaMax(), GetIPDistance() );
-  }
   //TRU parameters. These parameters values are not the final ones.
   fNTRU    = 3 ;
   fNTRUEta = 3 ;
   fNTRUPhi = 1 ;
+
 }
 
-//______________________________________________________________________
+void AliEMCALGeometry::PrintGeometry()
+{
+  // Separate routine is callable from broswer; Nov 7,2006
+  printf("\nInit: geometry of EMCAL named %s is as follows:\n", fGeoName.Data());
+  printf("Granularity: %d in eta and %d in phi\n", GetNZ(), GetNPhi()) ;
+  printf("Layout: phi = (%7.1f, %7.1f), eta = (%5.2f, %5.2f), IP = %7.2f -> 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<phi<190)\n");
+      printf(" phi SM boundaries \n"); 
+      for(int i=0; i<fPhiBoundariesOfSM.GetSize()/2.; i++) {
+       printf(" %i : %7.5f(%7.2f) -> %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<fCentersOfCellsEtaDir.GetSize(); i++) {
+        printf(" ind %2.2i : z %8.3f : x %8.3f \n", i, 
+              fCentersOfCellsEtaDir.At(i),fCentersOfCellsXDir.At(i));
+        int ind=0; // Nov 21,2006
+        for(Int_t iphi=0; iphi<fCentersOfCellsPhiDir.GetSize(); iphi++) {
+          ind = iphi*fCentersOfCellsEtaDir.GetSize() + i;
+          printf("%6.4f ", fEtaCentersOfCells[ind]);
+          if((iphi+1)%12 == 0) printf("\n");
+       }
+        printf("\n");
+      }
+
+      printf("\n Cells grid in phi directions : size %i\n", fCentersOfCellsPhiDir.GetSize());
+      for(Int_t i=0; i<fCentersOfCellsPhiDir.GetSize(); i++) {
+        double phi=fPhiCentersOfCells.At(i);
+        printf(" ind %2.2i : y %8.3f : phi %7.5f(%6.2f) \n", i, fCentersOfCellsPhiDir.At(i), 
+              phi, phi*TMath::RadToDeg());
+      }
+    }
+  }
+}
+
+void AliEMCALGeometry::PrintCellIndexes(Int_t absId, int pri, char *tit)
+{
+  // Service methods
+  Int_t nSupMod, nTower, nIphi, nIeta;
+  Int_t iphi, ieta;
+  TVector3 vg;
+
+  GetCellIndex(absId,  nSupMod, nTower, nIphi, nIeta);
+  printf(" %s | absId : %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 "<<absId<<" nSupMod "<<nSupMod << " iphi "<<iphi<<" ieta "<<ieta;
-    //    cout<< " xr " << xr << " yr " << yr << " zr " << zr <<endl;
+    yr = fCentersOfCellsPhiDir.At(iphi + phiIndexShift);
   }
+  if(absId==1104 || absId==1105) 
+  AliDebug(2,Form("absId %i nSupMod %i iphi %i ieta %i xr %f yr %f zr %f ",absId,nSupMod,iphi,ieta,xr,yr,zr));
 
   return kTRUE;
 }
@@ -794,15 +930,16 @@ Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, TVector3 &vloc) const
 
 void AliEMCALGeometry::CreateListOfTrd1Modules()
 {
-  //Generate the list of Trd1 modules
-  //which will make up the EMCAL
-  //geometry
+  // Generate the list of Trd1 modules
+  // which will make up the EMCAL
+  // geometry
 
   AliDebug(2,Form(" AliEMCALGeometry::CreateListOfTrd1Modules() started "));
 
   AliEMCALShishKebabTrd1Module *mod=0, *mTmp=0; // current module
   if(fShishKebabTrd1Modules == 0) {
     fShishKebabTrd1Modules = new TList;
+    fShishKebabTrd1Modules->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; it<fNZ; it++) { // array index
+  Double_t xr, zr, theta, phi, eta, r, x,y;
+  TVector3 vglob;
+  Double_t ytCenterModule, ytCenterCell;
+
+  fCentersOfCellsPhiDir.Set(fNPhi*fNPHIdiv);
+  fPhiCentersOfCells.Set(fNPhi*fNPHIdiv);
+
+  Double_t R0 = GetIPDistance() + GetLongModuleSize()/2.;
+  for(Int_t it=0; it<fNPhi; it++) { // cycle on modules
+    ytCenterModule = -fParSM[1] + fPhiModuleSize*(2*it+1)/2;  // center of module
+    for(Int_t ic=0; ic<fNPHIdiv; ic++) { // cycle on cells in module
+      if(fNPHIdiv==2) {
+        ytCenterCell = ytCenterModule + fPhiTileSize *(2*ic-1)/2.;
+      } else if(fNPHIdiv==3){
+        ytCenterCell = ytCenterModule + fPhiTileSize *(ic-1);
+      }
+      fCentersOfCellsPhiDir.AddAt(ytCenterCell,ind);
+      // Define grid on phi direction
+      // Grid is not the same for different eta bin;
+      // Effect is small but is still here
+      phi = TMath::ATan2(ytCenterCell, R0);
+      fPhiCentersOfCells.AddAt(phi, ind);
+
+      AliDebug(2,Form(" ind %2.2i : y %8.3f ", ind, fCentersOfCellsPhiDir.At(ind))); 
+      ind++;
+    }
+  }
+
+  fCentersOfCellsEtaDir.Set(fNZ *fNETAdiv);
+  fCentersOfCellsXDir.Set(fNZ *fNETAdiv);
+  fEtaCentersOfCells.Set(fNZ *fNETAdiv * fNPhi*fNPHIdiv);
+  AliDebug(2,Form(" Cells grid in eta directions : size %i\n", fCentersOfCellsEtaDir.GetSize()));
+  for(Int_t it=0; it<fNZ; it++) {
     AliEMCALShishKebabTrd1Module *trd1 = GetShishKebabModule(it);
     nTower = fNPhi*it;
-    for(Int_t ic=0; ic<fNETAdiv; ic++) { // array index
-      trd1->GetCenterOfCellInLocalCoordinateofSM(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; ic<fNETAdiv; ic++) {
+      if(fNPHIdiv==2) {
+        trd1->GetCenterOfCellInLocalCoordinateofSM(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; iphi<fCentersOfCellsPhiDir.GetSize(); iphi++) {
+        x = xr + trd1->GetRadius();
+        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; i<fEtaCentersOfCells.GetSize(); i++) {
+  for(Int_t i=0; i<fCentersOfCellsEtaDir.GetSize(); i++) {
     AliDebug(2,Form(" ind %2.2i : z %8.3f : x %8.3f", i+1, 
-                    fEtaCentersOfCells.At(i),fXCentersOfCells.At(i)));
+                    fCentersOfCellsEtaDir.At(i),fCentersOfCellsXDir.At(i)));
   }
 
- // define grid for cells in phi(y) direction in local coordinates system of SM
-  //  fPhiCentersOfCells = new TArrayD(fNPhi*fNPHIdiv);
-  fPhiCentersOfCells.Set(fNPhi*fNPHIdiv);
-  AliDebug(2,Form(" Cells grid in phi directions : size %i\n", fPhiCentersOfCells.GetSize()));
-  Int_t ind=0;
-  for(Int_t it=0; it<fNPhi; it++) { // array index
-    Float_t ytLeftCenterModule = -fParSM[1] + fPhiModuleSize*(2*it+1)/2;         // module
-    for(Int_t ic=0; ic<fNPHIdiv; ic++) { // array index
-      Float_t ytLeftCenterCell = ytLeftCenterModule + fPhiTileSize *(2*ic-1)/2.; // tower(cell) 
-      fPhiCentersOfCells.AddAt(ytLeftCenterCell,ind);
-      AliDebug(2,Form(" ind %2.2i : y %8.3f ", ind, fPhiCentersOfCells.At(ind))); 
-      ind++;
-    }
-  }
 }
 
 void  AliEMCALGeometry::GetTransformationForSM()
@@ -963,16 +1142,124 @@ void AliEMCALGeometry::GetGlobal(const AliRecPoint *rp, TVector3 &vglob) const
   GetGlobal(vloc, vglob, nSupMod);
 }
 
-void AliEMCALGeometry::EtaPhiFromIndex(Int_t absId,Float_t &eta,Float_t &phi) const
+void AliEMCALGeometry::EtaPhiFromIndex(Int_t absId,Double_t &eta,Double_t &phi) const
 {
-  // Jun 03, 2006 - version for TRD1
+  // Nov 16, 2006- float to double
+  // version for TRD1 only
   static TVector3 vglob;
   GetGlobal(absId, vglob);
   eta = vglob.Eta();
   phi = vglob.Phi();
 }
 
-AliEMCALShishKebabTrd1Module* AliEMCALGeometry::GetShishKebabModule(Int_t neta=0)
+void AliEMCALGeometry::EtaPhiFromIndex(Int_t absId,Float_t &eta,Float_t &phi) const
+{
+  // Nov 16,2006 - should be discard in future
+  static TVector3 vglob;
+  GetGlobal(absId, vglob);
+  eta = float(vglob.Eta());
+  phi = float(vglob.Phi());
+}
+
+Bool_t AliEMCALGeometry::GetPhiBoundariesOfSM(Int_t nSupMod, Double_t &phiMin, Double_t &phiMax) const
+{
+  // 0<= nSupMod <=11; phi in rad
+  static int i;
+  if(nSupMod<0 || nSupMod >11) 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; i<nphi; i++) {
+      d = TMath::Abs(fPhiCentersOfCells[i] - phiLoc);
+      if(d < dmin) {
+        dmin = d;
+        iphi = i;
+      }
+      //      printf(" i %i : d %f : dmin %f : fPhiCentersOfCells[i] %f \n", i, d, dmin, fPhiCentersOfCells[i]);
+    }
+    // odd SM are turned with respect of even SM - reverse indexes
+    AliDebug(2,Form(" iphi %i : dmin %f (phi %f, phiLoc %f ) ", iphi, dmin, phi, phiLoc));
+    // eta index
+    absEta   = TMath::Abs(eta);
+    etaShift = iphi*fCentersOfCellsEtaDir.GetSize();
+    dmin     = TMath::Abs(fEtaCentersOfCells[etaShift]-absEta);
+    ieta     = 0;
+    for(i=1; i<fCentersOfCellsEtaDir.GetSize(); i++) {
+      d = TMath::Abs(fEtaCentersOfCells[i+etaShift] - absEta);
+      if(d < dmin) {
+        dmin = d;
+        ieta = i;
+      }
+    }
+    AliDebug(2,Form(" ieta %i : dmin %f (eta=%f) : nSupMod %i ", ieta, dmin, eta, nSupMod));
+
+    if(eta<0) iphi = (nphi-1) - iphi;
+    absId = GetAbsCellIdFromCellIndexes(nSupMod, iphi, ieta);
+
+    return kTRUE;
+  }
+  return kFALSE;
+}
+
+AliEMCALShishKebabTrd1Module* AliEMCALGeometry::GetShishKebabModule(Int_t neta)
 {
   //This method was too long to be
   //included in the header file - the
@@ -986,3 +1273,14 @@ AliEMCALShishKebabTrd1Module* AliEMCALGeometry::GetShishKebabModule(Int_t neta=0
   } else trd1 = 0;
   return trd1;
 }
+
+void AliEMCALGeometry::Browse(TBrowser* b)
+{
+  if(fShishKebabTrd1Modules) b->Add(fShishKebabTrd1Modules);
+}
+
+Bool_t AliEMCALGeometry::IsFolder() const
+{
+  if(fShishKebabTrd1Modules) return kTRUE;
+  else                       return kFALSE;
+}
index c86f922..db79b9e 100644 (file)
@@ -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
index cb39912..2ceb877 100644 (file)
@@ -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; iDigit<fMulDigit; iDigit++) {
+  Double_t etai , phii, w; 
+  int nSupMod=0, nTower=0, nIphi=0, nIeta=0;
+  int iphi=0, ieta=0;
+  for(Int_t iDigit=0; iDigit<fMulDigit; iDigit++) {
     digit = (AliEMCALDigit *) digits->At(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; ic<GetMultiplicity(); ic++) e += double(fEnergyList[ic]);
+  return e;
+}
index 64955d8..fad2b82 100644 (file)
@@ -59,6 +59,7 @@ class AliEMCALRecPoint : public AliRecPoint {
   virtual void    GetElipsAxis(Float_t * lambda)const {lambda[0] = fLambda[0]; lambda[1] = fLambda[1];};
   
   Float_t *   GetEnergiesList() const {return fEnergyList ;}       // gets the list of energies making this recpoint
+  Double_t    GetPointEnergy() const;                              // gets point energy (sum of energy list)
   Float_t *   GetTimeList() const {return fTimeList ;}       // gets the list of digit times in this recpoint
   Float_t     GetMaximalEnergy(void) const ;                       // get the highest energy in the cluster
   Int_t       GetMaximumMultiplicity() const {return fMaxDigit ;}  // gets the maximum number of digits allowed
index fc26e45..5847ef1 100644 (file)
 //_________________________________________________________________________
 // Main class for TRD1 geometry of Shish-Kebab case.
 // Author: Aleksei Pavlinov(WSU).
-// Sep 20004.
+// Sep 20004 - Nov 2006
 // See web page with description of Shish-Kebab geometries:
 // http://pdsfweb01.nersc.gov/~pavlinov/ALICE/SHISHKEBAB/RES/shishkebabALICE.html
+// Nov 9,2006 - added cas of 3X3
 //_________________________________________________________________________
 
 #include "AliLog.h"
@@ -53,7 +54,8 @@ AliEMCALShishKebabTrd1Module::AliEMCALShishKebabTrd1Module(Double_t theta, AliEM
     fOK2(),
     fOB(),
     fOB1(),
-    fOB2()
+    fOB2(),
+    fOK3X3()
 { 
   // theta in radians ; first object shold be with theta=pi/2.
   if(fgGeometry==0) {
@@ -79,11 +81,14 @@ AliEMCALShishKebabTrd1Module::AliEMCALShishKebabTrd1Module(AliEMCALShishKebabTrd
     fOK2(),
     fOB(),
     fOB1(),
-    fOB2()
+    fOB2(),
+    fOK3X3()
 { 
   //  printf("** Left Neighbor : %s **\n", leftNeighbor.GetName());
-  TObject::SetUniqueID(leftNeighbor.GetUniqueID()+1);
   fTheta  = leftNeighbor.GetTheta() - fgangle; 
+
+  TObject::SetUniqueID(leftNeighbor.GetUniqueID()+1);
+
   Init(leftNeighbor.GetA(),leftNeighbor.GetB());
 }
 
@@ -99,7 +104,8 @@ AliEMCALShishKebabTrd1Module::AliEMCALShishKebabTrd1Module(const AliEMCALShishKe
     fOK2(mod.fOK2),
     fOB(mod.fOB),
     fOB1(mod.fOB1),
-    fOB2(mod.fOB2)
+    fOB2(mod.fOB2),
+    fOK3X3(mod.fOK3X3)
 {
   //copy ctor
 }
@@ -130,8 +136,13 @@ void AliEMCALShishKebabTrd1Module::Init(Double_t A, Double_t B)
   fA = TMath::Tan(fThetaA); // !!
   fB = yA - fA*xA;
 
+  DefineAllStaff();
+}
+
+void AliEMCALShishKebabTrd1Module::DefineAllStaff()
+{
   DefineName(fTheta);
-  // Centers of module
+  // Centers of module - 2X2 case
   Double_t kk1 = (fga+fga2)/(2.*4.); // kk1=kk2 
 
   Double_t xk1 = fOK.X() - kk1*TMath::Sin(fTheta);
@@ -142,10 +153,24 @@ void AliEMCALShishKebabTrd1Module::Init(Double_t A, Double_t B)
   Double_t yk2 = fOK.Y() - kk1*TMath::Cos(fTheta) - fgr;
   fOK2.Set(xk2,yk2);
 
+  // Centers of module - 3X3 case; Nov 9,2006
+  fOK3X3[1].Set(fOK.X(), fOK.Y()-fgr); // coincide with module center
+
+  kk1 = ((fga+fga2)/4. + fga/6.)/2.; 
+
+  xk1 = fOK.X() - kk1*TMath::Sin(fTheta);
+  yk1 = fOK.Y() + kk1*TMath::Cos(fTheta) - fgr;
+  fOK3X3[0].Set(xk1,yk1);
+
+  xk2 = fOK.X() + kk1*TMath::Sin(fTheta);
+  yk2 = fOK.Y() - kk1*TMath::Cos(fTheta) - fgr;
+  fOK3X3[2].Set(xk2,yk2);
+
   // May 15, 2006; position of cell face of cells 
   fOB.Set(fOK.X()-fgb/2.*TMath::Cos(fTheta),  fOK.Y()-fgb/2.*TMath::Sin(fTheta)-fgr);
   fOB1.Set(fOB.X()-fga/4.*TMath::Sin(fTheta), fOB.Y()+fga/4.*TMath::Cos(fTheta));
   fOB2.Set(fOB.X()+fga/4.*TMath::Sin(fTheta), fOB.Y()-fga/4.*TMath::Cos(fTheta));
+
 }
 
 //_____________________________________________________________________________
@@ -160,15 +185,9 @@ void AliEMCALShishKebabTrd1Module::DefineFirstModule()
   Double_t xA = fga/2. + fga2/2., yA = fgr;
   fB = yA - fA*xA;
 
-  Double_t kk1 = (fga+fga2)/(2.*4.); // kk1=kk2 
-  fOK1.Set(fOK.X() - kk1, fOK.Y()-fgr);
-  fOK2.Set(fOK.X() + kk1, fOK.Y()-fgr);
-
-  fOB.Set(fOK.X(),fOK.Y()-fgb/2.-fgr);
-  fOB1.Set(fOB.X()-fga/4., fOB.Y());
-  fOB2.Set(fOB.X()+fga/4., fOB.Y());
-
   TObject::SetUniqueID(1); //
+
+  DefineAllStaff();
 }
 
 //_____________________________________________________________________________
@@ -195,7 +214,9 @@ Bool_t AliEMCALShishKebabTrd1Module::GetParameters()
   fgangle    = Double_t(fgGeometry->GetTrd1Angle())*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;
+}
index 8ee884c..0e89a47 100644 (file)
@@ -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
index 366d63b..6a9efd9 100644 (file)
@@ -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; iy<iyMax; iy++) { // flat in phi
@@ -862,13 +866,17 @@ void AliEMCALv0::CreateEmod(const char* mother, const char* child)
       }
     }
   }
-  AliDebug(2,Form(" Number of modules in Super Module %i \n", nr));
+  AliDebug(2,Form(" Number of modules in Super Module(%s) %i \n", mother, nr));
 }
 
 // 8-dec-04 by PAI
 void AliEMCALv0::Trd1Tower3X3(const double parSCM0[4])
 {
-  // PB should be for whole SCM0 - ?
+  // Fixed Nov 13,2006
+  printf(" AliEMCALv0::Trd1Tower3X3() : parSCM0");
+  for(int i=0; i<4; i++) printf(" %7.4f ", parSCM0[i]);
+  printf("\n"); 
+  // Nov 10, 2006 - different name of SCMX
   double parTRAP[11], *dummy=0;
   AliEMCALGeometry * g = GetGeometry(); 
   TString gn(g->GetName()), 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++) {
index 40c09ce..1e8715b 100644 (file)
@@ -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)