Tuning of EMCAL cluster finder for TRD1 geometry
authorpavlinov <pavlinov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 1 Mar 2006 21:15:40 +0000 (21:15 +0000)
committerpavlinov <pavlinov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 1 Mar 2006 21:15:40 +0000 (21:15 +0000)
12 files changed:
EMCAL/AliEMCALClusterizerv1.cxx
EMCAL/AliEMCALClusterizerv1.h
EMCAL/AliEMCALDigitizer.cxx
EMCAL/AliEMCALGeometry.cxx
EMCAL/AliEMCALGeometry.h
EMCAL/AliEMCALGeometryOfflineTrd1.cxx
EMCAL/AliEMCALRecPoint.cxx
EMCAL/AliEMCALRecPoint.h
EMCAL/AliEMCALSDigitizer.cxx
EMCAL/AliEMCALShishKebabTrd1Module.cxx
EMCAL/AliEMCALShishKebabTrd1Module.h
EMCAL/AliEMCALv0.cxx

index 3834137..c629f9c 100644 (file)
@@ -1,4 +1,3 @@
-
 /**************************************************************************
  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  *                                                                        *
 
 // --- ROOT system ---
 
-#include "TROOT.h" 
-#include "TFile.h" 
-#include "TFolder.h" 
-#include "TMath.h" 
-#include "TMinuit.h"
-#include "TTree.h" 
-#include "TSystem.h" 
-#include "TBenchmark.h"
-
+#include <TROOT.h>
+#include <TH1.h>
+#include <TFile.h> 
+#include <TFolder.h> 
+#include <TMath.h> 
+#include <TMinuit.h>
+#include <TTree.h> 
+#include <TSystem.h> 
+#include <TBenchmark.h>
+#include <TBrowser.h>
 // --- Standard library ---
 
 
 #include "AliEMCALDigitizer.h"
 #include "AliEMCAL.h"
 #include "AliEMCALGeometry.h"
+#include "AliEMCALHistoUtilities.h"
 
 ClassImp(AliEMCALClusterizerv1)
 
-Int_t addOn[20][60][60]; 
-
+AliEMCALGeometry * geom = 0;
 //____________________________________________________________________________
   AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer()
 {
@@ -81,16 +81,8 @@ Int_t addOn[20][60][60];
   
   InitParameters() ; 
   fDefaultInit = kTRUE ;
-  for(Int_t is=0;is<20;is++){ 
-    for(Int_t i=0;i<60;i++){ 
-      for(Int_t j=0;j<60;j++){ 
-       addOn[is][i][j]=0;
-      }
-    }
-  }
-//PH   cout<<"file to read 1"<<endl;
-  ReadFile();
-//PH   cout<<"file read 1"<<endl;
+  geom = AliEMCALGeometry::GetInstance();
+  geom->GetTransformationForSM(); // Global <-> Local
 }
 
 //____________________________________________________________________________
@@ -102,31 +94,18 @@ AliEMCALClusterizerv1::AliEMCALClusterizerv1(const TString alirunFileName, const
   InitParameters() ; 
   Init() ;
   fDefaultInit = kFALSE ; 
-  for(Int_t is=0;is<20;is++){ 
-    for(Int_t i=0;i<60;i++){ 
-      for(Int_t j=0;j<60;j++){ 
-       addOn[is][i][j]=0;
-      }
-    }
-  }
-//PH   cout<<"file to read 2"<<endl;
-  ReadFile();
-//PH   cout<<"file read 2"<<endl;
-
 }
 
 //____________________________________________________________________________
   AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
 {
   // dtor
-  
 }
 
 //____________________________________________________________________________
 const TString AliEMCALClusterizerv1::BranchName() const 
 { 
    return GetName();
-
 }
 
 //____________________________________________________________________________
@@ -185,11 +164,9 @@ void AliEMCALClusterizerv1::Exec(Option_t * option)
   if(strstr(option,"tim")){
     gBenchmark->Stop("EMCALClusterizer");
     printf("Exec took %f seconds for Clusterizing %f seconds per event", 
-        gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nEvents ) ;
+        gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nEvents );
   }
 
-  SaveHists();
-
 }
 
 //____________________________________________________________________________
@@ -226,8 +203,6 @@ Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * emcRP, AliEMCALDigit **
 
   Int_t iDigit ;
 
-  AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance(); //gime->EMCALGeometry() ; 
-
   for(iDigit = 0; iDigit < nDigits; iDigit++){
     digit = maxAt[iDigit]; 
 
@@ -309,18 +284,18 @@ void AliEMCALClusterizerv1::Init()
   // Make all memory allocations which can not be done in default constructor.
   // Attach the Clusterizer task to the list of EMCAL tasks
   
-  //AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance() ;
   AliRunLoader *rl = AliRunLoader::GetRunLoader();
-  AliEMCALGeometry *geom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
+  geom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
+  geom->GetTransformationForSM(); // Global <-> Local
   AliInfo(Form("geom 0x%x",geom));
 
 //Sub  fNTowers = geom->GetNZ() *  geom->GetNPhi() ;
-  fNTowers =400;
+  fNTowers =400; // ?
   if(!gMinuit) 
     gMinuit = new TMinuit(100) ;
  //Sub if ( !gime->Clusterizer() ) 
  //Sub   gime->PostClusterizer(this); 
- BookHists();
+ fHists = BookHists();
 }
 
 //____________________________________________________________________________
@@ -330,74 +305,46 @@ void AliEMCALClusterizerv1::InitParameters()
   fNumberOfECAClusters = 0;
 
   fECAClusteringThreshold   = 0.5;  // value obtained from Alexei
-  fECALocMaxCut = 0.03;
+  fECALocMaxCut = 0.03; // ??
 
   fECAW0     = 4.5 ;
   fTimeGate = 1.e-8 ; 
   fToUnfold = kFALSE ;
   fRecPointsInRun  = 0 ;
-  fMinECut = 0.3;
+  fMinECut = 0.01; // have to be tune
 }
 
 //____________________________________________________________________________
-Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2)const
+Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2) const
 {
   // Gives the neighbourness of two digits = 0 are not neighbour but continue searching 
   //                                       = 1 are neighbour
-  //                                       = 2 are not neighbour but do not continue searching
+  //                                       = 2 is in different SM, quit from 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  
 
-  //AliEMCALGeometry * geom = AliEMCALGetter::Instance()->EMCALGeometry() ;
-   AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance() ;
-
-  Int_t rv = 0 ; 
-
-  Int_t relid1[2] ; 
- //Sub geom->AbsToRelNumbering(d1->GetId(), relid1) ; 
-    Int_t nSupMod=0;
-    Int_t nTower=0;
-    Int_t nIphi=0;
-    Int_t nIeta=0;
-    Int_t iphi=0;
-    Int_t ieta=0;
-   geom->GetCellIndex(d1->GetId(), nSupMod,nTower,nIphi,nIeta);
-   geom->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta);
-   relid1[0]=ieta;
-   relid1[1]=iphi;
-
-
-    Int_t nSupMod1=0;
-    Int_t nTower1=0;
-    Int_t nIphi1=0;
-    Int_t nIeta1=0;
-    Int_t iphi1=0;
-    Int_t ieta1=0;
-   geom->GetCellIndex(d2->GetId(), nSupMod1,nTower1,nIphi1,nIeta1);
-   geom->GetCellPhiEtaIndexInSModule(nSupMod1, nTower1,nIphi1,nIeta1, iphi1,ieta1);
-   Int_t relid2[2] ; 
-   relid2[0]=ieta1;
-   relid2[1]=iphi1;
-
-  //Sub  geom->AbsToRelNumbering(d2->GetId(), relid2) ; 
-  
+  static Int_t rv; 
+  static Int_t nSupMod1=0, nTower1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
+  static Int_t nSupMod2=0, nTower2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
+  static Int_t rowdiff, coldiff;
+  rv = 0 ; 
 
-  Int_t rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ;  
-  Int_t coldiff = TMath::Abs( relid1[1] - relid2[1] ) ;  
+  geom->GetCellIndex(d1->GetId(), nSupMod1,nTower1,nIphi1,nIeta1);
+  geom->GetCellIndex(d2->GetId(), nSupMod2,nTower2,nIphi2,nIeta2);
+  if(nSupMod1 != nSupMod2) return 2; // different SM
+
+  geom->GetCellPhiEtaIndexInSModule(nSupMod1,nTower1,nIphi1,nIeta1, iphi1,ieta1);
+  geom->GetCellPhiEtaIndexInSModule(nSupMod2,nTower2,nIphi2,nIeta2, iphi2,ieta2);
+
+  rowdiff = TMath::Abs(iphi1 - iphi2);  
+  coldiff = TMath::Abs(ieta1 - ieta2) ;  
   
-  if (( coldiff <= 1 )  && ( rowdiff <= 1 )){
-      rv = 1 ; 
-  }
-  else {
-    if((relid2[0] > relid1[0]) && (relid2[1] > relid1[1]+1)) 
-      rv = 2; //  Difference in row numbers is too large to look further 
-  }
+  if (( coldiff <= 1 )  && ( rowdiff <= 1 )) rv = 1;  // neighbours with at least commom vertex
  
-  if (gDebug == 2 ) 
-if(rv==1)printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n", 
-        rv, d1->GetId(), relid1[0], relid1[1],
-        d2->GetId(), relid2[0], relid2[1]) ;   
+  if (gDebug == 2 && rv==1) 
+  printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n", 
+        rv, d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2);   
   
   return rv ; 
 }
@@ -437,8 +384,10 @@ void AliEMCALClusterizerv1::WriteRecPoints()
   
   aECARecPoints->Sort() ;
 
-  for(index = 0; index < aECARecPoints->GetEntries(); index++)
+  for(index = 0; index < aECARecPoints->GetEntries(); index++) {
     (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->SetIndexInList(index) ;
+    (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->Print();
+  }
 
   Int_t bufferSize = 32000 ;    
   Int_t splitlevel = 0 ; 
@@ -467,53 +416,52 @@ void AliEMCALClusterizerv1::MakeClusters()
 
   TObjArray * aECARecPoints = emcalLoader->RecPoints() ; 
     
-  AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
   if (geom==0) 
      AliFatal("Did not get geometry from EMCALLoader");
 
   aECARecPoints->Clear();
 
-  TClonesArray * digits = emcalLoader->Digits() ; 
-  TClonesArray * digitsC =  dynamic_cast<TClonesArray*>(digits->Clone()) ;
+  TClonesArray * digits  = emcalLoader->Digits() ; 
+  TClonesArray * digitsC =  dynamic_cast<TClonesArray*>(digits->Clone());
 
-  // Clusterization starts    
   TIter nextdigit(digitsC) ;
   AliEMCALDigit * digit;
+  double e=0.0, ehs = 0.0;
+  while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // clean up digits
+    e = Calibrate(digit->GetAmp());
+    AliEMCALHistoUtilities::FillH1(fHists, 10, digit->GetAmp());
+    AliEMCALHistoUtilities::FillH1(fHists, 11, e);
+    if(e < fMinECut ) digitsC->Remove(digit);
+    else              ehs += e;
+  }  
+  cout << " Number of digits " << digits->GetEntries() << " -> (e>" <<fMinECut <<")";
+  cout << digitsC->GetEntries()<< " ehs "<<ehs<<endl; 
 
-  cout << "Outer Loop" << endl;
+  // Clusterization starts    
+  //  cout << "Outer Loop" << endl;
+  int ineb=0;
+  nextdigit.Reset();
+  AliEMCALRecPoint * recPoint = 0 ; 
   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
-    AliEMCALRecPoint * clu = 0 ; 
-    
-    TArrayI clusterECAdigitslist(50000);   
-
-////////////////////////temp solution, embedding///////////////////
-   int nSupMod=0, nTower=0, nIphi=0, nIeta=0;
-   int iphi=0, ieta=0;
-   geom->GetCellIndex(digit->GetId(), nSupMod,nTower,nIphi,nIeta);
-   geom->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta);
-   
-   //cout << "Some checking" << endl;
-   if (geom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmp()+addOn[nSupMod-1][ieta-1][iphi-1]) > fECAClusteringThreshold  ) ){
-      //if(addOn[nSupMod-1][ieta-1][iphi-1]>0)cout<<"11 digit, add "<<ieta<<" "<<iphi<<" "<<addOn[nSupMod-1][ieta-1][iphi-1]<<" "<<digit->GetAmp()<<endl;
-//         cout<<"crossed the threshold "<<endl;
-      Int_t iDigitInECACluster = 0;
-      // start a new Tower RecPoint
-      if(fNumberOfECAClusters >= aECARecPoints->GetSize()) 
-         aECARecPoints->Expand(2*fNumberOfECAClusters+1) ;
+    recPoint = 0 ; 
+    TArrayI clusterECAdigitslist(1000); // what is this   
+
+    if(geom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmp()) > fECAClusteringThreshold  ) ){
+      Int_t iDigitInECACluster = 0; // start a new Tower RecPoint
+      if(fNumberOfECAClusters >= aECARecPoints->GetSize()) aECARecPoints->Expand(2*fNumberOfECAClusters+1) ;
       AliEMCALRecPoint * rp = new  AliEMCALRecPoint("") ; 
-      //rp->SetGeom(geom);
       aECARecPoints->AddAt(rp, fNumberOfECAClusters) ;
-      clu = dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(fNumberOfECAClusters)) ; 
+      recPoint = dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(fNumberOfECAClusters)) ; 
       fNumberOfECAClusters++ ; 
-      clu->AddDigit(*digit, Calibrate(digit->GetAmp()+addOn[nSupMod-1][ieta-1][iphi-1])) ; 
+      recPoint->AddDigit(*digit, Calibrate(digit->GetAmp())) ; 
       clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ;     
       iDigitInECACluster++ ; 
-//    cout<<" 1st setp:cluno, digNo "<<fNumberOfECAClusters<<" "<<iDigitInECACluster<<endl;
       digitsC->Remove(digit) ; 
-      AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , thre = %f \n", digit->GetId(),Calibrate(digit->GetAmp()), fECAClusteringThreshold));  
-      nextdigit.Reset() ;
+      AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , thre = %f \n", digit->GetId(),
+      Calibrate(digit->GetAmp()), fECAClusteringThreshold));  
+      nextdigit.Reset(); // will start from beggining
       
-      AliEMCALDigit * digitN ; 
+      AliEMCALDigit * digitN = 0; // digi neighbor
       Int_t index = 0 ;
 
       // Find the neighbours
@@ -521,41 +469,26 @@ void AliEMCALClusterizerv1::MakeClusters()
        digit =  (AliEMCALDigit*)digits->At(clusterECAdigitslist[index]);
        index++ ; 
         while ( (digitN = (AliEMCALDigit *)nextdigit())) { // scan over the reduced list of digits 
-         // cout<<"we have new digit "<<endl;
-          // check that the digit is above the min E Cut
-          geom->GetCellIndex(digitN->GetId(), nSupMod,nTower,nIphi,nIeta);
-          geom->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta);
-          if( Calibrate(digitN->GetAmp()+addOn[nSupMod-1][ieta-1][iphi-1]) < fMinECut  )  digitsC->Remove(digitN);
-         //cout<<" new digit above ECut "<<endl;
-         Int_t ineb = AreNeighbours(digit, digitN);       // call (digit,digitN) in THAT oder !!!!! 
-//     cout<<" new digit neighbour?? "<<ineb<<endl;
-         switch (ineb ) {
-          case 0 :   // not a neighbour
-           break ;
-         case 1 :   // are neighbours 
-           //if(addOn[nSupMod-1][ieta-1][iphi-1]>0)cout<<"22 digit, add "<<nSupMod<<" "<<ieta<<" "<<iphi<<" "<<addOn[nSupMod-1][ieta-1][iphi-1]<<" "<<digit->GetAmp()<<endl;
-           clu->AddDigit(*digitN, Calibrate( digitN->GetAmp()+addOn[nSupMod-1][ieta-1][iphi-1]) ) ;
-           clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ; 
-           iDigitInECACluster++ ; 
-//    cout<<"when neighbour: cluno, digNo "<<digit->GetId()<<" "<<fNumberOfECAClusters<<" "<<iDigitInECACluster<<endl;
-           digitsC->Remove(digitN) ;
-//         break ;
-//          case 2 :   // too far from each other
-//Subh     goto endofloop1;   
-//         cout<<"earlier go to end of loop"<<endl;   
+         //          if( Calibrate(digitN->GetAmp()) < fMinECut  )  digitsC->Remove(digitN);
+         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()) ) ;
+             clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ; 
+             iDigitInECACluster++ ; 
+             digitsC->Remove(digitN) ;
+              break ;
+             case 2 :   // different SM
+              break ;
          } // switch
-    //cout<<"in nextDigit loop "<<fNumberOfECAClusters<<" "<<iDigitInECACluster<<endl;
-       } // while digitN
-       //cout << "Done neighbout searching" << endl;
-       nextdigit.Reset() ; 
-      } // loop over ECA cluster
-    } // energy threshold
-    else if(Calibrate(digit->GetAmp()+addOn[nSupMod-1][ieta-1][iphi-1]) < fMinECut  ){
-      //if(addOn[nSupMod-1][ieta-1][iphi-1]>0)cout<<"33 digit, add  "<<ieta<<" "<<iphi<<" "<<addOn[nSupMod-1][ieta-1][iphi-1]<<" "<<digit->GetAmp()<<endl;
-      digitsC->Remove(digit);
+        } // scan over the reduced list of digits
+      } // scan over digits already in cluster
+      nextdigit.Reset() ;  // will start from beggining
     }
-    //cout<<"after endofloop: cluno, digNo "<<fNumberOfECAClusters<<endl;
   } // while digit  
+  if(recPoint) cout << "cl.e " << recPoint->GetEnergy() << endl; 
   delete digitsC ;
   AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,digits->GetEntriesFast())); 
 }
@@ -704,8 +637,10 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
     printf("\n-----------------------------------------------------------------------\n");
   }
 }
-  void AliEMCALClusterizerv1::BookHists()
+TList* AliEMCALClusterizerv1::BookHists()
 {
+  gROOT->cd();
+
        pointE = new TH1F("pointE","point energy", 2000, 0.0, 150.);
        pointL1 = new TH1F("pointL1","point L1", 1000, 0.0, 3.);
        pointL2 = new TH1F("pointL2","point L2", 1000, 0.0, 3.);
@@ -715,35 +650,21 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
        MaxE = new TH1F("maxE","Max point energy", 2000, 0.0, 150.);
        MaxL1 = new TH1F("maxL1","Max point L1", 1000, 0.0, 3.);
        MaxL2 = new TH1F("maxL2","Max point L2", 1000, 0.0, 3.);
-       MaxDis = new TH1F("maxDis","Max point Dis", 1000, 0.0, 3.);
+       MaxDis = new TH1F("maxDis","Max point Dis", 1000, 0.0, 3.); // 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
+
+  return AliEMCALHistoUtilities::MoveHistsToList("EmcalClusterizerv1ControlHists", kFALSE);
 }
-void AliEMCALClusterizerv1::SaveHists()
+
+void AliEMCALClusterizerv1::SaveHists(const char *fn)
 {
- recofile=new TFile("reco.root","RECREATE"); 
-       pointE->Write();
-       pointL1->Write();
-       pointL2->Write();
-       pointDis->Write();
-       pointMult->Write();
-       digitAmp->Write();
-      MaxE->Write();
-      MaxL1->Write();
-      MaxL2->Write();
-      MaxDis->Write();
-recofile->Close();
+  AliEMCALHistoUtilities::SaveListOfHists(fHists, fn, kTRUE);
 }
 
-void AliEMCALClusterizerv1::ReadFile()
+void AliEMCALClusterizerv1::Browse(TBrowser* b)
 {
-  return; // 3-jan-05
-       FILE *fp = fopen("hijing1.dat","r");
-        for(Int_t line=0;line<9113;line++){
-         Int_t eg,l1,l2,sm;
-         Int_t ncols0;
-         ncols0 = fscanf(fp,"%d %d %d %d",&sm,&l1,&l2,&eg);
-       // cout<<eg<<" "<<l1<<" "<<l2<<endl;
-        addOn[sm-1][l1-1][l2-1]=eg;
-        //addOn[sm-1][l1-1][l2-1]=0;
-       }
+  if(fHists) b->Add(fHists);
+  TTask::Browse(b);
 }
-
index 6567d67..21113a2 100644 (file)
@@ -24,7 +24,7 @@
 // --- AliRoot header files ---
 
 #include "AliEMCALClusterizer.h"
-#include "TH1F.h"
+class TH1F;
 class AliEMCALRecPoint ; 
 class AliEMCALDigit ;
 class AliEMCALDigitizer ;
@@ -38,6 +38,7 @@ public:
   AliEMCALClusterizerv1() ;         
   AliEMCALClusterizerv1(const TString alirunFileNameFile, const TString eventFolderName = AliConfig::GetDefaultEventFolderName());
   virtual ~AliEMCALClusterizerv1()  ;
+  virtual void Browse(TBrowser* b);
   
   virtual Int_t   AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2)const ; 
                                // Checks if digits are in neighbour cells 
@@ -71,13 +72,14 @@ public:
   void Unload() ; 
   virtual const char * Version() const { return "clu-v1" ; }  
  
-  void BookHists();
-  void SaveHists();
+  TList* BookHists();
+  void   SaveHists(const char *fn="reco.root");  //*MENU*
 protected:
 
   void           WriteRecPoints() ;
   virtual void   MakeClusters( ) ;            
 ///////////////////// 
+   TList  *fHists;   //!
    TH1F* pointE;
    TH1F* pointL1;
    TH1F* pointL2;
@@ -106,7 +108,7 @@ private:
                               AliEMCALDigit ** /*maxAt*/,
                               Float_t * /*maxAtEnergy*/ ) const; //Unfolds cluster using TMinuit package
   void           PrintRecPoints(Option_t * option) ;
-  TFile* recofile;
+  //  TFile* recofile;
 private:
 
   Bool_t  fDefaultInit;              //! Says if the task was created by defaut ctor (only parameters are initialized)
@@ -129,8 +131,6 @@ private:
   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
 
-  void ReadFile();
-    
   ClassDef(AliEMCALClusterizerv1,4)   // Clusterizer implementation version 1
 
 };
index d862f90..72e01ee 100644 (file)
@@ -525,6 +525,8 @@ Bool_t AliEMCALDigitizer::Init()
   
   //to prevent cleaning of this object while GetEvent is called
   emcalLoader->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
+
+  Print();
   
   return fInit ;    
 }
index cb061ed..26a7240 100644 (file)
@@ -1,4 +1,4 @@
-/**************************************************************************
+ /**************************************************************************
  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  *                                                                        *
  * Author: The ALICE Off-line Project.                                    *
 
 // --- AliRoot header files ---
 #include <assert.h>
+#include "Riostream.h"
+
 #include <TMath.h>
 #include <TVector3.h>
+#include <TArrayD.h>
 #include <TRegexp.h>
 #include <TObjArray.h>
 #include <TObjString.h>
-#include <assert.h>
+#include <TGeoManager.h>
+#include <TGeoNode.h>
+#include <TGeoMatrix.h>
 #include <TMatrixD.h>
 #include <TClonesArray.h>
 
@@ -45,6 +50,9 @@
 
 // --- EMCAL headers
 #include "AliEMCALGeometry.h"
+#include "AliEMCALShishKebabTrd1Module.h"
+//#include "AliRecPoint.h"
+#include "AliEMCALRecPoint.h"
 #include "AliEMCALDigit.h"
 
 ClassImp(AliEMCALGeometry)
@@ -63,19 +71,6 @@ int  nAdditionalOpts = sizeof(additionalOpts) / sizeof(char*);
 AliEMCALGeometry::~AliEMCALGeometry(void){
     // dtor
 }
-
-//______________________________________________________________________
-Bool_t AliEMCALGeometry::AreInSameTower(Int_t id1, Int_t id2) const {
-  // Find out whether two hits are in the same tower - have to be change
-  Int_t idmax = TMath::Max(id1, id2) ; 
-  Int_t idmin = TMath::Min(id1, id2) ;
-  if ( ((idmax - GetNZ() * GetNPhi()) == idmin ) || 
-       ((idmax - 2 * GetNZ() * GetNPhi()) == idmin ) )
-    return kTRUE ; 
-  else 
-    return kFALSE ; 
-}
-
 //______________________________________________________________________
 void AliEMCALGeometry::Init(void){
   // Initializes the EMCAL parameters
@@ -92,6 +87,8 @@ void AliEMCALGeometry::Init(void){
   name.ToUpper();
   fKey110DEG = 0;
   if(name.Contains("110DEG")) fKey110DEG = 1; // for GetAbsCellId
+  fShishKebabTrd1Modules = 0;
+  fTrd2AngleY = f2Trd2Dy2 = fEmptySpace = fTubsR = fTubsTurnAngle = 0;
 
   fNZ             = 114;       // granularity along Z (eta) 
   fNPhi           = 168;       // granularity in phi (azimuth)
@@ -101,6 +98,7 @@ void AliEMCALGeometry::Init(void){
   fArm1EtaMax     = +0.7;      // pseudorapidity, Ending EMCAL Eta position
   fIPDistance     = 454.0;      // cm, Radial distance to inner surface of EMCAL
   fPhiGapForSM    = 0.;         // cm, only for final TRD1 geometry
+  for(int i=0; i<12; i++) fMatrixOfSM[i] = 0;
 
   // geometry
   if(name.Contains("SHISH")){ // Only shahslyk now
@@ -226,6 +224,10 @@ void AliEMCALGeometry::Init(void){
     } else if(name.Contains("TRD")) { // 1-oct-04
       fShellThickness  = TMath::Sqrt(fLongModuleSize*fLongModuleSize + f2Trd1Dx2*f2Trd1Dx2);
       fShellThickness += fSteelFrontThick;
+      // Local coordinates
+      fParSM[0] = GetShellThickness()/2.;        
+      fParSM[1] = GetPhiModuleSize() * GetNPhi()/2.;
+      fParSM[2] = 350./2.;
     }
   }
 
@@ -238,7 +240,9 @@ void AliEMCALGeometry::Init(void){
   
   if (kTRUE) {
     printf("Init: geometry of EMCAL named %s is as follows:\n", name.Data());
-    printf( "               ECAL      : %d x (%f cm Pb, %f cm Sc) \n", GetNECLayers(), GetECPbRadThick(), GetECScintThick() ) ; 
+    printf( "               ECAL      : %d x (%f cm Pb, %f cm Sc) \n", 
+    GetNECLayers(), GetECPbRadThick(), GetECScintThick() ) ; 
+    printf("                fSampling %5.2f \n",  fSampling );
     if(name.Contains("SHISH")){
       printf(" fIPDistance       %6.3f cm \n", fIPDistance);
       if(fSteelFrontThick>0.) 
@@ -269,6 +273,8 @@ void AliEMCALGeometry::Init(void){
         printf(" fTubsTurnAngle  %7.4f\n", fTubsTurnAngle);
         printf(" fEmptySpace     %7.4f cm\n", fEmptySpace);
       } else if(name.Contains("TRD1") && name.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(name.Contains("110DEG"))printf(" Last two modules have size 10 degree in  phi (180<phi<190)\n");
       }
@@ -422,7 +428,8 @@ TClonesArray *  AliEMCALGeometry::FillTRU(const TClonesArray * digits) {
 AliEMCALGeometry *  AliEMCALGeometry::GetInstance(){ 
   // Returns the pointer of the unique instance
   
-  return static_cast<AliEMCALGeometry *>( fgGeom ) ; 
+  AliEMCALGeometry * rv = static_cast<AliEMCALGeometry *>( fgGeom );
+  return rv; 
 }
 
 //______________________________________________________________________
@@ -450,7 +457,7 @@ AliEMCALGeometry* AliEMCALGeometry::GetInstance(const Text_t* name,
          printf(name);  
        }else{
          rv = (AliEMCALGeometry *) fgGeom; 
-       } // end if
+       } // end 
     }  // end if fgGeom
     return rv; 
 }
@@ -741,7 +748,7 @@ Bool_t AliEMCALGeometry::IsInEMCAL(Double_t x, Double_t y, Double_t z) const {
 //
 // == Shish-kebab cases ==
 //
-Int_t AliEMCALGeometry::GetAbsCellId(Int_t nSupMod, Int_t nTower, Int_t nIphi, Int_t nIeta)
+Int_t AliEMCALGeometry::GetAbsCellId(Int_t nSupMod, Int_t nTower, Int_t nIphi, Int_t nIeta) const
 { // 27-aug-04; 
   // corr. 21-sep-04; 
   //       13-oct-05; 110 degree case
@@ -772,7 +779,7 @@ Int_t AliEMCALGeometry::GetAbsCellId(Int_t nSupMod, Int_t nTower, Int_t nIphi, I
   return id;
 }
 
-Bool_t  AliEMCALGeometry::CheckAbsCellId(Int_t ind)
+Bool_t  AliEMCALGeometry::CheckAbsCellId(Int_t ind) const
 { // 17-niv-04 - analog of IsInECA
    if(name.Contains("TRD")) {
      if(ind<=0 || ind > fNCells) return kFALSE;
@@ -780,7 +787,7 @@ Bool_t  AliEMCALGeometry::CheckAbsCellId(Int_t ind)
    } else return IsInECA(ind);
 }
 
-Bool_t AliEMCALGeometry::GetCellIndex(Int_t absId,Int_t &nSupMod,Int_t &nTower,Int_t &nIphi,Int_t &nIeta)
+Bool_t AliEMCALGeometry::GetCellIndex(Int_t absId,Int_t &nSupMod,Int_t &nTower,Int_t &nIphi,Int_t &nIeta) const
 { // 21-sep-04
   // 19-oct-05;
   static Int_t tmp=0, sm10=0;
@@ -805,7 +812,7 @@ Bool_t AliEMCALGeometry::GetCellIndex(Int_t absId,Int_t &nSupMod,Int_t &nTower,I
   return kTRUE;
 }
 
-void AliEMCALGeometry::GetTowerPhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower,  int &iphit, int &ietat)
+void AliEMCALGeometry::GetTowerPhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower,  int &iphit, int &ietat) const
 { // added nSupMod; have to check  - 19-oct-05 ! 
   static Int_t nphi;
 
@@ -817,7 +824,7 @@ void AliEMCALGeometry::GetTowerPhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower,
 }
 
 void AliEMCALGeometry::GetCellPhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower, Int_t nIphi, Int_t nIeta, 
-int &iphi, int &ieta)
+int &iphi, int &ieta) const
 { // added nSupMod; Nov 25, 05
   static Int_t iphit, ietat;
 
@@ -827,6 +834,161 @@ int &iphi, int &ieta)
   // iphi - have to change from 1 to fNPhi*fNPHIdiv
   iphi  = (iphit-1)*fNPHIdiv + nIphi;     // y(module) =  y(SM) 
 }
+
+Int_t  AliEMCALGeometry::GetSuperModuleNumber(Int_t absId)  const
+{
+  static Int_t nSupMod, nTower, nIphi, nIeta;
+  GetCellIndex(absId, nSupMod, nTower, nIphi, nIeta);
+  return nSupMod;
+} 
+
+// Methods for AliEMCALRecPoint - Feb 19, 2006
+Bool_t AliEMCALGeometry::RelPosCellInSModule(Int_t absId, Double_t &xr, Double_t &yr, Double_t &zr)
+{
+  static Int_t nSupMod, nTower, nIphi, nIeta, iphi, ieta;
+  if(!CheckAbsCellId(absId)) return kFALSE;
+
+  GetCellIndex(absId, nSupMod, nTower, nIphi, nIeta);
+  GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi, ieta); 
+  xr = fXCentersOfCells->At(ieta-1);
+  zr = fEtaCentersOfCells->At(ieta-1);
+
+  yr = fPhiCentersOfCells->At(iphi-1);
+
+  //  cout<<" absId "<<absId<<" iphi "<<iphi<<"ieta"<<ieta;
+  // cout<< " xr " << xr << " yr " << yr << " zr " << zr <<endl;
+  return kTRUE;
+}
+
+void AliEMCALGeometry::CreateListOfTrd1Modules()
+{
+  cout<< endl<< " AliEMCALGeometry::CreateListOfTrd1Modules() started " << endl;
+  AliEMCALShishKebabTrd1Module *mod=0, *mTmp=0; // current module
+  if(fShishKebabTrd1Modules == 0) {
+    fShishKebabTrd1Modules = new TList;
+    for(int iz=0; iz< GetNZ(); iz++) { 
+      if(iz==0) { 
+        mod  = new AliEMCALShishKebabTrd1Module(TMath::Pi()/2.,this);
+      } else {
+        mTmp  = new AliEMCALShishKebabTrd1Module(*mod);
+        mod   = mTmp;
+      }
+      fShishKebabTrd1Modules->Add(mod);
+    }
+  } else {
+    cout<<" Already exits : ";
+  }
+  cout<<" fShishKebabTrd1Modules "<< fShishKebabTrd1Modules << " has " 
+  << fShishKebabTrd1Modules->GetSize() << " modules" <<endl << endl;
+  // Feb 20,2006;
+  // 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);
+  printf(" Cells grid in eta directions : size %i\n", fEtaCentersOfCells->GetSize());
+  Int_t iphi=0, ieta=0, nTower=0;
+  Double_t xr, zr;
+  for(Int_t it=0; it<fNZ; it++) { // array index
+    AliEMCALShishKebabTrd1Module *trd1 = GetShishKebabModule(it);
+    nTower = fNPhi*it + 1;
+    for(Int_t ic=0; ic<fNETAdiv; ic++) { // array index
+      trd1->GetCenterOfCellInLocalCoordinateofSM(ic+1, xr, zr);
+      GetCellPhiEtaIndexInSModule(1, nTower, 1, ic+1, iphi, ieta); // don't depend from phi
+      fXCentersOfCells->AddAt(float(xr) - fParSM[0],ieta-1);
+      fEtaCentersOfCells->AddAt(float(zr) - fParSM[2],ieta-1);
+    }
+  }
+  for(Int_t i=0; i<fEtaCentersOfCells->GetSize(); i++) {
+    printf(" ind %2.2i : z %8.3f : x %8.3f", i+1, fEtaCentersOfCells->At(i),fXCentersOfCells->At(i));
+    if(i%2 != 0) printf("\n"); 
+  }
+  printf("\n"); 
+ // define grid for cells in phi(y) direction in local coordinates system of SM
+  fPhiCentersOfCells = new TArrayD(fNPhi*fNPHIdiv);
+  printf(" 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);
+      printf(" ind %2.2i : y %8.3f ", ind, fPhiCentersOfCells->At(ind++));
+      if(ic == fNPHIdiv-1) printf("\n"); 
+    }
+  }
+  printf("\n"); 
+}
+
+void  AliEMCALGeometry::GetTransformationForSM()
+{
+  static Bool_t transInit=kFALSE;
+  if(transInit) return;
+
+  int i=0;
+  if(gGeoManager == 0) {
+    Info("CreateTransformationForSM() "," Load geometry : TGeoManager::Import()");
+    assert(0);
+  }
+  TGeoNode *tn = gGeoManager->GetTopNode();
+  TGeoNode *node=0, *XEN1 = 0;
+  for(i=0; i<tn->GetNdaughters(); i++) {
+    node = tn->GetDaughter(i);
+    TString ns(node->GetName());
+    if(ns.Contains(GetNameOfEMCALEnvelope())) {
+      XEN1 = node;
+      break;
+    }
+  }
+  if(!XEN1) {
+    Info("CreateTransformationForSM() "," geometry has not EMCAL envelope with name %s", 
+    GetNameOfEMCALEnvelope());
+    assert(0);
+  }
+  printf(" i %i : EMCAL Envelope is %s : #SM %i \n", i, XEN1->GetName(), XEN1->GetNdaughters());
+  for(i=0; i<XEN1->GetNdaughters(); i++) {
+    TGeoNodeMatrix *sm = (TGeoNodeMatrix*)XEN1->GetDaughter(i);
+    fMatrixOfSM[i] = sm->GetMatrix();
+    printf(" %i : matrix %x \n", i, fMatrixOfSM[i]);
+  }
+  transInit = kTRUE;
+}
+
+void AliEMCALGeometry::GetGlobal(const Double_t *loc, Double_t *glob, int nsm) const
+{
+  //  if(fMatrixOfSM[0] == 0) GetTransformationForSM();
+  static int ind;
+  ind = nsm-1;
+  if(ind>=0 && ind < GetNumberOfSuperModules()) {
+    fMatrixOfSM[ind]->LocalToMaster(loc, glob);
+  }
+}
+
+void AliEMCALGeometry::GetGlobal(const Int_t absId, TVector3 &vglob) const
+{ // have to be defined  
+}
+
+void AliEMCALGeometry::GetGlobal(const TVector3 &vloc, TVector3 &vglob, int nsm) const
+{
+  static Double_t tglob[3], tloc[3];
+  vloc.GetXYZ(tloc);
+  GetGlobal(tloc, tglob, nsm);
+  vglob.SetXYZ(tglob[0], tglob[1], tglob[2]);
+}
+
+void AliEMCALGeometry::GetGlobal(const AliRecPoint *rp, TVector3 &vglob) const
+{
+  static TVector3 vloc;
+  static Int_t nSupMod, nTower, nIphi, nIeta;
+
+  AliRecPoint *rpTmp = (AliRecPoint*)rp; // const_cast ??
+  if(!rpTmp) return;
+  AliEMCALRecPoint *rpEmc = (AliEMCALRecPoint*)rpTmp;
+
+  GetCellIndex(rpEmc->GetAbsId(0), nSupMod, nTower, nIphi, nIeta);
+  rpTmp->GetLocalPosition(vloc);
+  GetGlobal(vloc, vglob, nSupMod);
+}
+
 // Service routine 
 int  AliEMCALGeometry::ParseString(const TString &topt, TObjArray &Opt)
 { // Feb 06, 2006
index d008031..45c4add 100644 (file)
 // --- ROOT system ---
 class TString ;
 class TObjArray;
-class TVector3 ;
+class TVector3;
+class TGeoMatrix;
+class TArrayD;
 class TParticle ; 
+class AliEMCALShishKebabTrd1Module;
+class AliEMCALRecPoint;
 class TClonesArray ;
 
 // --- AliRoot header files ---
@@ -40,12 +44,21 @@ public:
     return *(GetInstance()) ; 
   };
 
-  Bool_t AreInSameTower(Int_t id1, Int_t id2) const ;  
+  // Have to call GetTransformationForSM() before calculation global charachteristics 
+  void GetGlobal(const Double_t *loc, Double_t *glob, int nsm) const;
+  void GetGlobal(const TVector3 &vloc, TVector3 &vglob, int nsm) const;
+  void GetGlobal(const 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 index,Float_t &eta,Float_t &phi) const;
+
+  //  virtual void GetGlobal(const AliEMCALRecPoint *rp, TVector3 &vglob) const;
+
+  virtual void GetGlobal(const AliRecPoint *rp, TVector3 &vglob) const;
+  // Bool_t AreInSameTower(Int_t id1, Int_t id2) const ;  
 
   TClonesArray *  FillTRU(const TClonesArray * digits)  ;
-  
   virtual void GetGlobal(const AliRecPoint *, TVector3 &, TMatrixF &) const {}
-  virtual void GetGlobal(const AliRecPoint *, TVector3 &) const {}
+
   virtual Bool_t Impact(const TParticle *) const {return kTRUE;}
 
   Bool_t IsInEMCAL(Double_t x, Double_t y, Double_t z) const;
@@ -53,6 +66,7 @@ public:
   Bool_t  IsInitialized(void) const { return fgInit ; }
   // Return EMCAL geometrical parameters
   // geometry
+  Char_t* GetNameOfEMCALEnvelope() const {return "XEN1";}
   Float_t GetAlFrontThickness() const { return fAlFrontThick;}
   Float_t GetArm1PhiMin() const { return fArm1PhiMin ; }
   Float_t GetArm1PhiMax() const { return fArm1PhiMax ; }
@@ -78,8 +92,8 @@ public:
   Float_t GetSampling() const {return fSampling ; } 
   Bool_t IsInECA(Int_t index) const { if ( (index > 0 && (index <= GetNZ() * GetNPhi()))) return kTRUE; else return kFALSE ;}
 
-  Int_t   GetNumberOfSuperModules() {return fNumberOfSuperModules;}
-  Float_t GetfPhiGapForSuperModules() {return fPhiGapForSM;}
+  Int_t   GetNumberOfSuperModules() const {return fNumberOfSuperModules;}
+  Float_t GetfPhiGapForSuperModules() const {return fPhiGapForSM;}
   Float_t GetPhiModuleSize() const  {return fPhiModuleSize;}
   Float_t GetEtaModuleSize() const  {return fEtaModuleSize;}
   Float_t GetFrontSteelStrip() const {return fFrontSteelStrip;}
@@ -105,18 +119,38 @@ public:
   Float_t Get2Trd2Dy2()  const {return f2Trd2Dy2;}
   Float_t GetTubsR()     const {return fTubsR;}
   Float_t GetTubsTurnAngle() const {return fTubsTurnAngle;}
-  // Dabs id <-> indexes; Shish-kebab case 
-  Int_t   GetAbsCellId(Int_t nSupMod, Int_t nTower, Int_t nIphi, Int_t nIeta);
-  Bool_t  GetCellIndex(Int_t absId, Int_t &nSupMod, Int_t &nTower, Int_t &nIphi, Int_t &nIeta);
-  void    GetTowerPhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower, Int_t &iphit, Int_t &ietat);
+
+  // TRD1 staff
+  void    CreateListOfTrd1Modules();
+  TList  *GetShishKebabTrd1Modules() const {return fShishKebabTrd1Modules;}
+  AliEMCALShishKebabTrd1Module *GetShishKebabModule(const Int_t neta=0)
+  {
+    static AliEMCALShishKebabTrd1Module* trd1=0;
+    if(fShishKebabTrd1Modules && neta>=0 && neta<fShishKebabTrd1Modules->GetSize()) {
+      trd1 = (AliEMCALShishKebabTrd1Module*)fShishKebabTrd1Modules->At(neta);
+    } else trd1 = 0;
+    return trd1;
+  }
+  void     GetTransformationForSM();
+  Float_t *GetSuperModulesPars() {return fParSM;}
+  TGeoMatrix *GetTransformationForSM(int i) {
+  if(i>=0 && GetNumberOfSuperModules()) return fMatrixOfSM[i]; 
+                                        else return 0;}
+  // abs id <-> indexes; Shish-kebab case (TRD1 or TRD2)
+  Int_t   GetAbsCellId(Int_t nSupMod, Int_t nTower, Int_t nIphi, Int_t nIeta) const;
+  Bool_t  CheckAbsCellId(Int_t ind) const; // replace the IsInECA
+  Bool_t  GetCellIndex(Int_t absId, Int_t &nSupMod, Int_t &nTower, Int_t &nIphi, Int_t &nIeta) const;
+  void    GetTowerPhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower, Int_t &iphit, Int_t &ietat) const;
   void    GetCellPhiEtaIndexInSModule(Int_t nSupMod, Int_t nTower, Int_t nIphi, Int_t nIeta,
-                                      Int_t &iphi, Int_t &ieta);
-  Bool_t  CheckAbsCellId(Int_t ind); // replace the IsInECA
+                                      Int_t &iphi, Int_t &ieta) const ;
+  Int_t   GetSuperModuleNumber(Int_t absId)  const; 
+  // Methods for AliEMCALRecPoint - Frb 19, 2006
+  Bool_t   RelPosCellInSModule(Int_t absId, Double_t &xr, Double_t &yr, Double_t &zr);
   // ---
-  Float_t AngleFromEta(Float_t eta){ // returns theta in radians for a given pseudorapidity
+  Float_t AngleFromEta(Float_t eta) const { // returns theta in radians for a given pseudorapidity
     return 2.0*TMath::ATan(TMath::Exp(-eta));
   }
-  Float_t ZFromEtaR(Float_t r,Float_t eta){ // returns z in for a given
+  Float_t ZFromEtaR(Float_t r,Float_t eta) const { // returns z in for a given
     // pseudorapidity and r=sqrt(x*x+y*y).
     return r/TMath::Tan(AngleFromEta(eta));
   }
@@ -124,8 +158,6 @@ public:
   Int_t TowerIndex(Int_t iz,Int_t iphi) const; // returns tower index
        // returns tower indexs iz, iphi.
   void TowerIndexes(Int_t index,Int_t &iz,Int_t &iphi) const;
-       // for a given tower index it returns eta and phi of center of that tower.
-  void EtaPhiFromIndex(Int_t index,Float_t &eta,Float_t &phi) const;
        // returns x, y, and z (cm) on the inner surface of a given EMCAL Cell specified by relid.
   void XYZFromIndex(const Int_t *relid,Float_t &x,Float_t &y, Float_t &z) const;
   void XYZFromIndex(Int_t absid, TVector3 &v) const;
@@ -149,10 +181,12 @@ protected:
   AliEMCALGeometry(const Text_t* name, const Text_t* title="") :
     AliGeometry(name, title) {// ctor only for internal usage (singleton)
     Init();
+    CreateListOfTrd1Modules();
   };
   AliEMCALGeometry() :
     AliGeometry() {// ctor only for internal usage (singleton)
-    Init();
+    CreateListOfTrd1Modules();
+    //Init();
   };
   void Init(void);                             // initializes the parameters of EMCAL
   void CheckAditionalOptions();              //
@@ -216,6 +250,15 @@ private:
   // Super module as TUBS
   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)
+  // Move from AliEMCALv0 - Feb 19, 2006
+  TList *fShishKebabTrd1Modules; //! list of modules
+  // Local coordinates of SM for TRD1
+  Float_t     fParSM[3];       // SM sizes as in GEANT (TRD1)
+  TGeoMatrix* fMatrixOfSM[12]; //![fNumberOfSuperModules]; get from gGeoManager;
   // Service routine 
   static int ParseString(const TString &topt, TObjArray &Opt);
 
index 0e11e81..45988d9 100644 (file)
@@ -53,12 +53,12 @@ void AliEMCALGeometryOfflineTrd1::Init()
   fTrd1Modules[0] =  new AliEMCALShishKebabTrd1Module();
   fSMMaxEta = 2*fMaxInEta;
   fSMPositionEta = new TVector2[fSMMaxEta];
-  fSMPositionEta[0] = fTrd1Modules[0]->GetCenterOfCell(1);
-  fSMPositionEta[1] = fTrd1Modules[0]->GetCenterOfCell(2);
+  fSMPositionEta[0] = fTrd1Modules[0]->GetCenterOfCellInLocalCoordinateofSM(1);
+  fSMPositionEta[1] = fTrd1Modules[0]->GetCenterOfCellInLocalCoordinateofSM(2);
   for(Int_t i=1; i<fMaxInEta; i++) {
     fTrd1Modules[i] = new AliEMCALShishKebabTrd1Module(*(fTrd1Modules[i-1]));
-    fSMPositionEta[2*i]   = fTrd1Modules[i]->GetCenterOfCell(1);
-    fSMPositionEta[2*i+1] = fTrd1Modules[i]->GetCenterOfCell(2);
+    fSMPositionEta[2*i]   = fTrd1Modules[i]->GetCenterOfCellInLocalCoordinateofSM(1);
+    fSMPositionEta[2*i+1] = fTrd1Modules[i]->GetCenterOfCellInLocalCoordinateofSM(2);
   }
   // PHI direction
   fSMPositionPhi.Set(2*fGeometry->GetNPhi());
index db3d007..5d89dec 100644 (file)
 //*-- Author: Heather Gray (LBL) merged AliEMCALRecPoint and AliEMCALTowerRecPoint 02/04
 
 // --- ROOT system ---
-#include "TPad.h"
-#include "TGraph.h"
-#include "TPaveText.h"
-#include "TClonesArray.h"
-#include "TMath.h"
+#include <Riostream.h>
+#include <TPad.h>
+#include <TGraph.h>
+#include <TPaveText.h>
+#include <TClonesArray.h>
+#include <TMath.h>
 
 // --- Standard library ---
 
@@ -37,6 +38,8 @@
 
 ClassImp(AliEMCALRecPoint)
 
+AliEMCALGeometry * geom = 0;
+
 //____________________________________________________________________________
 AliEMCALRecPoint::AliEMCALRecPoint()
   : AliRecPoint()
@@ -49,10 +52,13 @@ AliEMCALRecPoint::AliEMCALRecPoint()
   fAmp   = 0. ;   
   fCoreEnergy = 0 ; 
   fEnergyList = 0 ;
+  fAbsIdList  = 0;
   fParentsList = 0;
   fTime = 0. ;
-  fLocPos.SetX(0.)  ;      //Local position should be evaluated
+  //  fLocPos.SetX(1.e+6)  ;      //Local position should be evaluated
   fCoreRadius = 10;        //HG Check this
+  geom = AliEMCALGeometry::GetInstance();
+  geom->GetTransformationForSM(); // Global <-> Local
 }
 
 //____________________________________________________________________________
@@ -66,10 +72,13 @@ AliEMCALRecPoint::AliEMCALRecPoint(const char * opt) : AliRecPoint(opt)
   fAmp   = 0. ;   
   fCoreEnergy = 0 ; 
   fEnergyList = 0 ;
+  fAbsIdList  = 0;
   fParentsList = new Int_t[fMaxParent];
   fTime = -1. ;
-  fLocPos.SetX(1000000.)  ;      //Local position should be evaluated
+  //fLocPos.SetX(1.e+6)  ;      //Local position should be evaluated
   fCoreRadius = 10;        //HG Check this
+  geom = AliEMCALGeometry::GetInstance();
+  geom->GetTransformationForSM(); // Global <-> Local
 }
 //____________________________________________________________________________
 AliEMCALRecPoint::~AliEMCALRecPoint()
@@ -77,6 +86,8 @@ AliEMCALRecPoint::~AliEMCALRecPoint()
   // dtor
   if ( fEnergyList )
     delete[] fEnergyList ; 
+  if ( fAbsIdList )
+    delete[] fAbsIdList ; 
    if ( fParentsList)
     delete[] fParentsList;
 }
@@ -89,16 +100,22 @@ void AliEMCALRecPoint::AddDigit(AliEMCALDigit & digit, Float_t Energy)
   
   if(fEnergyList == 0)
     fEnergyList =  new Float_t[fMaxDigit]; 
+  if(fAbsIdList == 0) {
+    fAbsIdList =  new Int_t[fMaxDigit];
+    fSuperModuleNumber = geom->GetSuperModuleNumber(digit.GetId());
+  }
 
   if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists 
     fMaxDigit*=2 ; 
     Int_t * tempo = new Int_t[fMaxDigit]; 
     Float_t * tempoE =  new Float_t[fMaxDigit];
+    Int_t * tempoId = new Int_t[fMaxDigit]; 
 
     Int_t index ;     
     for ( index = 0 ; index < fMulDigit ; index++ ){
-      tempo[index]  = fDigitsList[index] ;
-      tempoE[index] = fEnergyList[index] ; 
+      tempo[index]   = fDigitsList[index] ;
+      tempoE[index]  = fEnergyList[index] ; 
+      tempoId[index] = fAbsIdList[index] ; 
     }
     
     delete [] fDigitsList ; 
@@ -107,17 +124,23 @@ void AliEMCALRecPoint::AddDigit(AliEMCALDigit & digit, Float_t Energy)
     delete [] fEnergyList ;
     fEnergyList =  new Float_t[fMaxDigit];
 
+    delete [] fAbsIdList ;
+    fAbsIdList =  new Int_t[fMaxDigit];
+
     for ( index = 0 ; index < fMulDigit ; index++ ){
       fDigitsList[index] = tempo[index] ;
       fEnergyList[index] = tempoE[index] ; 
+      fAbsIdList[index]  = tempoId[index] ; 
     }
  
     delete [] tempo ;
     delete [] tempoE ; 
+    delete [] tempoId ; 
   } // if
   
   fDigitsList[fMulDigit]   = digit.GetIndexInList()  ; 
   fEnergyList[fMulDigit]   = Energy ;
+  fAbsIdList[fMulDigit]    = digit.GetId();
   fMulDigit++ ; 
   fAmp += Energy ; 
 
@@ -128,38 +151,27 @@ Bool_t AliEMCALRecPoint::AreNeighbours(AliEMCALDigit * digit1, AliEMCALDigit * d
   // Tells if (true) or not (false) two digits are neighbours
   // A neighbour is defined as being two digits which share a corner
   
-  Bool_t areNeighbours = kFALSE ;
-  
-  //AliEMCALGeometry * geom =  (AliEMCALGetter::Instance())->EMCALGeometry();
-  AliEMCALGeometry * geom =  AliEMCALGeometry::GetInstance();
+  static Bool_t areNeighbours = kFALSE ;
+  static Int_t nSupMod=0, nTower=0, nIphi=0, nIeta=0;
+  static int nSupMod1=0, nTower1=0, nIphi1=0, nIeta1=0;
+  static Int_t relid1[2] , relid2[2] ; // ieta, iphi
+  static Int_t rowdiff=0, coldiff=0;
 
-  Int_t relid1[2] ; 
-    //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;
-       geom->GetCellIndex(digit1->GetId(), nSupMod,nTower,nIphi,nIeta);
-       geom->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta);
-       relid1[0]=ieta;
-       relid1[1]=iphi;
-//  geom->AbsToRelNumbering(digit1->GetId(), relid1) ; 
-
-  Int_t relid2[2] ; 
-    //copied for shish-kebab geometry, ieta,iphi is cast as float as eta,phi conversion
-    // for this geometry does not exist
-    int nSupMod1=0, nTower1=0, nIphi1=0, nIeta1=0;
-    int iphi1=0, ieta1=0;
-       geom->GetCellIndex(digit2->GetId(), nSupMod1,nTower1,nIphi1,nIeta1);
-       geom->GetCellPhiEtaIndexInSModule(nSupMod1,nTower1,nIphi1,nIeta1, iphi1,ieta1);
-       relid2[0]=ieta1;
-       relid2[1]=iphi1;
-//  geom->AbsToRelNumbering(digit2->GetId(), relid2) ; 
+  areNeighbours = kFALSE ;
+
+  //  AliEMCALGeometry * geom =  AliEMCALGeometry::GetInstance();
+
+  geom->GetCellIndex(digit1->GetId(), nSupMod,nTower,nIphi,nIeta);
+  geom->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, relid1[0],relid1[1]);
+
+  geom->GetCellIndex(digit2->GetId(), nSupMod1,nTower1,nIphi1,nIeta1);
+  geom->GetCellPhiEtaIndexInSModule(nSupMod1,nTower1,nIphi1,nIeta1, relid2[0],relid2[1]);
   
-  Int_t rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ;  
-  Int_t coldiff = TMath::Abs( relid1[1] - relid2[1] ) ;  
+  rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ;  
+  coldiff = TMath::Abs( relid1[1] - relid2[1] ) ;  
 
   if (( coldiff <= 1 )  && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0)) 
-    areNeighbours = kTRUE ;
+  areNeighbours = kTRUE ;
   
   return areNeighbours;
 }
@@ -329,124 +341,99 @@ void  AliEMCALRecPoint::EvalDispersion(Float_t logWeight, TClonesArray * digits)
 {
   // Calculates the dispersion of the shower at the origin of the RecPoint
 
-  Float_t d    = 0. ;
-  Float_t wtot = 0. ;
-
+  Double_t d = 0., wtot = 0., w = 0., xyzi[3], diff=0.;
+  Int_t iDigit=0, nstat=0, i=0;
   AliEMCALDigit * digit ;
-  //AliEMCALGeometry * geom = (AliEMCALGetter::Instance())->EMCALGeometry();
-  AliEMCALGeometry * geom =  AliEMCALGeometry::GetInstance();
+  //  AliEMCALGeometry * geom =  AliEMCALGeometry::GetInstance();
 
   // Calculates the centre of gravity in the local EMCAL-module coordinates   
-  Int_t iDigit;
-
-  if (!fLocPos.X() || !fLocPos.Y()) 
+  if (!fLocPos.Mag()) 
     EvalLocalPosition(logWeight, digits) ;
   
-  //Sub const Float_t kDeg2Rad = TMath::DegToRad() ; 
-    
-  Float_t cluEta =  fLocPos.X() ; 
-  Float_t cluPhi =  fLocPos.Y() ; 
-  Float_t cluR =  fLocPos.Z() ; 
-  
-  if (gDebug == 2) 
-    printf("EvalDispersion: eta,phi,r = %f,%f,%f", cluEta, cluPhi, cluR) ;
-  
   // Calculates the dispersion in coordinates 
-  wtot = 0.;
   for(iDigit=0; iDigit < fMulDigit; iDigit++) {
     digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit])  ;
-    Float_t etai = 0.;
-    Float_t phii = 0.;
-    //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;
-       geom->GetCellIndex(digit->GetId(), nSupMod,nTower,nIphi,nIeta);
-       geom->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta);
-       etai=(Float_t)ieta;
-       phii=(Float_t)iphi;
-       //        printf("%f,%d,%d \n", fAmp, ieta, iphi) ;
-
-/* Sub
-  geom->EtaPhiFromIndex(digit->GetId(), etai, phii);
-    phii = phii * kDeg2Rad;
-*/
-///////////////////////////
-//  if(fAmp>0)printf("%f %d %f", fAmp,iDigit,fEnergyList[iDigit]) ;
-/////////////////////////
 
-    if (gDebug == 2) 
-      printf("EvalDispersion: id = %d, etai,phii = %f,%f", digit->GetId(), etai, phii) ;
+    geom->RelPosCellInSModule(digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
+    w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
 
-    Float_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
-    d += w * ( (etai-cluEta)*(etai-cluEta) + (phii-cluPhi)*(phii-cluPhi) ) ; 
-    wtot+=w ;
+    if(w>0.0) {
+      wtot += w;
+      nstat++;
+      for(i=0; i<3; i++ ) {
+        diff = xyzi[i] - double(fLocPos[i]);
+        d   += w * diff*diff;
+      }
+    }
   }
   
-  if ( wtot > 0 ) 
-    d /= wtot ;
-  else 
-    d = 0. ; 
+  if ( wtot > 0 && nstat>1) d /= wtot ;
+  else                      d = 0. ; 
 
   fDispersion = TMath::Sqrt(d) ;
-  //      printf("Dispersion: = %f", fDispersion) ;
 }
 
 //____________________________________________________________________________
 void AliEMCALRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digits)
 {
   // Calculates the center of gravity in the local EMCAL-module coordinates 
-  Float_t wtot = 0. ;
-  //  Int_t relid[3] ;
+  //  Info("Print", " logWeight %f : cluster energy %f ", logWeight, fAmp); // for testing
   
-  AliEMCALDigit * digit ;
-  //AliEMCALGeometry * geom  =  (AliEMCALGetter::Instance())->EMCALGeometry();
-  AliEMCALGeometry * geom =  AliEMCALGeometry::GetInstance();
-  Int_t iDigit;
-  Float_t cluEta = 0;
-  Float_t cluPhi = 0;
- //Sub const Float_t kDeg2Rad = TMath::DegToRad();
+  AliEMCALDigit * digit;
+  //  AliEMCALGeometry * geom =  AliEMCALGeometry::GetInstance();
+  Int_t i=0, nstat=0;
+  Double_t clXYZ[3]={0.,0.,0.}, clRmsXYZ[3]={0.,0.,0.}, xyzi[3], wtot=0., w=0.;
 
-  for(iDigit=0; iDigit<fMulDigit; iDigit++) {
+  for(Int_t iDigit=0; iDigit<fMulDigit; iDigit++) {
     digit = dynamic_cast<AliEMCALDigit *>(digits->At(fDigitsList[iDigit])) ;
 
-    Float_t etai ;
-    Float_t phii ;
-    //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;
-       geom->GetCellIndex(digit->GetId(), nSupMod,nTower,nIphi,nIeta);
-       geom->GetCellPhiEtaIndexInSModule(nSupMod, nTower,nIphi,nIeta, iphi,ieta); //19-oct-05
-       etai=(Float_t)ieta;
-       phii=(Float_t)iphi;
-//Sub    geom->EtaPhiFromIndex(digit->GetId(), etai, phii);
-//Sub    phii = phii * kDeg2Rad;
-    Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
-    cluEta += (etai * w) ;
-    cluPhi += (phii * w );
-    wtot += w ;
-  }
+    geom->RelPosCellInSModule(digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
+    // printf(" Id %i : Local x,y,z %f %f %f \n", digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
+
+    if(logWeight > 0.0) w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ));
+    else                w = fEnergyList[iDigit]; // just energy
 
+    if(w>0.0) {
+      wtot += w ;
+      nstat++;
+      for(i=0; i<3; i++ ) {
+        clXYZ[i]    += (w*xyzi[i]);
+        clRmsXYZ[i] += (w*xyzi[i]*xyzi[i]);
+      }
+    }
+  }
+  //  cout << " wtot " << wtot << endl;
   if ( wtot > 0 ) { 
-    cluEta /= wtot ;
-    cluPhi /= wtot ;
+    //    xRMS   = TMath::Sqrt(x2m - xMean*xMean);
+    for(i=0; i<3; i++ ) {
+      clXYZ[i] /= wtot;
+      if(nstat>1) {
+        clRmsXYZ[i] /= (wtot*wtot);  
+        clRmsXYZ[i]  =  clRmsXYZ[i] - clXYZ[i]*clXYZ[i];
+       if(clRmsXYZ[i] > 0.0) {
+          clRmsXYZ[i] =  TMath::Sqrt(clRmsXYZ[i]);
+        } else      clRmsXYZ[i] = 0;
+      } else        clRmsXYZ[i] = 0;
+    }    
   } else {
-    cluEta = -1 ; 
-    cluPhi = -1.; 
+    for(i=0; i<3; i++ ) {
+      clXYZ[i] = clRmsXYZ[i] = -1.;
+    }
   }
-  
-  fLocPos.SetX(cluEta);
-  fLocPos.SetY(cluPhi);
-  fLocPos.SetZ(geom->GetIP2ECASection());
+  // clRmsXYZ[i] ??
+  fLocPos.SetX(clXYZ[0]);
+  fLocPos.SetY(clXYZ[1]);
+  fLocPos.SetZ(clXYZ[2]);
     
 //  if (gDebug==2)
 //    printf("EvalLocalPosition: eta,phi,r = %f,%f,%f", fLocPos.X(), fLocPos.Y(), fLocPos.Z()) ; 
-  fLocPosM = 0 ;
+  fLocPosM = 0 ; // covariance matrix
 }
 
+//void AliEMCALRecPoint::EvalLocalPositionSimple()
+//{ // Weight is proportional of cell energy 
+//}
+
 //______________________________________________________________________________
 void AliEMCALRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits)
 {
@@ -457,12 +444,11 @@ void AliEMCALRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits)
 
   AliEMCALDigit * digit ;
   const Float_t kDeg2Rad = TMath::DegToRad() ;
-  //AliEMCALGeometry * geom = (AliEMCALGetter::Instance())->EMCALGeometry();    
-  AliEMCALGeometry * geom =  AliEMCALGeometry::GetInstance();
+  //  AliEMCALGeometry * geom =  AliEMCALGeometry::GetInstance();
 
   Int_t iDigit;
 
-  if (!fLocPos.X() || !fLocPos.Y() ) {
+  if (!fLocPos.Mag()) {
     EvalLocalPosition(logWeight, digits);
   }
   
@@ -494,8 +480,7 @@ void  AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
   const Float_t kDeg2Rad = TMath::DegToRad();
   AliEMCALDigit * digit ;
 
-  //AliEMCALGeometry * geom = (AliEMCALGetter::Instance())->EMCALGeometry();
-  AliEMCALGeometry * geom =  AliEMCALGeometry::GetInstance();
+  //  AliEMCALGeometry * geom =  AliEMCALGeometry::GetInstance();
   TString gn(geom->GetName());
 
   Int_t iDigit;
@@ -504,7 +489,7 @@ void  AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
     digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit])  ;
     Float_t etai = 0. ;
     Float_t phii = 0. ; 
-    if(gn.Contains("SHISH")) {
+    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;
@@ -668,7 +653,6 @@ void  AliEMCALRecPoint::EvalParents(TClonesArray * digits)
 void AliEMCALRecPoint::GetLocalPosition(TVector3 & lpos) const
 {
   // returns the position of the cluster in the local reference system of ALICE
-  // X = eta, Y = phi, Z = r (a constant for the EMCAL)
   
   lpos.SetX(fLocPos.X()) ;
   lpos.SetY(fLocPos.Y()) ;
@@ -680,11 +664,9 @@ void AliEMCALRecPoint::GetGlobalPosition(TVector3 & gpos) const
 {
   // returns the position of the cluster in the global reference system of ALICE
   // These are now the Cartesian X, Y and Z
-
-  //AliEMCALGeometry * geom = (AliEMCALGetter::Instance())->EMCALGeometry();
-  AliEMCALGeometry * geom =  AliEMCALGeometry::GetInstance();
-  Int_t absid = geom->TowerIndexFromEtaPhi(fLocPos.X(), TMath::RadToDeg()*fLocPos.Y());
-  geom->XYZFromIndex(absid, gpos);
+  geom = AliEMCALGeometry::GetInstance();
+  //  cout<<" geom "<<geom<<endl;
+  geom->GetGlobal(fLocPos, gpos, fSuperModuleNumber);
 }
 
 //____________________________________________________________________________
@@ -826,7 +808,7 @@ Float_t AliEMCALRecPoint::ThetaToEta(Float_t arg) const
 void AliEMCALRecPoint::Print(Option_t *) const
 {
   // Print the list of digits belonging to the cluster
-  
+  return;
   TString message ; 
   message  = "AliPHOSEmcRecPoint:\n" ;
   message +=  " digits # = " ; 
@@ -835,21 +817,29 @@ void AliEMCALRecPoint::Print(Option_t *) const
   Int_t iDigit;
   for(iDigit=0; iDigit<fMulDigit; iDigit++)
     printf(" %d ", fDigitsList[iDigit] ) ;  
-  
+  printf("\n");
+
   Info("Print", " Energies = ") ;
   for(iDigit=0; iDigit<fMulDigit; iDigit++) 
     printf(" %f ", fEnergyList[iDigit] ) ;
-  printf("\n") ; 
-   Info("Print", " Primaries  ") ;
+  printf("\n");
+
+  Info("Print", "\n Abs Ids = ") ;
+  for(iDigit=0; iDigit<fMulDigit; iDigit++) 
+    printf(" %i ", fAbsIdList[iDigit] ) ;
+  printf("\n");
+
+  Info("Print", " Primaries  ") ;
   for(iDigit = 0;iDigit < fMulTrack; iDigit++)
     printf(" %d ", fTracksList[iDigit]) ;
-  printf("\n") ;       
+
+  printf("\n Local x %6.2f y %7.2f z %7.1f \n", fLocPos[0], fLocPos[1], fLocPos[2]);
+
   message  = "       Multiplicity    = %d" ;
   message += "       Cluster Energy  = %f" ; 
   message += "       Core energy     = %f" ; 
   message += "       Core radius     = %f" ; 
   message += "       Number of primaries %d" ; 
   message += "       Stored at position %d" ; 
   Info("Print", message.Data(), fMulDigit, fAmp, fCoreEnergy, fCoreRadius, fMulTrack, GetIndexInList() ) ;  
 }
index a79d16b..6ca191e 100644 (file)
@@ -39,15 +39,16 @@ class AliEMCALRecPoint : public AliRecPoint {
 
   virtual void    EvalAll(Float_t logWeight, TClonesArray * digits);
   virtual void    EvalLocalPosition(Float_t logWeight, TClonesArray * digits) ;
+  //  void            EvalLocalPositionSimple(TClonesArray *digits); // ??
   virtual void    EvalPrimaries(TClonesArray * digits) ;
   virtual void    EvalParents(TClonesArray * digits) ;
 
   // virtual void    GetGlobalPosition(TVector3 & gpos, TMatrix & /*gmat*/) const; // return global position in ALICE
   virtual void    GetGlobalPosition(TVector3 & gpos) const; // return global position (x, y, z) in ALICE
-  virtual void    GetLocalPosition(TVector3 & lpos) const; // return local position (eta, phi, r) in EMCAL
+  virtual void    GetLocalPosition(TVector3 & lpos) const;  // return local position  (x, y, z) in EMCAL SM
   virtual Int_t * GetPrimaries(Int_t & number) const {number = fMulTrack ; 
                                                       return fTracksList ; }
-    virtual Int_t * GetParents(Int_t & number) const {number = fMulParent ; 
+  virtual Int_t * GetParents(Int_t & number) const {number = fMulParent ; 
                                                       return fParentsList ; }
   Float_t         GetCoreEnergy()const {return fCoreEnergy ;}
   virtual Float_t GetDispersion()const {return fDispersion ;}
@@ -58,7 +59,9 @@ class AliEMCALRecPoint : public AliRecPoint {
   Int_t       GetMaximumMultiplicity() const {return fMaxDigit ;}  // gets the maximum number of digits allowed
   Int_t       GetMultiplicity(void) const { return fMulDigit ; }   // gets the number of digits making this recpoint
   Int_t       GetMultiplicityAtLevel(Float_t level) const ;  // computes multiplicity of digits with 
-                                                                   // energy above relative level
+  Int_t *     GetAbsId() const {return fAbsIdList;}
+  Int_t       GetAbsId(int i) const {if(i>=0 && i<fMulDigit)return fAbsIdList[i]; else return -1;}
+  // energy above relative level
   virtual Int_t GetNumberOfLocalMax(AliEMCALDigit **  maxAt, Float_t * maxAtEnergy,
                                     Float_t locMaxCut,TClonesArray * digits ) const ; 
                                                                    // searches for the local maxima 
@@ -90,13 +93,15 @@ protected:
          Float_t fLambda[2] ;        // shower ellipse axes
          Float_t fDispersion ;       // shower dispersion
          Float_t *fEnergyList ;      //[fMulDigit] energy of digits
+          Int_t   *fAbsIdList;        //[fMulDigit] absId  of digits
          Float_t fTime ;             // Time of the digit with maximal energy deposition
          Float_t fCoreRadius;        // The radius in which the core energy is evaluated
           Int_t fMulParent;           // Multiplicity of the parents
           Int_t fMaxParent;           // Maximum number of parents allowed
           Int_t * fParentsList;       // [fMulParent] list of the parents of the digits
-
-  ClassDef(AliEMCALRecPoint,6) // RecPoint for EMCAL (Base Class)
+          Int_t   fSuperModuleNumber; //
+    
+  ClassDef(AliEMCALRecPoint,7) // RecPoint for EMCAL (Base Class)
  
 };
 
index 2d7fb96..676cd71 100644 (file)
@@ -92,6 +92,8 @@ AliEMCALSDigitizer::AliEMCALSDigitizer(const char * alirunFileName,
   InitParameters() ; 
   fDefaultInit = kFALSE ; 
   fHists = 0;
+  fControlHists = 1;
+  if(fControlHists) BookControlHists(1);
 }
 
 
@@ -155,6 +157,7 @@ void AliEMCALSDigitizer::InitParameters()
 
  // threshold for deposit energy of hit
   fECPrimThreshold  = 0.; // 24-nov-04 - was 1.e-6;
+  Print("");
 }
 
 //____________________________________________________________________________
@@ -323,7 +326,7 @@ void AliEMCALSDigitizer::Print1(Option_t * option)
 void AliEMCALSDigitizer::Print(Option_t *option) const
 { 
   // Prints parameters of SDigitizer
-  printf("Print: \n------------------- %s -------------\n", GetName() ) ; 
+  printf("Print: \n------------------- %s ------------- option %s\n", GetName() , option) ; 
   printf("   fInit                                 %i\n", int(fInit));
   printf("   fFirstEvent                           %i\n", fFirstEvent);
   printf("   fLastEvent                            %i\n", fLastEvent);
@@ -405,6 +408,7 @@ TList *AliEMCALSDigitizer::BookControlHists(int var)
   gROOT->cd();
   const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance() ;
   if(var>=1){
+    printf(" AliEMCALSDigitizer::BookControlHists() in action \n");
     new TH1F("HSDigiN",  "#EMCAL  sdigits ", 1001, -0.5, 1000.5);
     new TH1F("HSDigiSumEnergy","Sum.EMCAL energy", 1000, 0.0, 100.);
     new TH1F("HSDigiAmp",  "EMCAL sdigits amplitude", 1000, 0., 2.e+9);
@@ -412,8 +416,10 @@ TList *AliEMCALSDigitizer::BookControlHists(int var)
     new TH1F("HSDigiAbsId","EMCAL absID for sdigits",
     geom->GetNCells(), 0.5, Double_t(geom->GetNCells())+0.5);
   }
+
   fHists = AliEMCALHistoUtilities::MoveHistsToList("EmcalSDigiControlHists", kFALSE);
   fHists = 0;
+
   return fHists;
 }
 
index d988b3e..1ca7780 100644 (file)
 
 #include "AliEMCALShishKebabTrd1Module.h"
 #include <assert.h>
-#include <TMath.h>
 #include "AliEMCALGeometry.h"
 
+#include "Riostream.h"
+#include <TMath.h>
+
 ClassImp(AliEMCALShishKebabTrd1Module)
 
   AliEMCALGeometry *AliEMCALShishKebabTrd1Module::fgGeometry=0; 
@@ -29,13 +31,15 @@ ClassImp(AliEMCALShishKebabTrd1Module)
   Double_t AliEMCALShishKebabTrd1Module::fga2=0.; 
   Double_t AliEMCALShishKebabTrd1Module::fgb=0.; 
   Double_t AliEMCALShishKebabTrd1Module::fgr=0.; 
-  Double_t AliEMCALShishKebabTrd1Module::fgangle=0.;   // one degrre 
+  Double_t AliEMCALShishKebabTrd1Module::fgangle=0.;   // around one degree 
   Double_t AliEMCALShishKebabTrd1Module::fgtanBetta=0; //
 
-AliEMCALShishKebabTrd1Module::AliEMCALShishKebabTrd1Module(double theta) : TNamed()
+AliEMCALShishKebabTrd1Module::AliEMCALShishKebabTrd1Module(double theta, AliEMCALGeometry *g) : TNamed()
 { // theta in radians ; first object shold be with theta=pi/2.
+  cout<< " theta " << theta << " geometry " << g << endl;  
   fTheta = theta;
   if(fgGeometry==0) {
+    fgGeometry = g;
     if(GetParameters()) {
       DefineFirstModule();
     }
@@ -76,15 +80,15 @@ void AliEMCALShishKebabTrd1Module::Init(double A, double B)
   fB = yA - fA*xA;
 
   DefineName(fTheta);
-  // Centers of modules
+  // Centers of module
   Double_t kk1 = (fga+fga2)/(2.*4.); // kk1=kk2 
 
   Double_t xk1 = fOK.X() - kk1*TMath::Sin(fTheta);
-  Double_t yk1 = fOK.Y() + kk1*TMath::Cos(fTheta);
+  Double_t yk1 = fOK.Y() + kk1*TMath::Cos(fTheta) - fgr;
   fOK1.Set(xk1,yk1);
 
   Double_t xk2 = fOK.X() + kk1*TMath::Sin(fTheta);
-  Double_t yk2 = fOK.Y() - kk1*TMath::Cos(fTheta);
+  Double_t yk2 = fOK.Y() - kk1*TMath::Cos(fTheta) - fgr;
   fOK2.Set(xk2,yk2);
 }
 
@@ -99,8 +103,8 @@ void AliEMCALShishKebabTrd1Module::DefineFirstModule()
   fB = yA - fA*xA;
 
   Double_t kk1 = (fga+fga2)/(2.*4.); // kk1=kk2 
-  fOK1.Set(fOK.X() - kk1, fOK.Y());
-  fOK2.Set(fOK.X() + kk1, fOK.Y());
+  fOK1.Set(fOK.X() - kk1, fOK.Y()-fgr);
+  fOK2.Set(fOK.X() + kk1, fOK.Y()-fgr);
 
   TObject::SetUniqueID(1); //
 }
@@ -115,7 +119,7 @@ void AliEMCALShishKebabTrd1Module::DefineName(double theta)
 
 Bool_t AliEMCALShishKebabTrd1Module::GetParameters()
 {
-  fgGeometry = AliEMCALGeometry::GetInstance();
+  if(!fgGeometry) fgGeometry = AliEMCALGeometry::GetInstance();
   TString sn(fgGeometry->GetName()); // 2-Feb-05
   sn.ToUpper();
   if(!fgGeometry) {
@@ -138,18 +142,18 @@ Bool_t AliEMCALShishKebabTrd1Module::GetParameters()
 void AliEMCALShishKebabTrd1Module::PrintShish(int pri) const
 {
   if(pri>=0) {
-    Info("PrintShish()", "\n a %7.3f:%7.3f | b %7.2f | r %7.2f \n TRD1 angle %7.6f(%5.2f) | tanBetta %7.6f", 
+    Info("\n PrintShish()", "\n a %7.3f:%7.3f | b %7.2f | r %7.2f \n TRD1 angle %7.6f(%5.2f) | tanBetta %7.6f", 
     fga, fga2, fgb, fgr, fgangle, fgangle*TMath::RadToDeg(), fgtanBetta);
     printf(" fTheta %f : %5.2f : cos(theta) %f\n", 
     fTheta, GetThetaInDegree(),TMath::Cos(fTheta)); 
-    if(pri>0) {
+    if(pri>=1) {
       printf("\n%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());
-      fOK.Dump();
       printf(" fOK  : X %8.3f: Y %8.3f \n",  fOK.X(), fOK.Y());
-      printf(" fOK1 : X %8.3f: Y %8.3f \n", fOK1.X(), fOK1.Y());
-      printf(" fOK2 : X %8.3f: Y %8.3f \n", fOK2.X(), fOK2.Y());
+      printf(" fOK1 : X %8.3f: Y %8.3f (in SM, ieta=2)\n", fOK1.X(), fOK1.Y());
+      printf(" fOK2 : X %8.3f: Y %8.3f (in SM, ieta=1)\n\n", fOK2.X(), fOK2.Y());
+      fOK.Dump();
     }
   }
 }
index a97015c..7b96e7f 100644 (file)
@@ -16,7 +16,7 @@ class AliEMCALGeometry;
 
 class AliEMCALShishKebabTrd1Module : public TNamed {
  public:
-  AliEMCALShishKebabTrd1Module(double theta=TMath::Pi()/2.);
+  AliEMCALShishKebabTrd1Module(double theta=TMath::Pi()/2., AliEMCALGeometry *g=0);
   AliEMCALShishKebabTrd1Module(AliEMCALShishKebabTrd1Module &leftNeighbor);
   void Init(double A, double B);
 
@@ -36,9 +36,15 @@ class AliEMCALShishKebabTrd1Module : public TNamed {
   Double_t  GetA() {return fA;}
   Double_t  GetB() {return fB;}
   //  Additional offline staff 
-  TVector2& GetCenterOfCell(Int_t ieta)
-  { if(ieta<=1) return fOK1;
-    else        return fOK2;}
+  TVector2& GetCenterOfCellInLocalCoordinateofSM(Int_t ieta)
+  { if(ieta<=1) return fOK2;
+    else        return fOK1;}
+  void GetCenterOfCellInLocalCoordinateofSM(Int_t ieta, Double_t &xr, Double_t &zr)
+  { 
+    if(ieta<=1) {xr = fOK2.Y(); zr = fOK2.X();
+    } else      {xr = fOK1.Y(); zr = fOK1.X();
+    }
+  }
   // 
   Double_t GetTanBetta() {return fgtanBetta;}
   Double_t Getb()        {return fgb;}
@@ -61,15 +67,13 @@ class AliEMCALShishKebabTrd1Module : public TNamed {
   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
-  // position of towers with differents ieta (1 or 2) -  4-nov-04
-  TVector2 fOK1;
-  TVector2 fOK2;
+  // position of towers(cells) with differents ieta (1 or 2) in local coordinate of SM
+  // Nov 04,2004; Feb 19,2006 
+  TVector2 fOK1; // ieta=2
+  TVector2 fOK2; // ieta=1
 
  public:
   ClassDef(AliEMCALShishKebabTrd1Module,0) // Turned Shish-Kebab module 
 };
 
 #endif
-/* To do
- 1. Insert position the center of towers - 2 additional TVector2
- */
index e8c9350..677f3c5 100644 (file)
@@ -67,8 +67,9 @@ AliEMCALv0::AliEMCALv0(const char *name, const char *title):
   AliEMCAL(name,title)
 {
   // ctor : title is used to identify the layout
-  GetGeometry() ; 
-  fShishKebabModules = 0;
+  AliEMCALGeometry *geom = GetGeometry() ; 
+  //geom->CreateListOfTrd1Modules(); 
+  fShishKebabModules = geom->GetShishKebabTrd1Modules(); 
 }
 
 //______________________________________________________________________
@@ -78,7 +79,8 @@ void AliEMCALv0::BuildGeometry()
 
     const Int_t kColorArm1   = kBlue ;
 
-    AliEMCALGeometry * geom = GetGeometry() ; 
+    AliEMCALGeometry * geom = GetGeometry();
+
     TString gn(geom->GetName());
     gn.ToUpper(); 
 
@@ -319,7 +321,7 @@ void AliEMCALv0::CreateGeometry()
         printf("CreateGeometry: XU0 = %f, %f\n", envelopA[5], envelopA[6]); 
       }
     // Position the EMCAL Mother Volume (XEN1) in Alice (ALIC)  
-      gMC->Gspos("XEN1", 1, "ALIC", 0.0, 0.0, 0.0, idrotm, "ONLY") ;
+      gMC->Gspos(geom->GetNameOfEMCALEnvelope(), 1, "ALIC", 0.0, 0.0, 0.0, idrotm, "ONLY") ;
     }
 
     if(gn.Contains("SHISH")){
@@ -368,7 +370,7 @@ void AliEMCALv0::CreateShishKebabGeometry()
   //  idAL = 1602;
   Double_t par[10], xpos=0., ypos=0., zpos=0.;
 
-  CreateSmod("XEN1");        // 18-may-05 
+  CreateSmod(g->GetNameOfEMCALEnvelope());
 
   CreateEmod("SMOD","EMOD"); // 18-may-05
 
@@ -537,8 +539,12 @@ void AliEMCALv0::CreateSmod(const char* mother)
       par[2] += 0.4;      // for 27 division
     } else if(gn.Contains("TRD")) {
       par[2]  = 350./2.; // 11-oct-04 - for 26 division
+      if(gn.Contains("TRD1")) {
+        printf(" par[0] %7.2f (old) \n",  par[0]);  
+        Float_t *parSM = g->GetSuperModulesPars(); 
+        for(int i=0; i<3; i++) par[i] = parSM[i];
+      }
     }
-  //  par[2] = g->GetXYModuleSize()*g->GetNZ()/2.; 
     gMC->Gsvolu("SMOD", "BOX", idtmed[idAIR], par, 3);
     printf("tmed %i | dx %7.2f dy %7.2f dz %7.2f (SMOD, BOX)\n", 
     idtmed[idAIR], par[0],par[1],par[2]);
@@ -727,18 +733,7 @@ void AliEMCALv0::CreateEmod(const char* mother, const char* child)
   } else if(gn.Contains("TRD")) { // 30-sep-04; 27-jan-05 - as for TRD1 as for TRD2
     // X->Z(0, 0); Y->Y(90, 90); Z->X(90, 0)
     AliEMCALShishKebabTrd1Module *mod=0, *mTmp; // current module
-    if(fShishKebabModules == 0) {
-      fShishKebabModules = new TList;
-      for(int iz=0; iz<g->GetNZ(); iz++) { // 27-may-05; g->GetNZ() -> 26
-        if(iz==0) { 
-          mod  = new AliEMCALShishKebabTrd1Module();
-        } else {
-          mTmp  = new AliEMCALShishKebabTrd1Module(*mod);
-          mod   = mTmp;
-        }
-      fShishKebabModules->Add(mod);
-      }
-    }
+
     for(int iz=0; iz<g->GetNZ(); iz++) {
       Double_t  angle=90., phiOK=0;
       if(gn.Contains("TRD1")) {