"Version
authorpavlinov <pavlinov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 1 Aug 2005 16:11:19 +0000 (16:11 +0000)
committerpavlinov <pavlinov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 1 Aug 2005 16:11:19 +0000 (16:11 +0000)
35 files changed:
EMCAL/AliEMCAL.cxx
EMCAL/AliEMCAL.h
EMCAL/AliEMCALClusterizerv1.cxx
EMCAL/AliEMCALClusterizerv1.h
EMCAL/AliEMCALDigit.cxx
EMCAL/AliEMCALDigitizer.cxx
EMCAL/AliEMCALDigitizer.h
EMCAL/AliEMCALGeneratorFactory.cxx [new file with mode: 0644]
EMCAL/AliEMCALGeneratorFactory.h [new file with mode: 0644]
EMCAL/AliEMCALGeometry.cxx
EMCAL/AliEMCALGeometry.h
EMCAL/AliEMCALGeometryOfflineTrd1.cxx [new file with mode: 0644]
EMCAL/AliEMCALGeometryOfflineTrd1.h [new file with mode: 0644]
EMCAL/AliEMCALJetMicroDst.cxx
EMCAL/AliEMCALJetMicroDst.h
EMCAL/AliEMCALRecPoint.cxx
EMCAL/AliEMCALReconstruct.C
EMCAL/AliEMCALReconstructor.cxx
EMCAL/AliEMCALReconstructor.h
EMCAL/AliEMCALSDigitizer.cxx
EMCAL/AliEMCALSDigitizer.h
EMCAL/AliEMCALShishKebabModule.cxx [new file with mode: 0644]
EMCAL/AliEMCALShishKebabModule.h [new file with mode: 0644]
EMCAL/AliEMCALShishKebabTrd1Module.cxx [new file with mode: 0644]
EMCAL/AliEMCALShishKebabTrd1Module.h [new file with mode: 0644]
EMCAL/AliEMCALWsuCosmicRaySetUp.cxx [new file with mode: 0644]
EMCAL/AliEMCALWsuCosmicRaySetUp.h [new file with mode: 0644]
EMCAL/AliEMCALv0.cxx
EMCAL/AliEMCALv0.h
EMCAL/AliEMCALv1.cxx
EMCAL/AliEMCALv1.h
EMCAL/AliEMCALv2.cxx [new file with mode: 0644]
EMCAL/AliEMCALv2.h [new file with mode: 0644]
EMCAL/EMCALLinkDef.h
EMCAL/libEMCAL.pkg

index 2822ed5..044c664 100644 (file)
@@ -121,6 +121,12 @@ void AliEMCAL::CreateMaterials()
   AliMaterial(3, "Al$", 26.98, 13., 2.7, 8.9, 999., 0, 0) ;
   // ---         Absorption length is ignored ^
 
+  // 25-aug-04 by PAI - see  PMD/AliPMDv0.cxx for STEEL definition
+  Float_t asteel[4] = { 55.847,51.9961,58.6934,28.0855 };
+  Float_t zsteel[4] = { 26.,24.,28.,14. };
+  Float_t wsteel[4] = { .715,.18,.1,.005 };
+  AliMixture(4, "STAINLESS STEEL$", asteel, zsteel, 7.88, 4, wsteel);
+
   // DEFINITION OF THE TRACKING MEDIA
 
   // for EMCAL: idtmed[1599->1698] equivalent to fIdtmed[0->100]
@@ -145,6 +151,9 @@ void AliEMCAL::CreateMaterials()
   AliMedium(3, "Al parts     $", 3, 0,
              isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.001, 0.001, 0, 0) ;
 
+  // 25-aug-04 by PAI : see  PMD/AliPMDv0.cxx for STEEL definition                 -> idtmed[1603]
+  AliMedium(4, "S steel$", 4, 0, 
+             isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.001, 0.001, 0, 0) ;
 
 // --- Set decent energy thresholds for gamma and electron tracking
 
@@ -170,12 +179,21 @@ void AliEMCAL::CreateMaterials()
   gMC->Gstpar(idtmed[1601],"CUTELE",0.001) ;
   gMC->Gstpar(idtmed[1601],"BCUTE",0.0001) ;
 
+  // the same parameters as for Pb - 10-sep-04
+  gMC->Gstpar(idtmed[1603],"CUTGAM",0.00008) ;
+  gMC->Gstpar(idtmed[1603],"CUTELE",0.001) ;
+  gMC->Gstpar(idtmed[1603],"BCUTE", 0.0001); 
+  // --- Generate explicitly delta rays in Lead ---
+  gMC->Gstpar(idtmed[1603], "LOSS",3.);
+  gMC->Gstpar(idtmed[1603], "DRAY",1.);
+  gMC->Gstpar(idtmed[1603], "DCUTE",0.00001);
+  gMC->Gstpar(idtmed[1603], "DCUTM",0.00001);
+
   //set constants for Birk's Law implentation
   fBirkC0 =  1;
   fBirkC1 =  0.013/dP;
   fBirkC2 =  9.6e-6/(dP * dP);
 
-
 }
       
 //____________________________________________________________________________
@@ -365,27 +383,21 @@ void AliEMCAL::SetTreeAddress()
 { 
   // Linking Hits in Tree to Hits array
   TBranch *branch;
-  char branchname[20];
-  sprintf(branchname,"%s",GetName());
-  
+  //  char branchname[20];
+  //  sprintf(branchname,"%s",GetName());
   // Branch address for hit tree
   TTree *treeH = TreeH();
   if (treeH) {
-    branch = treeH->GetBranch(branchname);
-    if (branch) 
-      { 
+    //    treeH->Print();
+    branch = treeH->GetBranch(GetName());
+    if (branch) { 
        if (fHits == 0x0) 
-         fHits= new TClonesArray("AliEMCALHit",1000);
+       fHits= new TClonesArray("AliEMCALHit",1000);
        branch->SetAddress(&fHits);
-      }
-    else
-      {
-       Warning("SetTreeAddress","<%s> Failed",GetName());
-      }
+    } else {
+      Warning("SetTreeAddress","<%s> Failed",GetName());
+    }
+  } else {
+    //    Warning("SetTreeAddress"," no treeH ");
   }
 }
-
-
-
-
-
index 13a103e..36093e4 100644 (file)
@@ -65,7 +65,7 @@ class AliEMCAL : public AliDetector {
   virtual const TString Version() const {return TString(" ") ; }   
   AliEMCAL & operator = (const AliEMCAL & /*rvalue*/)  {
     Fatal("operator =", "not implemented") ;  return *this ; }
+
 protected:
   
   friend class AliEMCALGetter;
@@ -76,6 +76,7 @@ protected:
   Int_t fBirkC0;    // constants for Birk's Law implementation
   Double_t fBirkC1; // constants for Birk's Law implementation
   Double_t fBirkC2; // constants for Birk's Law implementation
+
   static Double_t fgCapa ;              // capacitor of the preamplifier for the raw RO signal
   Double_t fHighCharge ;                // high charge (to convert energy to charge) for the raw RO signal
   Double_t fHighGain ;                  // high gain for the raw RO signal
index 951c97a..bfb8d83 100644 (file)
 #include "AliEMCALGeometry.h"
 
 ClassImp(AliEMCALClusterizerv1)
-  
+
+Int_t addOn[20][60][60]; 
+
 //____________________________________________________________________________
   AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer()
 {
   // default ctor (to be used mainly by Streamer)
   
   InitParameters() ; 
-  fDefaultInit = kTRUE ; 
+  fDefaultInit = kTRUE ;
+  cout<<"file to read 1"<<endl;
+  ReadFile();
+  cout<<"file read 1"<<endl;
 }
 
 //____________________________________________________________________________
@@ -88,6 +93,16 @@ 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;
+   }
+   }
+ }
+  cout<<"file to read 2"<<endl;
+  ReadFile();
+  cout<<"file read 2"<<endl;
 
 }
 
@@ -129,14 +144,13 @@ void AliEMCALClusterizerv1::Exec(Option_t * option)
   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
 
   if (fLastEvent == -1) 
-    fLastEvent = gime->MaxEvent() - 1;
-
+    //fLastEvent = gime->MaxEvent() - 1;
+    fLastEvent = 1000 - 1;
   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 
   Int_t ievent ;
   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
     gime->Event(ievent,"D") ;
-
     GetCalibrationParameters() ;
 
     fNumberOfECAClusters = 0;
@@ -162,6 +176,9 @@ void AliEMCALClusterizerv1::Exec(Option_t * option)
     printf("Exec took %f seconds for Clusterizing %f seconds per event", 
         gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nEvents ) ;
   }
+
+  SaveHists();
+
 }
 
 //____________________________________________________________________________
@@ -267,6 +284,7 @@ void AliEMCALClusterizerv1::GetCalibrationParameters()
 
   fADCchannelECA   = dig->GetECAchannel() ;
   fADCpedestalECA  = dig->GetECApedestal();
+  cout<<"ChannelECA, peds "<<fADCchannelECA<<" "<<fADCpedestalECA<<endl;
 }
 
 //____________________________________________________________________________
@@ -280,13 +298,16 @@ void AliEMCALClusterizerv1::Init()
     gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data());
 
   AliEMCALGeometry * geom = gime->EMCALGeometry() ;
+  cout<<"gime,geom "<<gime<<" "<<geom<<endl;
 
-  fNTowers = geom->GetNZ() *  geom->GetNPhi() ;
+//Sub  fNTowers = geom->GetNZ() *  geom->GetNPhi() ;
+  fNTowers =400;
   if(!gMinuit) 
     gMinuit = new TMinuit(100) ;
-
- if ( !gime->Clusterizer() ) 
-    gime->PostClusterizer(this); 
+ //Sub if ( !gime->Clusterizer() ) 
+ //Sub   gime->PostClusterizer(this); 
+ BookHists();
+  cout<<"hists booked "<<endl;
 }
 
 //____________________________________________________________________________
@@ -294,8 +315,10 @@ void AliEMCALClusterizerv1::InitParameters()
 { 
   // Initializes the parameters for the Clusterizer
   fNumberOfECAClusters = 0;
-  fECAClusteringThreshold   = 1.0;  // must be adjusted according to the noise level set by digitizer
-  fECALocMaxCut = 0.03 ;
+
+  fECAClusteringThreshold   = 0.5;  // value obtained from Alexei
+  fECALocMaxCut = 0.03;
+
   fECAW0     = 4.5 ;
   fTimeGate = 1.e-8 ; 
   fToUnfold = kFALSE ;
@@ -318,10 +341,32 @@ Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d
   Int_t rv = 0 ; 
 
   Int_t relid1[2] ; 
-  geom->AbsToRelNumbering(d1->GetId(), relid1) ; 
-
-  Int_t relid2[2] ; 
-  geom->AbsToRelNumbering(d2->GetId(), relid2) ; 
+ //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(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(nTower1,nIphi1,nIeta1, iphi1,ieta1);
+   Int_t relid2[2] ; 
+   relid2[0]=ieta1;
+   relid2[1]=iphi1;
+
+  //Sub  geom->AbsToRelNumbering(d2->GetId(), relid2) ; 
   
 
   Int_t rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ;  
@@ -336,7 +381,7 @@ Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d
   }
  
   if (gDebug == 2 ) 
-    printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d ", 
+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]) ;   
   
@@ -388,8 +433,8 @@ void AliEMCALClusterizerv1::WriteRecPoints()
 
   branchECA->Fill() ;
 
-  gime->WriteRecPoints("OVERWRITE");
-  gime->WriteClusterizer("OVERWRITE");
+//Sub  gime->WriteRecPoints("OVERWRITE");
+//Sub  gime->WriteClusterizer("OVERWRITE");
 }
 
 //____________________________________________________________________________
@@ -412,18 +457,38 @@ void AliEMCALClusterizerv1::MakeClusters()
   // Clusterization starts    
   TIter nextdigit(digitsC) ;
   AliEMCALDigit * digit;
-  
+
+//just for hist
+/*
+  while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
+         digitAmp->Fill(digit->GetAmp());
+  }
+  */
+/////////// 
+
   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
     AliEMCALRecPoint * clu = 0 ; 
     
     TArrayI clusterECAdigitslist(50000);   
    
-    if (gDebug == 2) { 
-      printf("MakeClusters: id = %d, ene = %f , thre = %f", digit->GetId(),Calibrate(digit->GetAmp()),
-         fECAClusteringThreshold) ;  
-    }
-
-    if ( geom->IsInECA(digit->GetId()) && (Calibrate(digit->GetAmp()) > fECAClusteringThreshold  ) ){
+//Sub    if (gDebug == 2) { 
+   //  printf("MakeClusters: id = %d, ene = %f , thre = %f \n", digit->GetId(),Calibrate(digit->GetAmp()),
+//         fECAClusteringThreshold) ;  
+//    }
+////////////////////////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(nTower,nIphi,nIeta, iphi,ieta);
+
+         /////////////////////////////////// 
+
+//cout<<ieta<<" "<<iphi<<endl;
+
+//    if ( geom->IsInECA(digit->GetId()) && (Calibrate(digit->GetAmp()) > fECAClusteringThreshold  ) ){
+    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()) 
@@ -432,12 +497,13 @@ void AliEMCALClusterizerv1::MakeClusters()
       aECARecPoints->AddAt(rp, fNumberOfECAClusters) ;
       clu = dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(fNumberOfECAClusters)) ; 
       fNumberOfECAClusters++ ; 
-      clu->AddDigit(*digit, Calibrate(digit->GetAmp())) ; 
+      clu->AddDigit(*digit, Calibrate(digit->GetAmp()+addOn[nSupMod-1][ieta-1][iphi-1])) ; 
       clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ;     
       iDigitInECACluster++ ; 
+//    cout<<" 1st setp:cluno, digNo "<<fNumberOfECAClusters<<" "<<iDigitInECACluster<<endl;
       digitsC->Remove(digit) ; 
       if (gDebug == 2 ) 
-       printf("MakeClusters: OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp()), fECAClusteringThreshold) ;  
+       printf("MakeClusters: OK id = %d, ene = %f , thre = %f \n", digit->GetId(),Calibrate(digit->GetAmp()), fECAClusteringThreshold) ;  
       nextdigit.Reset() ;
       
       AliEMCALDigit * digitN ; 
@@ -448,32 +514,46 @@ 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
-          if(Calibrate(digitN->GetAmp()) < fMinECut  )  digitsC->Remove(digitN);
+         ////////////////////
+          geom->GetCellIndex(digitN->GetId(), nSupMod,nTower,nIphi,nIeta);
+          geom->GetCellPhiEtaIndexInSModule(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 
-           clu->AddDigit(*digitN, Calibrate( digitN->GetAmp()) ) ;
+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
-           goto endofloop1;   
+//         break ;
+//          case 2 :   // too far from each other
+//Subh     goto endofloop1;   
+//         cout<<"earlier go to end of loop"<<endl;   
          } // switch
+    //cout<<"in nextDigit loop "<<fNumberOfECAClusters<<" "<<iDigitInECACluster<<endl;
        } // while digitN
        
-      endofloop1: ;
+//Sub      endofloop1: ;
        nextdigit.Reset() ; 
       } // loop over ECA cluster
     } // energy threshold
-    else if(Calibrate(digit->GetAmp()) < fMinECut  ){
+    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);
     }
+    //cout<<"after endofloop: cluno, digNo "<<fNumberOfECAClusters<<endl;
   } // while digit  
   delete digitsC ;
+cout<<"total no of clusters "<<fNumberOfECAClusters<<"from "<<digits->GetEntriesFast()<<" digits"<<endl; 
 }
 
 //____________________________________________________________________________
@@ -572,11 +652,15 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
     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") ;      
-    
+   Float_t maxE=0; 
+   Float_t maxL1=0; 
+   Float_t maxL2=0; 
+   Float_t maxDis=0; 
+
     for (index = 0 ; index < aECARecPoints->GetEntries() ; index++) {
       AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(aECARecPoints->At(index)) ; 
       TVector3  globalpos;  
-      rp->GetGlobalPosition(globalpos);
+      //rp->GetGlobalPosition(globalpos);
       TVector3  localpos;  
       rp->GetLocalPosition(localpos);
       Float_t lambda[2]; 
@@ -588,10 +672,73 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
             rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
             globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), 
             rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ; 
+  /////////////
+      if(rp->GetEnergy()>maxE){
+             maxE=rp->GetEnergy();
+             maxL1=lambda[0];
+             maxL2=lambda[1];
+             maxDis=rp->GetDispersion();
+      }
+      pointE->Fill(rp->GetEnergy());
+      pointL1->Fill(lambda[0]);
+      pointL2->Fill(lambda[1]);
+      pointDis->Fill(rp->GetDispersion());
+      pointMult->Fill(rp->GetMultiplicity());
+      ///////////// 
       for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
        printf("%d ", primaries[iprimary] ) ; 
       } 
     }
+
+      MaxE->Fill(maxE);
+      MaxL1->Fill(maxL1);
+      MaxL2->Fill(maxL2);
+      MaxDis->Fill(maxDis);
+
+
     printf("\n-----------------------------------------------------------------------\n");
   }
 }
+  void AliEMCALClusterizerv1::BookHists()
+{
+       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.);
+       pointDis = new TH1F("pointDis","point Dis", 1000, 0.0, 3.);
+       pointMult = new TH1F("pointMult","point Mult", 100, 0.0, 100.);
+       digitAmp = new TH1F("digitAmp","Digit Amplitude", 2000, 0.0, 5000.);
+       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.);
+}
+void AliEMCALClusterizerv1::SaveHists()
+{
+ 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();
+}
+
+void AliEMCALClusterizerv1::ReadFile()
+{
+  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;
+       }
+}
+
index e2648d9..49c260d 100644 (file)
@@ -24,6 +24,7 @@
 // --- AliRoot header files ---
 
 #include "AliEMCALClusterizer.h"
+#include "TH1F.h"
 class AliEMCALRecPoint ; 
 class AliEMCALDigit ;
 class AliEMCALDigitizer ;
@@ -69,12 +70,27 @@ public:
                                             // Chi^2 of the fit. Should be static to be passes to MINUIT
   void Unload() ; 
   virtual const char * Version() const { return "clu-v1" ; }  
-
+  void BookHists();
+  void SaveHists();
 protected:
 
   void           WriteRecPoints() ;
   virtual void   MakeClusters( ) ;            
-  
+///////////////////// 
+   TH1F* pointE;
+   TH1F* pointL1;
+   TH1F* pointL2;
+   TH1F* pointDis;
+   TH1F* pointMult;
+   TH1F* digitAmp;
+   TH1F* MaxE;
+   TH1F* MaxL1;
+   TH1F* MaxL2;
+   TH1F* MaxDis;
+///////////////////////
+
+
 private:
 
   const TString BranchName() const ; 
@@ -90,7 +106,7 @@ private:
                               AliEMCALDigit ** /*maxAt*/,
                               Float_t * /*maxAtEnergy*/ ) const; //Unfolds cluster using TMinuit package
   void           PrintRecPoints(Option_t * option) ;
-
+  TFile* recofile;
 private:
 
   Bool_t  fDefaultInit;              //! Says if the task was created by defaut ctor (only parameters are initialized)
@@ -111,6 +127,7 @@ 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 a5664c3..602eab8 100644 (file)
@@ -153,7 +153,7 @@ Int_t AliEMCALDigit::Compare(const TObject * obj) const
 
 //____________________________________________________________________________
 Float_t AliEMCALDigit::GetEta() const
-{
+{ // should be change in EMCALGeometry - 19-nov-04
   Float_t eta=-10., phi=-10.;
   Int_t id = GetId();
   const AliEMCALGeometry *g = AliEMCALGetter::Instance()->EMCALGeometry();
@@ -163,7 +163,7 @@ Float_t AliEMCALDigit::GetEta() const
 
 //____________________________________________________________________________
 Float_t AliEMCALDigit::GetPhi() const
-{
+{ // should be change in EMCALGeometry - 19-nov-04
   Float_t eta=-10., phi=-10.;
   Int_t id = GetId();
   const AliEMCALGeometry *g = AliEMCALGetter::Instance()->EMCALGeometry();
index ac37421..54167f0 100644 (file)
 // Modif: 
 //  August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
 //                           of new  IO (à la PHOS)
-///_________________________________________________________________________________
+//  November 2003 Aleksei Pavlinov : adopted for Shish-Kebab geometry 
+//_________________________________________________________________________________
 
 // --- ROOT system ---
+#include "TROOT.h"
 #include "TTree.h"
 #include "TSystem.h"
 #include "TBenchmark.h"
-
-// --- Standard library ---
+#include "TList.h"
+#include "TH1.h"
+#include "TBrowser.h"
+#include "TObjectTable.h"
 
 // --- AliRoot header files ---
 #include "AliRunDigitizer.h"
@@ -73,6 +77,7 @@
 #include "AliEMCALSDigitizer.h"
 #include "AliEMCALGeometry.h"
 #include "AliEMCALTick.h"
+#include "AliEMCALJetMicroDst.h"
 
 ClassImp(AliEMCALDigitizer)
 
@@ -160,20 +165,35 @@ void AliEMCALDigitizer::Digitize(Int_t event)
   // after that adds contributions from SDigits. This design 
   // helps to avoid scanning over the list of digits to add 
   // contribution of any new SDigit.
+  static int isTrd1Geom = -1; // -1 - mean undefined 
+  static int nEMC=0; //max number of digits possible
 
+  //<<<<<<< AliEMCALDigitizer.cxx
+  int deb=0;
+  AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName) ; 
+  /* =======
   AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle()) ; 
+  >>>>>>> 1.59 */
   Int_t ReadEvent = event ; 
+  // fManager is data member from AliDigitizer
   if (fManager) 
     ReadEvent = dynamic_cast<AliStream*>(fManager->GetInputStream(0))->GetCurrentEventNumber() ; 
-  Info("Digitize", "Adding event %d from input stream 0 %s %s", ReadEvent, GetTitle(), fEventFolderName.Data()) ; 
+  if(deb>0) Info("Digitize", "Adding event %d from input stream 0 %s %s", 
+  ReadEvent, GetTitle(), fEventFolderName.Data()) ; 
   gime->Event(ReadEvent, "S") ;
   TClonesArray * digits = gime->Digits() ; 
   digits->Clear() ;
 
-  const AliEMCALGeometry *geom = gime->EMCALGeometry() ; 
+  const AliEMCALGeometry *geom = gime->EMCALGeometry() ;
+  if(isTrd1Geom < 0) { 
+    TString ng(geom->GetName());
+    isTrd1Geom = 0;
+    if(ng.Contains("SHISH") &&  ng.Contains("TRD1")) isTrd1Geom = 1;
 
-  Int_t nEMC = geom->GetNPhi()*geom->GetNZ(); //max number of digits possible
-  
+    if(isTrd1Geom == 0) nEMC = geom->GetNPhi()*geom->GetNZ();
+    else                nEMC = geom->GetNCells();
+    printf(" nEMC %i (number cells in EMCAL) | %s | isTrd1Geom %i\n", nEMC, geom->GetName(), isTrd1Geom);
+  }
   Int_t absID ;
 
   digits->Expand(nEMC) ;
@@ -200,7 +220,7 @@ void AliEMCALDigitizer::Digitize(Int_t event)
     gime->Event(ReadEvent,"S");
     sdigArray->AddAt(gime->SDigits(), i) ;
   }
-  
   //Find the first tower with signal
   Int_t nextSig = nEMC + 1 ; 
   TClonesArray * sdigits ;  
@@ -211,7 +231,9 @@ void AliEMCALDigitizer::Digitize(Int_t event)
     Int_t curNext = dynamic_cast<AliEMCALDigit *>(sdigits->At(0))->GetId() ;
      if(curNext < nextSig) 
        nextSig = curNext ;
+     if(deb>0) printf("input %i : #sdigits %i \n", i, sdigits->GetEntriesFast());
   }
+  if(deb>0) printf("FIRST tower with signal %i \n", nextSig);
 
   TArrayI index(fInput) ;
   index.Reset() ;  //Set all indexes to zero
@@ -221,7 +243,7 @@ void AliEMCALDigitizer::Digitize(Int_t event)
 
   TClonesArray * ticks = new TClonesArray("AliEMCALTick",1000) ;
 
-  //Put Noise contribution
+  //Put Noise contributiony
   for(absID = 1; absID <= nEMC; absID++){
     Float_t amp = 0 ;
     // amplitude set to zero, noise will be added later
@@ -265,7 +287,7 @@ void AliEMCALDigitizer::Digitize(Int_t event)
 
          index[i]++ ;
          if( dynamic_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
-           curSDigit = dynamic_cast<AliEMCALDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;            
+           curSDigit = dynamic_cast<AliEMCALDigit*>(dynamic_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
          else
            curSDigit = 0 ;
        }
@@ -292,7 +314,8 @@ void AliEMCALDigitizer::Digitize(Int_t event)
     // add the noise now
     amp += TMath::Abs(gRandom->Gaus(0., fPinNoise)) ;
     digit->SetAmp(sDigitizer->Digitize(amp)) ;  
-  }
+    if(deb>=10) printf(" absID %5i amp %f nextSig %5i\n", absID, amp, nextSig);
+  } // for(absID = 1; absID <= nEMC; absID++)
   
   ticks->Delete() ;
   delete ticks ;
@@ -314,13 +337,20 @@ void AliEMCALDigitizer::Digitize(Int_t event)
   Int_t ndigits = digits->GetEntriesFast() ; 
   digits->Expand(ndigits) ;
   
-  //Set indexes in list of digits
+  //Set indexes in list of digits and fill hists.
+  sv::FillH1(fHists, 0, Double_t(ndigits));
+  Float_t energy=0., esum=0.;
   for (i = 0 ; i < ndigits ; i++) { 
     digit = dynamic_cast<AliEMCALDigit *>( digits->At(i) ) ; 
     digit->SetIndexInList(i) ; 
-    Float_t energy = sDigitizer->Calibrate(digit->GetAmp()) ;
-    digit->SetAmp(DigitizeEnergy(energy) ) ;
+    energy = sDigitizer->Calibrate(digit->GetAmp()) ;
+    esum += energy;
+    digit->SetAmp(DigitizeEnergy(energy) ) ; // for what ??
+    sv::FillH1(fHists, 2, double(digit->GetAmp()));
+    sv::FillH1(fHists, 3, double(energy));
+    sv::FillH1(fHists, 4, double(digit->GetId()));
   }
+  sv::FillH1(fHists, 1, esum);
 }
 
 //____________________________________________________________________________
@@ -364,14 +394,17 @@ void AliEMCALDigitizer::Exec(Option_t *option)
   if (fLastEvent == -1) 
     fLastEvent = gime->MaxEvent() - 1 ;
   else if (fManager) 
-    fLastEvent = fFirstEvent ; 
+    fLastEvent = fFirstEvent ; // what is this ??
 
   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
-  
-  Int_t ievent ;
-
-  for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
-  
+  Int_t ievent, nfr=50;
+
+  fLastEvent = TMath::Min(fLastEvent, gime->MaxEvent());
+  printf("AliEMCALDigitizer::Exec : option: %s | %i -> %i events : Max events %i \n", 
+  option, fFirstEvent, fLastEvent,  gime->MaxEvent());
+  for (ievent = fFirstEvent; ievent < fLastEvent; ievent++) {
+    if(ievent%nfr==0 || ievent==fLastEvent-1);
+    printf(" processed event %i\n", ievent);
     gime->Event(ievent,"S") ; 
 
     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
@@ -380,11 +413,13 @@ void AliEMCALDigitizer::Exec(Option_t *option)
 
     if(strstr(option,"deb"))
       PrintDigits(option);
+    if(strstr(option,"table")) gObjectTable->Print();
 
     //increment the total number of Digits per run 
     fDigitsInRun += gime->Digits()->GetEntriesFast() ;  
   }
-  
+  //  gime->WriteDigitizer("OVERWRITE");
+
   gime->EmcalLoader()->CleanDigitizer() ;
 
   if(strstr(option,"tim")){
@@ -459,21 +494,24 @@ Bool_t AliEMCALDigitizer::Init()
 
 //____________________________________________________________________________ 
 void AliEMCALDigitizer::InitParameters()
-{
-  fMeanPhotonElectron = 18200 ; // electrons per GeV
-  fPinNoise           = 0.003 ; // noise equivalent GeV (random choice)
+{ // Tune parameters - 24-nov-04
+
+  fMeanPhotonElectron = 3300 ; // electrons per GeV 
+  fPinNoise           = 0.004; 
   if (fPinNoise == 0. ) 
     Warning("InitParameters", "No noise added\n") ; 
-  fDigitThreshold     = fPinNoise * 3; //2 sigma
-  fTimeResolution     = 1.0e-9 ;
+  fDigitThreshold     = fPinNoise * 3; // 3 * sigma
+  fTimeResolution     = 0.3e-9 ; // 300 psc
   fTimeSignalLength   = 1.0e-9 ;
 
-  fADCchannelEC    = 0.00305;                     // width of one ADC channel in GeV - HG fix so that we see 200 GeV gammas
-  fADCpedestalEC   = 0.005 ;                       // GeV
-  fNADCEC          = (Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC
+  fADCchannelEC    = 0.00305; // 200./65536 - width of one ADC channel in GeV
+  fADCpedestalEC   = 0.009 ;  // GeV
+  fNADCEC          = (Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC - 65536
 
-  fTimeThreshold      = 0.001*10000000 ; //Means 1 MeV in terms of SDigits amplitude
+  fTimeThreshold      = 0.001*10000000 ; // Means 1 MeV in terms of SDigits amplitude ??
+  // hists. for control; no hists on default
+  fControlHists = 0;
+  fHists        = 0;
 }
 
 //__________________________________________________________________
@@ -527,7 +565,13 @@ void AliEMCALDigitizer::MixWith(TString alirunFileName, TString eventFolderName)
     fEventNames[fInput]     = eventFolderName ;
     fInput++ ;
 }  
+
+void AliEMCALDigitizer::Print1(Option_t * option)
+{ // 19-nov-04 - just for convinience
+  Print(); 
+  PrintDigits(option);
+}
+
 //__________________________________________________________________
 void AliEMCALDigitizer::Print()const 
 {
@@ -569,11 +613,11 @@ void AliEMCALDigitizer::Print()const
 void AliEMCALDigitizer::PrintDigits(Option_t * option){
     
   AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName) ; 
-  TClonesArray * digits = gime->Digits() ;
+  TClonesArray * digits  = gime->Digits() ;
+  TClonesArray * sdigits = gime->SDigits() ;
   
-  printf("PrintDigits: %d", digits->GetEntriesFast()) ; 
-  printf("\nevent %d", gAlice->GetEvNumber()) ;
-  printf("\n       Number of entries in Digits list %d", digits->GetEntriesFast() )  ;  
+  printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ; 
+  printf("\n event %d", gAlice->GetEvNumber()) ;
   
   if(strstr(option,"all")){  
     //loop over digits
@@ -591,6 +635,7 @@ void AliEMCALDigitizer::PrintDigits(Option_t * option){
       }
     }   
   }
+  printf("\n");
 }
 
 //__________________________________________________________________
@@ -629,7 +674,12 @@ void AliEMCALDigitizer::WriteDigits()
   //      and branch "AliEMCALDigitizer", with the same title to keep all the parameters
   //      and names of files, from which digits are made.
 
+  //<<<<<<< AliEMCALDigitizer.cxx
+  static Int_t writeSdigitizer=0;
+  AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName) ; 
+  /* =======
   AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle()) ; 
+  >>>>>>> 1.59*/
   const TClonesArray * digits = gime->Digits() ; 
   TTree * treeD = gime->TreeD(); 
 
@@ -640,9 +690,41 @@ void AliEMCALDigitizer::WriteDigits()
   digitsBranch->Fill() ;
   
   gime->WriteDigits("OVERWRITE");
-  gime->WriteDigitizer("OVERWRITE");
+  if(writeSdigitizer==0) {
+    gime->WriteDigitizer("OVERWRITE");
+    writeSdigitizer = 1;
+  }
 
   Unload() ; 
 
 }
 
+void AliEMCALDigitizer::Browse(TBrowser* b)
+{
+  if(fHists) b->Add(fHists);
+  TTask::Browse(b);
+}
+
+TList *AliEMCALDigitizer::BookControlHists(const int var)
+{ // 22-nov-04
+  Info("BookControlHists"," started ");
+  gROOT->cd();
+  AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName);
+  const AliEMCALGeometry *geom = gime->EMCALGeometry() ;
+  if(var>=1){
+    new TH1F("hDigiN",  "#EMCAL digits with fAmp > fDigitThreshold", 
+    fNADCEC+1, -0.5, Double_t(fNADCEC));
+    new TH1F("HDigiSumEnergy","Sum.EMCAL energy from digi", 1000, 0.0, 200.);
+    new TH1F("hDigiAmp",  "EMCAL digital amplitude", fNADCEC+1, -0.5, Double_t(fNADCEC));
+    new TH1F("hDigiEnergy","EMCAL cell energy", 2000, 0.0, 200.);
+    new TH1F("hDigiAbsId","EMCAL absId cells with fAmp > fDigitThreshold ",
+    geom->GetNCells(), 0.5, Double_t(geom->GetNCells())+0.5);
+  }
+  fHists = sv::MoveHistsToList("EmcalDigiControlHists", kFALSE);
+  return fHists;
+}
+
+void AliEMCALDigitizer::SaveHists(const char* name, Bool_t kSingleKey, const char* opt)
+{
+  sv::SaveListOfHists(fHists, name, kSingleKey, opt); 
+}
index 970b82c..61e2874 100644 (file)
@@ -10,7 +10,8 @@
 //                  
 //*-- Author: Sahal Yacoob (LBL)
 // based on : AliPHOSDigit
-// July 2003 Yves Schutz : NewIO 
+// July     2003 Yves Schutz : NewIO 
+// November 2003 Aleksei Pavlinov : Shish-Kebab geometry 
 //_________________________________________________________________________ 
 
 
@@ -18,6 +19,8 @@
 #include "TObjString.h"
 class TArrayI ;
 class TClonesArray ; 
+class TList;
+class TBrowser;
 
 // --- Standard library ---
 
@@ -54,7 +57,8 @@ public:
   Int_t   GetDigitsInRun()  const { return fDigitsInRun; } 
   void  MixWith(TString alirunFileName, 
                TString eventFolderName = AliConfig::GetDefaultEventFolderName()) ; // Add another one file to mix
-  void  Print()const ;
+  void  Print() const ;
+  void  Print1(Option_t * option); // *MENU*
  
   AliEMCALDigitizer & operator = (const AliEMCALDigitizer & /*rvalue*/)  {
     // assignement operator requested by coding convention but not needed
@@ -62,6 +66,15 @@ public:
    return *this ; 
   }
 
+  virtual void Browse(TBrowser* b);
+  // hists
+  void   SetControlHists(const Int_t var=0) {fControlHists=var;}
+  Int_t  GetControlHist() const {return fControlHists;}
+  TList *GetListOfHists() {return fHists;}
+  TList* BookControlHists(const int var=0);
+  void   SaveHists(const char* name="RF/TRD1/Digitizations/DigiVar?",
+  Bool_t kSingleKey=kTRUE, const char* opt="RECREATE"); // *MENU*
+
 private:
 
   Bool_t  Init();                   
@@ -100,9 +113,11 @@ private:
   TString fEventFolderName;         // skowron: name of EFN to read data from in stand alone mode
   Int_t   fFirstEvent;        // first event to process
   Int_t   fLastEvent;         // last  event to process
+  // Control hists
+  Int_t   fControlHists;          //!
+  TList  *fHists;                 //!
 
   ClassDef(AliEMCALDigitizer,5)  // description 
-
 };
 
 
diff --git a/EMCAL/AliEMCALGeneratorFactory.cxx b/EMCAL/AliEMCALGeneratorFactory.cxx
new file mode 100644 (file)
index 0000000..012fbbb
--- /dev/null
@@ -0,0 +1,669 @@
+/**************************************************************************
+ * Copyright(c) 1998-2002, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/*
+$Log$
+*/
+#include "assert.h"
+#include "TString.h"
+#include "AliEMCALGeneratorFactory.h"
+#include "AliGenerator.h"
+#include "AliGenFixed.h"
+#include "AliGenBox.h"
+#include "AliGenHIJINGpara.h"
+#include "AliGenHIJINGparaBa.h"
+#include "AliGenHijing.h"
+#include "AliGenCocktail.h"
+#include "AliGenPythia.h"
+
+//*-- Authors: Aleksei Pavlinov (WSU)
+//*   Initial variant is in Config.C for EMCAL production.
+ClassImp(AliEMCALGeneratorFactory)
+
+AliEMCALGeneratorFactory::AliEMCALGeneratorFactory(PprRunFact_t run, PprRadFact_t rad) 
+{
+    fGenerator = fBgGenerator = fSignalGenerator = 0;
+    fRunType   = run;
+    fRadiation = rad;
+    fMomentum  = 0;
+
+    Int_t isw = 3;
+    if (rad == kNoGluonRadiation) isw = 0;    
+    Float_t thmin=0, thmax=0;
+
+    AliGenFixed *genGG=0;
+    AliGenBox   *genB=0;
+    AliGenHIJINGparaBa *genHijParaB=0, *bg=0;
+    AliGenHIJINGpara   *genHijPara=0;
+    AliGenHijing       *genHij=0;
+    AliGenCocktail     *genCoct=0;
+    AliGenPythia       *genPy=0, *jets=0;
+    //    AliPythia          *aliPy = 0;
+
+    fComment = new TString;
+
+    switch (run) {
+    case kGammaGun:
+        genGG = new AliGenFixed(1);
+       genGG->SetMomentum(100.);
+       genGG->SetPhi(240.);
+       genGG->SetTheta(91.);
+       genGG->SetPart(kGamma);
+       fGenerator = (AliGenerator*)genGG;
+     break;
+
+    case kGammaBox:
+       genB = new AliGenBox(100);
+       genB->SetMomentumRange(100., 100.1);
+       genB->SetPhiRange(0,360);
+       genB->SetThetaRange(45., 135.);
+       genB->SetPart(kGamma);
+        fGenerator = (AliGenerator*)genB;
+       break;
+       
+    case kTest50:
+       genHijParaB = new AliGenHIJINGparaBa(50);
+       genHijParaB->SetMomentumRange(0, 999999.);
+       genHijParaB->SetPhiRange(-180., 180.);
+       // Set pseudorapidity range from -8 to 8.
+       thmin = EtaToTheta(8);   // theta min. <---> eta max
+       thmax = EtaToTheta(-8);  // theta max. <---> eta min 
+       genHijParaB->SetThetaRange(thmin,thmax);
+        fGenerator = (AliGenerator*)genHijParaB;
+       break;
+
+    case kParam_8000:
+       //coment= fComment.Append(":HIJINGparam N=8000");
+       genHijPara = new AliGenHIJINGpara(86030);
+       genHijPara->SetMomentumRange(0, 999999.);
+       genHijPara->SetPhiRange(-180., 180.);
+       // Set pseudorapidity range from -8 to 8.
+       thmin = EtaToTheta(8);   // theta min. <---> eta max
+       thmax = EtaToTheta(-8);  // theta max. <---> eta min 
+       genHijPara->SetThetaRange(thmin,thmax);
+        fGenerator = (AliGenerator*)genHijPara;
+       break;
+    case kParam_4000:
+       genHijPara = new AliGenHIJINGpara(43015);
+       genHijPara->SetMomentumRange(0, 999999.);
+       genHijPara->SetPhiRange(-180., 180.);
+       // Set pseudorapidity range from -8 to 8.
+       thmin = EtaToTheta(8);   // theta min. <---> eta max
+       thmax = EtaToTheta(-8);  // theta max. <---> eta min 
+       genHijPara->SetThetaRange(thmin,thmax);
+        fGenerator = (AliGenerator*)genHijPara;
+       break;
+    case kParam_2000:
+        (*fComment) = "HIJINGparam N=2000";
+       genHijPara = new AliGenHIJINGpara(21507);
+       genHijPara->SetMomentumRange(0, 999999.);
+       genHijPara->SetPhiRange(-180., 180.);
+       // Set pseudorapidity range from -8 to 8.
+       thmin = EtaToTheta(8);   // theta min. <---> eta max
+       thmax = EtaToTheta(-8);  // theta max. <---> eta min 
+       genHijPara->SetThetaRange(thmin,thmax);
+        fGenerator = (AliGenerator*)genHijPara;
+       break;
+
+    case kParam_8000_Ecal:
+       genHijParaB = new AliGenHIJINGparaBa(82534);
+       genHijParaB->SetMomentumRange(0, 999999.);
+       genHijParaB->SetPhiRange(-180., 180.);
+       // Set pseudorapidity range from -8 to 8.
+       thmin = EtaToTheta( 5);   // theta min. <---> eta max
+       thmax = EtaToTheta(-5);  // theta max. <---> eta min 
+       genHijParaB->SetThetaRange(thmin,thmax);
+        fGenerator = (AliGenerator*)genHijParaB;
+       break;
+
+    case kParam_4000_Ecal:
+       genHijParaB = new AliGenHIJINGparaBa(82534/2);
+       genHijParaB->SetMomentumRange(0, 999999.);
+       genHijParaB->SetPhiRange(-180., 180.);
+       // Set pseudorapidity range from -8 to 8.
+       thmin = EtaToTheta( 5);   // theta min. <---> eta max
+       thmax = EtaToTheta(-5);  // theta max. <---> eta min 
+       genHijParaB->SetThetaRange(thmin,thmax);
+        fGenerator = (AliGenerator*)genHijParaB;
+       break;
+//
+//  Hijing Central
+//
+    case kHijing_cent1:
+       genHij = HijingStandard();
+// impact parameter range
+       genHij->SetImpactParameterRange(0., 5.);
+        fGenerator = (AliGenerator*)genHij;
+       break;
+    case kHijing_cent2:
+       genHij = HijingStandard();
+// impact parameter range
+       genHij->SetImpactParameterRange(0., 2.);
+        fGenerator = (AliGenerator*)genHij;
+       break;
+//
+// Hijing Peripheral 
+//
+    case kHijing_per1:
+       genHij = HijingStandard();
+// impact parameter range
+       genHij->SetImpactParameterRange(5., 8.6);
+        fGenerator = (AliGenerator*)genHij;
+       break;
+    case kHijing_per2:
+       //coment= comment.Append("HIJING per2");
+       genHij = HijingStandard();
+// impact parameter range
+       genHij->SetImpactParameterRange(8.6, 11.2);
+        fGenerator = (AliGenerator*)genHij;
+       break;
+    case kHijing_per3:
+       genHij = HijingStandard();
+// impact parameter range
+       genHij->SetImpactParameterRange(11.2, 13.2);
+        fGenerator = (AliGenerator*)genHij;
+       break;
+    case kHijing_per4:
+       genHij = HijingStandard();
+// impact parameter range
+       genHij->SetImpactParameterRange(13.2, 15.);
+        fGenerator = (AliGenerator*)genHij;
+       break;
+    case kHijing_per5:
+       //coment= comment.Append("HIJING per5");
+       genHij = HijingStandard();
+// impact parameter range
+       genHij->SetImpactParameterRange(15., 100.);
+        fGenerator = (AliGenerator*)genHij;
+       break;
+//
+//  Jet-Jet
+//
+    case kHijing_jj25:
+       //coment= comment.Append("HIJING Jet 25 GeV");
+       genHij = HijingStandard();
+// impact parameter range
+       genHij->SetImpactParameterRange(0., 5.);
+       // trigger
+       genHij->SetTrigger(1);
+       genHij->SetPtJet(25.);
+       genHij->SetSimpleJets(1);
+       genHij->SetRadiation(isw);
+       genHij->SetJetEtaRange(-0.3,0.3);
+       genHij->SetJetPhiRange(15.,105.);   
+        fGenerator = (AliGenerator*)genHij;
+       break;
+
+    case kHijing_jj50:
+       //coment= comment.Append("HIJING Jet 50 GeV");
+       genHij = HijingStandard();
+// impact parameter range
+       genHij->SetImpactParameterRange(0., 5.);
+       // trigger
+       genHij->SetTrigger(1);
+       genHij->SetPtJet(50.);
+       genHij->SetSimpleJets(1);
+       genHij->SetRadiation(isw);
+       genHij->SetJetEtaRange(-0.3,0.3);
+       genHij->SetJetPhiRange(15.,105.);   
+        fGenerator = (AliGenerator*)genHij;
+       break;
+
+    case kHijing_jj75:
+       //coment= comment.Append("HIJING Jet 75 GeV");
+       genHij = HijingStandard();
+// impact parameter range
+       genHij->SetImpactParameterRange(0., 5.);
+       // trigger
+       genHij->SetTrigger(1);
+       genHij->SetPtJet(75.);
+       genHij->SetSimpleJets(1);
+       genHij->SetRadiation(isw);
+       genHij->SetJetEtaRange(-0.3,0.3);
+       genHij->SetJetPhiRange(15.,105.);   
+        fGenerator = (AliGenerator*)genHij;
+       break;
+
+    case kHijing_jj100:
+       //coment= comment.Append("HIJING Jet 100 GeV");
+       genHij = HijingStandard();
+// impact parameter range
+       genHij->SetImpactParameterRange(0., 5.);
+       // trigger
+       genHij->SetTrigger(1);
+       genHij->SetPtJet(100.);
+       genHij->SetSimpleJets(1);
+       genHij->SetRadiation(isw);
+       genHij->SetJetEtaRange(-0.3,0.3);
+       genHij->SetJetPhiRange(15.,105.);   
+        fGenerator = (AliGenerator*)genHij;
+       break;
+
+    case kHijing_jj125:
+       //coment= comment.Append("HIJING Jet 125 GeV");
+       genHij = HijingStandard();
+// impact parameter range
+       genHij->SetImpactParameterRange(0., 5.);
+       // trigger
+       genHij->SetTrigger(1);
+       genHij->SetPtJet(125.);
+       genHij->SetSimpleJets(1);
+       genHij->SetRadiation(isw);
+       genHij->SetJetEtaRange(-0.3,0.3);
+       genHij->SetJetPhiRange(15.,105.);   
+        fGenerator = (AliGenerator*)genHij;
+       break;
+//
+// Gamma-Jet
+//
+    case kHijing_gj25:
+       //coment= comment.Append("HIJING Gamma 25 GeV");
+       genHij = HijingStandard();
+// impact parameter range
+       genHij->SetImpactParameterRange(0., 5.);
+       // trigger
+       genHij->SetTrigger(2);
+       genHij->SetPtJet(25.);
+       genHij->SetSimpleJets(1);
+       genHij->SetRadiation(isw);
+       genHij->SetJetEtaRange(-0.3,0.3);
+       genHij->SetJetPhiRange(15.,105.);   
+        fGenerator = (AliGenerator*)genHij;
+       break;
+
+    case kHijing_gj50:
+       //coment= comment.Append("HIJING Gamma 50 GeV");
+       genHij = HijingStandard();
+// impact parameter range
+       genHij->SetImpactParameterRange(0., 5.);
+       // trigger
+       genHij->SetTrigger(2);
+       genHij->SetPtJet(50.);
+       genHij->SetSimpleJets(1);
+       genHij->SetRadiation(isw);
+       genHij->SetJetEtaRange(-0.3,0.3);
+       genHij->SetJetPhiRange(15.,105.);   
+        fGenerator = (AliGenerator*)genHij;
+       break;
+
+    case kHijing_gj75:
+       //coment= comment.Append("HIJING Gamma 75 GeV");
+       genHij = HijingStandard();
+// impact parameter range
+       genHij->SetImpactParameterRange(0., 5.);
+       // trigger
+       genHij->SetTrigger(2);
+       genHij->SetPtJet(75.);
+       genHij->SetSimpleJets(1);
+       genHij->SetRadiation(isw);
+       genHij->SetJetEtaRange(-0.3,0.3);
+       genHij->SetJetPhiRange(15.,105.);   
+        fGenerator = (AliGenerator*)genHij;
+       break;
+
+    case kHijing_gj100:
+       //coment= comment.Append("HIJING Gamma 100 GeV");
+       genHij = HijingStandard();
+// impact parameter range
+       genHij->SetImpactParameterRange(0., 5.);
+       // trigger
+       genHij->SetTrigger(2);
+       genHij->SetPtJet(100.);
+       genHij->SetSimpleJets(1);
+       genHij->SetRadiation(isw);
+       genHij->SetJetEtaRange(-0.3,0.3);
+       genHij->SetJetPhiRange(15.,105.);   
+        fGenerator = (AliGenerator*)genHij;
+       break;
+
+    case kHijing_gj125:
+       //coment= comment.Append("HIJING Gamma 125 GeV");
+       genHij = HijingStandard();
+// impact parameter range
+       genHij->SetImpactParameterRange(0., 5.);
+       // trigger
+       genHij->SetTrigger(2);
+       genHij->SetPtJet(125.);
+       genHij->SetSimpleJets(1);
+       genHij->SetRadiation(isw);
+       genHij->SetJetEtaRange(-0.3,0.3);
+       genHij->SetJetPhiRange(15.,105.);   
+        fGenerator = (AliGenerator*)genHij;
+       break;
+    case kJetPlusBg:
+       genCoct = new AliGenCocktail();
+       genCoct->SetMomentumRange(0, 999999.);
+       genCoct->SetPhiRange(-180., 180.);
+       // Set pseudorapidity range from -8 to 8.
+       thmin = EtaToTheta( 5.);   // theta min. <---> eta max
+       thmax = EtaToTheta(-5.);  // theta max. <---> eta min 
+       genCoct->SetThetaRange(thmin,thmax);
+
+//
+//      Underlying Event
+//
+//     AliGenHIJINGparaBa *bg = new AliGenHIJINGparaBa(82534);
+       bg = new AliGenHIJINGparaBa(10);
+        fBgGenerator = (AliGenerator*)bg;
+//
+//      Jets from Pythia
+//
+       jets = new AliGenPythia(-1);
+        fSignalGenerator = (AliGenerator*)jets; 
+//   Centre of mass energy 
+       jets->SetEnergyCMS(5500.);
+//   Process type
+       jets->SetProcess(kPyJets);
+//   final state kinematic cuts
+       jets->SetJetEtaRange(-0.3, 0.3);
+       jets->SetJetPhiRange(15., 105.);
+//   Structure function
+       jets->SetStrucFunc(kGRVLO98);
+//   
+//   Pt transfer of the hard scattering
+       jets->SetPtHard(100.,100.1);
+//   Decay type (semielectronic, semimuonic, nodecay)
+       jets->SetForceDecay(kAll);
+//
+//      Add all to cockail ...    
+//
+       genCoct->AddGenerator(jets,"Jets",1);
+       genCoct->AddGenerator(bg,"Underlying Event", 1);
+        fGenerator = (AliGenerator*)genCoct;
+
+     break;    
+    case kGammaPlusBg:
+       genCoct = new AliGenCocktail();
+       genCoct->SetMomentumRange(0, 999999.);
+       genCoct->SetPhiRange(-180., 180.);
+       // Set pseudorapidity range from -8 to 8.
+       thmin = EtaToTheta( 5.);   // theta min. <---> eta max
+       thmax = EtaToTheta(-5.);  // theta max. <---> eta min 
+       genCoct->SetThetaRange(thmin,thmax);
+//
+//      Underlying Event
+//
+       bg = new AliGenHIJINGparaBa(82534);
+        fBgGenerator = (AliGenerator*)bg;
+//
+//      Jets from Pythia
+//
+       jets = new AliGenPythia(-1);
+        fSignalGenerator = (AliGenerator*)jets; 
+//   Centre of mass energy 
+       jets->SetEnergyCMS(5500.);
+//   Process type
+       jets->SetProcess(kPyDirectGamma);
+//   final state kinematic cuts
+       jets->SetJetEtaRange(-0.3, 0.3);
+       jets->SetJetPhiRange(15., 105.);
+       jets->SetGammaEtaRange(-0.12, 0.12);
+       jets->SetGammaPhiRange(220., 320.);
+//   Structure function
+       jets->SetStrucFunc(kGRVLO98);
+//   
+//   Pt transfer of the hard scattering
+       jets->SetPtHard(100.,100.1);
+//   Decay type (semielectronic, semimuonic, nodecay)
+       jets->SetForceDecay(kAll);
+//
+//      Add all to cockail ...    
+//
+       genCoct->AddGenerator(jets,"Jets",1);
+       genCoct->AddGenerator(bg,"Underlying Event", 1);
+        fGenerator = (AliGenerator*)genCoct;
+
+       break;
+    case kJets_50:
+//  50 GeV Jets
+        genPy = PythiaJets(50.);
+        fGenerator = (AliGenerator*)genPy;
+        break; 
+    case kJets_75:
+//  75 GeV Jets
+        genPy = PythiaJets(75.);
+        fGenerator = (AliGenerator*)genPy;
+        break; 
+    case kJets_100:
+//  100 GeV Jets  
+       genPy = PythiaJets(100.);
+        fGenerator = (AliGenerator*)genPy;
+       break; 
+    case kJets_200:
+//  200 GeV Jets
+        genPy = PythiaJets(200.);
+        fGenerator = (AliGenerator*)genPy;
+        break;
+
+    case kJets_100RadOn:
+//  100 GeV Jets with radiation on - 22-mar-2002 
+//  See AliPythia.cxx for default
+       genPy = PythiaJets(100.);
+       //        genPy->SetKeyPartonJets(1);   // for jet partons
+       //        genPy->DefineParametersForPartonsJets(); 
+
+        fGenerator = (AliGenerator*)genPy;
+       break; 
+
+    case kGammaJets_50:
+//  50 GeV Jets + Gamma
+        genPy = PythiaJets(-1);
+        genPy->SetEnergyCMS(5500.);
+        genPy->SetProcess(kPyDirectGamma);
+        genPy->SetJetEtaRange(-0.3,+0.3);
+        genPy->SetJetPhiRange(15.,105.);
+        genPy->SetGammaEtaRange(-0.12, 0.12);
+        genPy->SetGammaPhiRange(220., 320.);
+        genPy->SetStrucFunc(kGRVLO98);
+        genPy->SetPtHard(50.,50.001);
+        genPy->SetForceDecay(kAll);
+        fGenerator = (AliGenerator*)genPy;
+        break;
+    case kGammaJets_75:
+//  75 GeV Jets + Gamma 
+        genPy = PythiaJets(-1);
+        genPy->SetEnergyCMS(5500.);
+        genPy->SetProcess(kPyDirectGamma);
+        genPy->SetJetEtaRange(-0.3,+0.3);
+        genPy->SetJetPhiRange(15.,105.);
+        genPy->SetGammaEtaRange(-0.12, 0.12);
+        genPy->SetGammaPhiRange(220., 320.);
+        genPy->SetStrucFunc(kGRVLO98);
+        genPy->SetPtHard(75.,75.001);
+        genPy->SetForceDecay(kAll);
+        fGenerator = (AliGenerator*)genPy;
+        break; 
+    case kGammaJets_100:
+// 100 GeV Jets + Gamma
+        genPy = PythiaJets(-1);
+        genPy->SetEnergyCMS(5500.);
+        genPy->SetProcess(kPyDirectGamma);
+        genPy->SetJetEtaRange(-0.3,+0.3);
+        genPy->SetJetPhiRange(15.,105.);
+        genPy->SetGammaEtaRange(-0.12, 0.12);
+        genPy->SetGammaPhiRange(220., 320.);
+        genPy->SetStrucFunc(kGRVLO98);
+        genPy->SetPtHard(100.,100.001);
+        genPy->SetForceDecay(kAll);
+        fGenerator = (AliGenerator*)genPy;
+        break; 
+    case kGammaJets_200:
+//  200 GeV Jets + Gamma
+        genPy = PythiaJets(-1);
+        genPy->SetEnergyCMS(5500.);
+        genPy->SetProcess(kPyDirectGamma);
+        genPy->SetJetEtaRange(-0.3,+0.3);
+        genPy->SetJetPhiRange(15.,105.);
+        genPy->SetGammaEtaRange(-0.12, 0.12);
+        genPy->SetGammaPhiRange(220., 320.);
+        genPy->SetStrucFunc(kGRVLO98);
+        genPy->SetPtHard(200.,200.001);
+        genPy->SetForceDecay(kAll);
+        fGenerator = (AliGenPythia*)genPy;
+        break; 
+    case kGammaJets_250:
+//  250 GeV Jets + Gamma
+        genPy = PythiaJets(-1);
+        genPy->SetEnergyCMS(5500.);
+        genPy->SetProcess(kPyDirectGamma);
+        genPy->SetJetEtaRange(-0.3,+0.3);
+        genPy->SetJetPhiRange(15.,105.);
+        genPy->SetGammaEtaRange(-0.12, 0.12);
+        genPy->SetGammaPhiRange(220., 320.);
+        genPy->SetStrucFunc(kGRVLO98);
+        genPy->SetPtHard(250.,250.001);
+        genPy->SetForceDecay(kAll);
+        fGenerator = (AliGenerator*)genPy;
+        break; 
+    case kGammaJets_300:
+//  300 GeV Jets + Gamma
+        genPy = PythiaJets(-1);
+        genPy->SetEnergyCMS(5500.);
+        genPy->SetProcess(kPyDirectGamma);
+        genPy->SetJetEtaRange(-0.3,+0.3);
+        genPy->SetJetPhiRange(15.,105.);
+        genPy->SetGammaEtaRange(-0.12, 0.12);
+        genPy->SetGammaPhiRange(220., 320.);
+        genPy->SetStrucFunc(kGRVLO98);
+        genPy->SetPtHard(300.,300.001);
+        genPy->SetForceDecay(kAll);
+        fGenerator = (AliGenerator*)genPy;
+        break;
+    default:
+        printf("<I> wrong parameter for generator run %i rad %i\n", run, rad);
+        assert(0);
+    }
+    if(fGenerator) fGenerator->SetPtRange(0.,1.e10); // discard the limit on pT
+
+}
+
+AliEMCALGeneratorFactory::AliEMCALGeneratorFactory(PprRunFact_t run, Float_t p) 
+{
+    fGenerator = fBgGenerator = fSignalGenerator = 0;
+    fRunType   = run;
+    fMomentum  = p;
+
+    AliGenBox   *genB=0;
+
+    switch (run) {
+    case kGammaBoxOne:
+        genB = OneParticleWithFixedEnergy(kGamma, p);  
+    break;
+    case kPi0BoxOne:
+        genB = OneParticleWithFixedEnergy(kPi0, p);  
+    break;
+    default:
+        printf("<I> wrong parameter for generator run %i \n",run);
+        assert(0);
+    }
+    if(genB) fGenerator = (AliGenerator*)genB;
+    //    if(fGenerator) fGenerator->SetPtRange(0.,1.e10); // discard the limit on pT - 23-aug-04
+}
+
+AliGenHijing* AliEMCALGeneratorFactory::HijingStandard()
+{
+    AliGenHijing *gener = new AliGenHijing(-1);
+// centre of mass energy 
+    gener->SetEnergyCMS(5500.);
+// reference frame
+    gener->SetReferenceFrame("CMS");
+// projectile
+    gener->SetProjectile("A", 208, 82);
+    gener->SetTarget    ("A", 208, 82);
+// tell hijing to keep the full parent child chain
+    gener->KeepFullEvent();
+// enable jet quenching
+    gener->SetJetQuenching(1);
+// enable shadowing
+    gener->SetShadowing(1);
+// neutral pion and heavy particle decays switched off
+    gener->SetDecaysOff(1);
+// Don't track spectators
+    gener->SetSpectators(0);
+// kinematic selection
+    gener->SetSelectAll(0);
+    return gener;
+}
+
+AliGenPythia* AliEMCALGeneratorFactory::PythiaJets(Float_t energy)
+{
+    AliGenPythia *gener = new AliGenPythia(-1);
+//   Centre of mass energy 
+    gener->SetEnergyCMS(5500.);
+//   Process type
+    gener->SetProcess(kPyJets);
+//   final state kinematic cuts
+    gener->SetJetEtaRange(-0.3, 0.3);
+    gener->SetJetPhiRange(15., 105.);
+//   Structure function
+    gener->SetStrucFunc(kGRVLO98);
+//   
+//   Pt transfer of the hard scattering
+    gener->SetPtHard(energy, energy+0.1);
+//   Decay type (semielectronic, semimuonic, nodecay)
+    gener->SetForceDecay(kAll);
+//
+    return gener;
+}
+
+
+AliGenPythia* AliEMCALGeneratorFactory::PythiaGamma(Float_t energy)
+{
+    AliGenPythia *gener = new AliGenPythia(-1);
+//   Centre of mass energy 
+    gener->SetEnergyCMS(5500.);
+//   Process type
+    gener->SetProcess(kPyDirectGamma);
+//   final state kinematic cuts
+    gener->SetJetEtaRange(-0.3, 0.3);
+    gener->SetJetPhiRange(15., 105.);
+    gener->SetGammaEtaRange(-0.12, 0.12);
+    gener->SetGammaPhiRange(220., 320.);
+//   Structure function
+    gener->SetStrucFunc(kGRVLO98);
+//   
+//   Pt transfer of the hard scattering
+    gener->SetPtHard(energy, energy+0.1);
+//   Decay type (semielectronic, semimuonic, nodecay)
+    gener->SetForceDecay(kAll);
+//
+    return gener;
+}
+
+//
+// Staff of Aleksei Pavlinov.
+//
+AliGenBox* AliEMCALGeneratorFactory::OneParticleWithFixedEnergy(Int_t type, Float_t p)
+{// one particle in EMCAL acceptance
+   Float_t thmin  = EtaToTheta(0.7), thmax  = EtaToTheta(-0.7); 
+   Float_t phimin=0, phimax=120;
+   Float_t pmin=p, pmax=p+0.01;
+
+   AliGenBox *gen = new AliGenBox(1);
+   gen->SetPart(type);
+   gen->SetNumberParticles(1); 
+   gen->SetThetaRange(thmin, thmax);
+   gen->SetPhiRange(phimin, phimax);
+   gen->SetMomentumRange(pmin, pmax);
+
+   printf("<I> AliEMCALGeneratorFactory::OneParticleWithFixedEnergy \n");
+   printf("       type of particle -> %i \n", type);
+   printf("    %6.4f <    eta   < %6.4f\n", thmin,  thmax);
+   printf("    %6.4f <    phi   < %6.4f\n", phimin, phimax);
+   printf("    %7.2f < Momentum < %7.2f\n", pmin, pmax);
+   printf("     TestBit(kPtRange) %i | TestBit(kMomentumRange) %i\n", 
+   gen->TestBit(BIT(17)), gen->TestBit(BIT(19)));
+   return gen;
+}
diff --git a/EMCAL/AliEMCALGeneratorFactory.h b/EMCAL/AliEMCALGeneratorFactory.h
new file mode 100644 (file)
index 0000000..252882d
--- /dev/null
@@ -0,0 +1,79 @@
+#ifndef ALIEMCALGENERATORFACTORY_H
+#define ALIEMCALGENERATORFACTORY_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice     */
+
+/* $Id$ */
+
+//_________________________________________________________________________
+//  Class for generator factory which used in production for EMCAL.
+//  Based on Config.C file 
+//*-- Author: Aleksei Pavlinov (WSU)
+#include <TObject.h>
+#include <TMath.h>
+#include "AliPDG.h"
+// 23-aug-04 for kGamma and kPi0
+#include <TPDGCode.h>
+
+class AliGenerator;
+class AliGenHijing;
+class AliGenPythia;
+class AliGenBox;
+class TString;
+
+enum PprRunFact_t 
+{
+    kTest50,
+    kParam_8000,   kParam_4000,  kParam_2000, 
+    kHijing_cent1, kHijing_cent2, 
+    kHijing_per1,  kHijing_per2, kHijing_per3, kHijing_per4,  kHijing_per5,
+    kHijing_jj25,  kHijing_jj50, kHijing_jj75, kHijing_jj100, kHijing_jj125, 
+    kHijing_gj25,  kHijing_gj50, kHijing_gj75, kHijing_gj100, kHijing_gj125, 
+    kJetPlusBg,    kGammaPlusBg, 
+    kParam_8000_Ecal, kParam_4000_Ecal, 
+    kJets_50,       kJets_75,      kJets_100,      kJets_200,
+    kJets_100RadOn, 
+    kGammaJets_50,  kGammaJets_75, kGammaJets_100, kGammaJets_200,
+    kGammaJets_250, kGammaJets_300,
+    kGammaGun, kGammaBox,
+    kGammaBoxOne, kPi0BoxOne 
+};
+
+enum PprRadFact_t
+{ // Concern only HIJING 
+    kGluonRadiation, kNoGluonRadiation
+};
+
+class AliEMCALGeneratorFactory : public TObject{
+
+ public:
+  explicit AliEMCALGeneratorFactory
+  (PprRunFact_t run=kJets_50, PprRadFact_t rad = kGluonRadiation);
+  explicit AliEMCALGeneratorFactory(PprRunFact_t run, Float_t p);
+  AliGenHijing* HijingStandard();
+  AliGenPythia* PythiaJets(Float_t energy);
+  AliGenPythia* PythiaGamma(Float_t energy);
+  AliGenBox*    OneParticleWithFixedEnergy(Int_t type=kGamma, Float_t p=1.0);
+
+  AliGenerator* GetGenerator() {return fGenerator;}
+  AliGenerator* GetBgGenerator() {return fBgGenerator;}
+  AliGenerator* GetSignalGenerator() {return fSignalGenerator;}
+  PprRunFact_t  GetRunType()   {return fRunType;}
+  PprRadFact_t  GetRadiation() {return fRadiation;}
+  TString*      GetComment() {return fComment;}
+  static Float_t EtaToTheta(Float_t arg) {return (180./TMath::Pi())*2.*atan(exp(-arg));}
+
+protected:
+  AliGenerator* fGenerator;        //! 
+  AliGenerator* fBgGenerator;      //! 
+  AliGenerator* fSignalGenerator;  //! 
+  PprRunFact_t  fRunType;
+  PprRadFact_t  fRadiation;
+  Float_t       fMomentum;
+  TString      *fComment;          //!
+
+  ClassDef(AliEMCALGeneratorFactory,1) // Generator Factory for EMCAL production
+
+};
+#endif // ALIEMCALGENERATORFACTORY_H
index f6f0816..a641e94 100644 (file)
@@ -26,6 +26,7 @@
 //*-- Author: Sahal Yacoob (LBL / UCT)
 //     and  : Yves Schutz (SUBATECH)
 //     and  : Jennifer Klay (LBL)
+//     SHASHLYK : Aleksei Pavlinov (WSU)
 
 // --- AliRoot header files ---
 #include <TMath.h>
@@ -41,6 +42,7 @@ ClassImp(AliEMCALGeometry)
 
 AliEMCALGeometry *AliEMCALGeometry::fgGeom = 0;
 Bool_t            AliEMCALGeometry::fgInit = kFALSE;
+TString name; // contains name of geometry
 
 //______________________________________________________________________
 AliEMCALGeometry::~AliEMCALGeometry(void){
@@ -66,9 +68,21 @@ void AliEMCALGeometry::Init(void){
   // WX inform about the composition of the EM calorimeter section: 
   //   thickness in mm of Pb radiator (W) and of scintillator (X), and number of scintillator layers (N)
   // New geometry: EMCAL_55_25
-
+  // 24-aug-04 for shish-kebab
+  // SHISH_25 or SHISH_62
   fgInit = kFALSE; // Assume failed until proven otherwise.
-  TString name(GetName()) ; 
+  name   = GetName();
+  name.ToUpper(); 
+
+  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
+  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
+
+  // geometry
   if (name == "EMCAL_55_25") {
     fECPbRadThickness  = 0.5;  // cm, Thickness of the Pb radiators
     fECScintThick      = 0.5;  // cm, Thickness of the scintillator
@@ -82,21 +96,114 @@ void AliEMCALGeometry::Init(void){
   else if( name == "G56_2_55_19" || name == "EMCAL_5655_21" || name == "G56_2_55_19_104_14"|| name == "G65_2_64_19" || name == "EMCAL_6564_21"){
     Fatal("Init", "%s is an old geometry! Please update your Config file", name.Data()) ;
   }
+  else if(name.Contains("SHISH")){
+    fNumberOfSuperModules = 12; // 12 = 6 * 2 (6 in phi, 2 in Z);
+    fSteelFrontThick = 2.54;    //  9-sep-04
+    fIPDistance      = 460.0;
+    fFrontSteelStrip = fPassiveScintThick = 0.0; // 13-may-05
+    fLateralSteelStrip = 0.025; // before MAY 2005 
+    fPhiModuleSize   = fEtaModuleSize   = 11.4;
+    fPhiTileSize = fEtaTileSize      = 5.52; // (11.4-5.52*2)/2. = 0.18 cm (wall thickness)
+    fNPhi            = 14;
+    fNZ              = 30;
+    fAlFrontThick    = fGap2Active = 0;
+    fNPHIdiv = fNETAdiv = 2;
+
+    fNECLayers       = 62;
+    fECScintThick    = fECPbRadThickness = 0.2;
+    fSampling        = 1.;  // 30-aug-04 - should be calculated
+    if(name.Contains("TWIST")) { // all about EMCAL module
+      fNZ             = 27;  // 16-sep-04
+    } else if(name.Contains("TRD")) {
+      fIPDistance      = 428.0;  //  11-may-05
+      fSteelFrontThick = 0.0;    // 3.17 -> 0.0; 28-mar-05 : no stell plate
+      fNPhi            = 12;
+      fSampling       = 12.327;
+      fPhiModuleSize = fEtaModuleSize = 12.26;
+      fNZ            = 26;     // 11-oct-04
+      fTrd1Angle     = 1.3;    // in degree
+// 18-nov-04; 1./0.08112=12.327
+// http://pdsfweb01.nersc.gov/~pavlinov/ALICE/SHISHKEBAB/RES/linearityAndResolutionForTRD1.html
+      if(name.Contains("TRD1")) {       // 30-jan-05
+       // for final design
+        if(name.Contains("MAY05") || name.Contains("WSUC")){
+          fNumberOfSuperModules = 12; // 20-may-05
+          if(name.Contains("WSUC")) fNumberOfSuperModules = 1; // 27-may-05
+          fNECLayers     = 77;       // (13-may-05 from V.Petrov)
+          fPhiModuleSize = 12.5;     // 20-may-05 - rectangular shape
+          fEtaModuleSize = 11.9;
+          fECScintThick  = fECPbRadThickness = 0.16;// (13-may-05 from V.Petrov)
+          fFrontSteelStrip   = 0.025;// 0.025cm = 0.25mm  (13-may-05 from V.Petrov)
+          fLateralSteelStrip = 0.01; // 0.01cm  = 0.1mm   (13-may-05 from V.Petrov) - was 0.025
+          fPassiveScintThick = 0.8;  // 0.8cm   = 8mm     (13-may-05 from V.Petrov)
+          fNZ                = 24;
+          fTrd1Angle         = 1.5;  // 1.3 or 1.5
+       }
+      } else if(name.Contains("TRD2")) {       // 30-jan-05
+        fSteelFrontThick = 0.0;         // 11-mar-05
+        fIPDistance+= fSteelFrontThick; // 1-feb-05 - compensate absence of steel plate
+        fTrd1Angle  = 1.64;             // 1.3->1.64
+        fTrd2AngleY = fTrd1Angle;       //  symmetric case now
+        fEmptySpace    = 0.2; // 2 mm
+        fTubsR         = fIPDistance; // 31-jan-05 - as for Fred case
+
+        fPhiModuleSize  = fTubsR*2.*TMath::Tan(fTrd2AngleY*TMath::DegToRad()/2.);
+        fPhiModuleSize -= fEmptySpace/2.; // 11-mar-05  
+        fEtaModuleSize  = fPhiModuleSize; // 20-may-05 
+        fTubsTurnAngle  = 3.;
+      }
+      fNPHIdiv = fNETAdiv  = 2;   // 13-oct-04 - division again
+      if(name.Contains("3X3")) {   // 23-nov-04
+        fNPHIdiv = fNETAdiv  = 3;
+      } else if(name.Contains("4X4")) {
+        fNPHIdiv = fNETAdiv  = 4;
+      }
+    }
+    fPhiTileSize = fPhiModuleSize/2. - fLateralSteelStrip; // 13-may-05 
+    fEtaTileSize = fEtaModuleSize/2. - fLateralSteelStrip; // 13-may-05 
+
+    if(name.Contains("25")){
+      fNECLayers     = 25;
+      fECScintThick  = fECPbRadThickness = 0.5;
+    }
+    if(name.Contains("WSUC")){ // 18-may-05 - about common structure
+      fShellThickness = 30.; // should be change 
+      fNPhi = fNZ = 4; 
+    }
+    // constant for transition absid <--> indexes
+    fNCellsInTower  = fNPHIdiv*fNETAdiv;
+    fNCellsInSupMod = fNCellsInTower*fNPhi*fNZ;
+    fNCells         = fNCellsInSupMod*fNumberOfSuperModules;
+
+    fLongModuleSize = fNECLayers*(fECScintThick + fECPbRadThickness);
+    if(name.Contains("MAY05")) fLongModuleSize += (fFrontSteelStrip + fPassiveScintThick);
+
+    // 30-sep-04
+    if(name.Contains("TRD")) {
+      f2Trd1Dx2 = fEtaModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd1Angle*TMath::DegToRad()/2.);
+      if(name.Contains("TRD2")) {  // 27-jan-05
+        f2Trd2Dy2 = fPhiModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd2AngleY*TMath::DegToRad()/2.);
+      }
+    }
+  }
   else
     Fatal("Init", "%s is an undefined geometry!", name.Data()) ; 
                 
-  // geometry
-  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
-  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
 
+  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(name.Contains("SHISH")) {
+    fShellThickness = fSteelFrontThick + fLongModuleSize;
+    if(name.Contains("TWIST")) { // 13-sep-04
+      fShellThickness  = TMath::Sqrt(fLongModuleSize*fLongModuleSize + fPhiModuleSize*fEtaModuleSize);
+      fShellThickness += fSteelFrontThick;
+    } else if(name.Contains("TRD")) { // 1-oct-04
+      fShellThickness  = TMath::Sqrt(fLongModuleSize*fLongModuleSize + f2Trd1Dx2*f2Trd1Dx2);
+      fShellThickness += fSteelFrontThick;
+    }
+  }
 
   fZLength        = 2.*ZFromEtaR(fIPDistance+fShellThickness,fArm1EtaMax); // Z coverage
   fEnvelop[0]     = fIPDistance; // mother volume inner radius
@@ -105,12 +212,41 @@ void AliEMCALGeometry::Init(void){
   
   fgInit = kTRUE; 
   
-  if (gDebug) {
-    printf("Init: geometry of EMCAL named %s is as follows:", name.Data());
+  if (kTRUE) {
+    printf("Init: geometry of EMCAL named %s is as follows:\n", name.Data());
     printf( "               ECAL      : %d x (%f mm Pb, %f mm Sc) \n", GetNECLayers(), GetECPbRadThick(), GetECScintThick() ) ; 
+    if(name.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);
+      if(name.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(" fLongModuleSize   %6.3f cm \n", fLongModuleSize);
+      printf(" #supermodule in phi direction %i \n", fNPhiSuperModule );
+    }
+    if(name.Contains("TRD")) {
+      printf(" fTrd1Angle %7.4f\n", fTrd1Angle);
+      printf(" f2Trd1Dx2  %7.4f\n",  f2Trd1Dx2);
+      if(name.Contains("TRD2")) {
+        printf(" fTrd2AngleY     %7.4f\n", fTrd2AngleY);
+        printf(" f2Trd2Dy2       %7.4f\n", f2Trd2Dy2);
+        printf(" fTubsR          %7.2f\n", fTubsR);
+        printf(" fTubsTurnAngle  %7.4f\n", fTubsTurnAngle);
+        printf(" fEmptySpace     %7.4f\n", fEmptySpace);
+      }
+    }
     printf("Granularity: %d in eta and %d in phi\n", GetNZ(), GetNPhi()) ;
-    printf("Layout: phi = (%f, %f), eta = (%f, %f), y = %f\n",  
-          GetArm1PhiMin(), GetArm1PhiMax(),GetArm1EtaMin(), GetArm1EtaMax(), GetIPDistance() ) ;    
+    printf("Layout: phi = (%7.1f, %7.1f), eta = (%5.2f, %5.2f), IP = %7.2f\n",  
+          GetArm1PhiMin(), GetArm1PhiMax(),GetArm1EtaMin(), GetArm1EtaMax(), GetIPDistance() );
   }
 }
 
@@ -431,3 +567,66 @@ Bool_t AliEMCALGeometry::IsInEMCAL(Double_t x, Double_t y, Double_t z) const {
   }
   return 0;
 }
+
+//
+// == Shish-kebab cases ==
+//
+Int_t AliEMCALGeometry::GetAbsCellId(const int nSupMod, const int nTower, const int nIphi, const int nIeta)
+{ // 27-aug-04; corr. 21-sep-04
+  static Int_t id; // have to change from 1 to fNCells
+  id  = fNCellsInSupMod*(nSupMod-1);
+  id += fNCellsInTower *(nTower-1);
+  id += fNPHIdiv *(nIphi-1);
+  id += nIeta;
+  if(id<=0 || id > fNCells) {
+    printf(" wrong numerations !!\n");
+    printf("    id      %6i(will be force to -1)\n", id);
+    printf("    fNCells %6i\n", fNCells);
+    printf("    nSupMod %6i\n", nSupMod);
+    printf("    nTower  %6i\n", nTower);
+    printf("    nIphi   %6i\n", nIphi);
+    printf("    nIeta   %6i\n", nIeta);
+    id = -1;
+  }
+  return id;
+}
+
+Bool_t  AliEMCALGeometry::CheckAbsCellId(Int_t ind)
+{ // 17-niv-04 - analog of IsInECA
+   if(name.Contains("TRD")) {
+     if(ind<=0 || ind > fNCells) return kFALSE;
+     else                        return kTRUE;
+   } else return IsInECA(ind);
+}
+
+Bool_t AliEMCALGeometry::GetCellIndex(const Int_t absId,Int_t &nSupMod,Int_t &nTower,Int_t &nIphi,Int_t &nIeta)
+{ // 21-sep-04
+  static Int_t tmp=0;
+  if(absId<=0 || absId>fNCells) {
+    Info("GetCellIndex"," wrong abs Id %i !! \n", absId); 
+    return kFALSE;
+  }
+  nSupMod = (absId-1) / fNCellsInSupMod + 1;
+  tmp     = (absId-1) % fNCellsInSupMod;
+
+  nTower  = tmp / fNCellsInTower + 1;
+  tmp     = tmp % fNCellsInTower;
+
+  nIphi     = tmp / fNPHIdiv + 1;
+  nIeta     = tmp % fNPHIdiv + 1;
+
+  return kTRUE;
+}
+
+void AliEMCALGeometry::GetCellPhiEtaIndexInSModule(const int nTower, const int nIphi, const int nIeta, 
+int &iphi, int &ieta)
+{ // don't check validity of nTower, nIphi and nIeta index
+  // have to change  - 1-nov-04 ?? 
+  static Int_t iphit, ietat;
+
+  ietat = (nTower-1)/fNPhi;
+  ieta  = ietat*fNETAdiv + nIeta; // change from 1 to fNZ*fNETAdiv
+
+  iphit = (nTower-1)%fNPhi;
+  iphi  = iphit*fNPHIdiv + nIphi;  // change from 1 to fNPhi*fNPHIdiv
+}
index 3dd0be0..cb17c68 100644 (file)
@@ -1,6 +1,6 @@
 #ifndef ALIEMCALGEOMETRY_H
 #define ALIEMCALGEOMETRY_H
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+/* Copyright(c) 1998-2004, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
 
 /* $Id$ */
@@ -78,7 +78,34 @@ public:
   Float_t GetECScintThick() const {return fECScintThick;}
   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 GetPhiModuleSize() const  {return fPhiModuleSize;}
+  Float_t GetEtaModuleSize() const  {return fEtaModuleSize;}
+  Float_t GetFrontSteelStrip() const {return fFrontSteelStrip;}
+  Float_t GetLateralSteelStrip() const {return fLateralSteelStrip;}
+  Float_t GetPassiveScintThick() const {return fPassiveScintThick;}
+  Float_t GetPhiTileSize() const {return fPhiTileSize;}
+  Float_t GetEtaTileSize() const {return fEtaTileSize;}
+  Int_t   GetNPhiSuperModule() const {return fNPhiSuperModule;}
+  Int_t   GetNPHIdiv() const {return fNPHIdiv ;}
+  Int_t   GetNETAdiv() const {return fNETAdiv ;}
+  Int_t   GetNCells()  const {return fNCells;}
+  Float_t GetSteelFrontThickness() const { return fSteelFrontThick;}
+  Float_t GetLongModuleSize() const {return fLongModuleSize;}
+
+  Float_t GetTrd1Angle() const {return fTrd1Angle;}
+  Float_t Get2Trd1Dx2()  const {return f2Trd1Dx2;}
+  Float_t GetTrd2AngleY()const {return fTrd2AngleY;}
+  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(const Int_t nSupMod, const Int_t nTower, const Int_t nIphi, const Int_t nIeta);
+  Bool_t  GetCellIndex(const Int_t absId, Int_t &nSupMod, Int_t &nTower, Int_t &nIphi, Int_t &nIeta);
+  void    GetCellPhiEtaIndexInSModule(const Int_t nTower, const Int_t nIphi, const Int_t nIeta, Int_t &iphi, Int_t &ieta);
+  Bool_t  CheckAbsCellId(Int_t ind); // replace the IsInECA
+  // ---
   Float_t AngleFromEta(Float_t eta){ // returns theta in radians for a given pseudorapidity
     return 2.0*TMath::ATan(TMath::Exp(-eta));
   }
@@ -132,11 +159,39 @@ private:
   Float_t fZLength;                    // Total length in z direction
   Float_t fGap2Active;                 // Gap between the envelop and the active material
   Int_t   fNZ;                         // Number of Towers in the Z direction
-  Int_t   fNPhi;                       // Number of Towers in the Phi Direction
+  Int_t   fNPhi;                       // Number of Towers in the PHI direction
   Float_t fSampling;                   // Sampling factor
-  
-  ClassDef(AliEMCALGeometry,8) // EMCAL geometry class 
-    
+
+  // Shish-kebab option - 23-aug-04 by PAI; COMPACT, TWIST, TRD1 and TRD2
+  Int_t   fNumberOfSuperModules;         // default is 12 = 6 * 2 
+  Float_t fSteelFrontThick;             // Thickness of the front stell face of the support box - 9-sep-04
+  Float_t fFrontSteelStrip;              // 13-may-05
+  Float_t fLateralSteelStrip;            // 13-may-05
+  Float_t fPassiveScintThick;            // 13-may-05
+  Float_t fPhiModuleSize;                // Phi -> X 
+  Float_t fEtaModuleSize;                // Eta -> Y
+  Float_t fPhiTileSize;                  // 
+  Float_t fEtaTileSize;                  // 
+  Float_t fLongModuleSize;               // 
+  Int_t   fNPhiSuperModule;              // 6 - number supermodule in phi direction
+  Int_t   fNPHIdiv;                      // number phi dvizion
+  Int_t   fNETAdiv;                      // number eta divizion
+  //
+  Int_t   fNCells;                       // number of cells in calo
+  Int_t   fNCellsInSupMod;               // number cell in super module
+  Int_t   fNCellsInTower;                // number cell in tower
+  // TRD1 options - 30-sep-04
+  Float_t fTrd1Angle;                    // angle in x-z plane (in degree) 
+  Float_t f2Trd1Dx2;                     // 2*dx2 for TRD1
+  // TRD2 options - 27-jan-07
+  Float_t fTrd2AngleY;                   // angle in y-z plane (in degree) 
+  Float_t f2Trd2Dy2;                     // 2*dy2 for TRD2
+  Float_t fEmptySpace;                   // 2mm om fred drawing
+  // Sumper module as TUBS
+  Float_t fTubsR;                        // radius of tubs 
+  Float_t fTubsTurnAngle;                // turn angle of tubs in degree
+
+  ClassDef(AliEMCALGeometry,9) // EMCAL geometry class 
     };
 
 #endif // AliEMCALGEOMETRY_H
diff --git a/EMCAL/AliEMCALGeometryOfflineTrd1.cxx b/EMCAL/AliEMCALGeometryOfflineTrd1.cxx
new file mode 100644 (file)
index 0000000..82706c4
--- /dev/null
@@ -0,0 +1,197 @@
+/**************************************************************************
+ * Copyright(c) 1998-2004, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+// GeometryOfflineTrd1 class  for EMCAL : singleton
+//*-- Author: Aleksei Pavlinov (WSU)
+
+/* $Id$*/
+
+#include <TBrowser.h>
+#include <TMath.h>
+#include <TVector2.h>
+#include <TVector3.h>
+#include <TObjArray.h>
+
+#include "AliEMCALGeometryOfflineTrd1.h"
+#include "AliEMCALShishKebabTrd1Module.h"
+#include "AliEMCALGeometry.h"
+
+ClassImp(AliEMCALGeometryOfflineTrd1);
+
+AliEMCALGeometryOfflineTrd1* AliEMCALGeometryOfflineTrd1::fgGeomOfflineTrd1=0;
+
+AliEMCALGeometryOfflineTrd1* AliEMCALGeometryOfflineTrd1::GetInstance()
+{
+  if(fgGeomOfflineTrd1==0) {
+    fgGeomOfflineTrd1 = new AliEMCALGeometryOfflineTrd1();
+  }
+  return fgGeomOfflineTrd1;
+}
+
+AliEMCALGeometryOfflineTrd1::AliEMCALGeometryOfflineTrd1() : TNamed("geomTRD1","")
+{ // this private constarctor
+  fGeometry = AliEMCALGeometry::GetInstance("SHISH_62_TRD1");
+  Init();
+}
+
+void AliEMCALGeometryOfflineTrd1::Init()
+{
+  // Super module
+  // ETA direction
+  fMaxInEta = fGeometry->GetNZ();
+  fTrd1Modules[0] =  new AliEMCALShishKebabTrd1Module();
+  fSMMaxEta = 2*fMaxInEta;
+  fSMPositionEta = new TVector2[fSMMaxEta];
+  fSMPositionEta[0] = fTrd1Modules[0]->GetCenterOfCell(1);
+  fSMPositionEta[1] = fTrd1Modules[0]->GetCenterOfCell(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);
+  }
+  // PHI direction
+  fSMPositionPhi.Set(2*fGeometry->GetNPhi());
+  fShiftOnPhi = -fGeometry->GetPhiModuleSize()*fGeometry->GetNPhi()/2;
+  for(Int_t i=0; i<fGeometry->GetNPhi(); i++) {
+    fSMPositionPhi[2*i]   = fGeometry->GetPhiModuleSize() * (double(i) + 0.25);
+    fSMPositionPhi[2*i+1] = fGeometry->GetPhiTileSize()   + fSMPositionPhi[2*i];
+  }
+  
+  // Super Module rotations
+  fNPhiSuperModule = fGeometry->GetNPhiSuperModule(); // see AliEMCALv0
+  double dphi = (fGeometry->GetArm1PhiMax() - fGeometry->GetArm1PhiMin()) / fNPhiSuperModule;
+  double phi, phiRad;
+  fSuperModuleRotationX.RotateX(TMath::Pi()); // matrix looks not so nice
+  for(Int_t i=0; i<fNPhiSuperModule; i++){
+    // rotations arround Z
+    phi    = fGeometry->GetArm1PhiMin() + dphi*(2*i+1)/2.; // phi= 70, 90, 110, 130, 150, 170
+    phiRad = phi*TMath::DegToRad();
+    if(i==1) phiRad = TMath::PiOver2();
+    fSuperModuleRotationZ[i].RotateZ(phiRad);
+    TString ntmp("rotationZ_");
+    ntmp += int(phi);
+    fNameSuperModuleRotationZ[i] = ntmp;
+    // Super Module rotation
+    fSuperModuleRotation[2*i]   = fSuperModuleRotationZ[i]; // Z
+    fSuperModuleRotation[2*i+1] = fSuperModuleRotationZ[i] * fSuperModuleRotationX; // Z*X
+  }
+  // Fill fXYZofCells
+  fXYZofCells = new TObjArray(fGeometry->GetNCells());
+  fXYZofCells->SetName("CellsInGC"); 
+  Int_t nSupMod, nTower, nIphi, nIeta, iphi, ieta;
+  for(Int_t absId=1; absId<=fGeometry->GetNCells(); absId++){
+    if(fGeometry->GetCellIndex(absId, nSupMod,nTower,nIphi,nIeta)){
+      fGeometry->GetCellPhiEtaIndexInSModule(nTower,nIphi,nIeta, iphi,ieta);
+      TVector3 *v = new TVector3;
+      v->SetX(fSMPositionEta[ieta-1].Y()); 
+      v->SetZ(fSMPositionEta[ieta-1].X()); 
+      v->SetY(fSMPositionPhi[iphi-1] + fShiftOnPhi);
+      v->Transform(fSuperModuleRotation[nSupMod-1]);
+      fXYZofCells->AddAt(v,absId-1);
+    }
+  }
+}
+
+TVector3& AliEMCALGeometryOfflineTrd1::PosInSuperModule(const Int_t nTower,const Int_t nIphi,const Int_t nIeta)
+{ // 10-nov-04
+  static Int_t iphi, ieta;
+  static TVector3 v;
+  fGeometry->GetCellPhiEtaIndexInSModule(nTower,nIphi,nIeta, iphi,ieta);
+
+  // x-radius; y-phi; eta-z;
+  v.SetXYZ(fSMPositionEta[ieta].Y(), fSMPositionPhi[iphi], fSMPositionEta[ieta].X());
+  return v;
+} 
+
+void AliEMCALGeometryOfflineTrd1::PositionInSuperModule(const Int_t iphi, const Int_t ieta, 
+double &lphi, double &leta)
+{ 
+  static Int_t ie=0;
+  lphi = fSMPositionPhi[iphi-1];
+  ie = ieta - 1;
+  if(ie<0) ie = 0;
+  if(ie>fSMMaxEta-1) ie = fSMMaxEta-1;
+  leta = fSMPositionEta[ie].X();
+}
+
+void AliEMCALGeometryOfflineTrd1::PositionInSuperModule(const Int_t nTower, const Int_t nIphi, const Int_t nIeta,
+double &lphi, double &leta)
+{
+  static Int_t iphi,ieta;
+  fGeometry->GetCellPhiEtaIndexInSModule(nTower,nIphi,nIeta,iphi,ieta);
+  PositionInSuperModule(iphi,ieta, lphi,leta);
+}
+
+TRotation* AliEMCALGeometryOfflineTrd1::Rotation(Int_t module)
+{ // module chabge from 1 to 12
+  if(module<1)  module=1;
+  if(module>12) module=12;
+  return &fSuperModuleRotation[module];
+}
+
+TVector3* AliEMCALGeometryOfflineTrd1::CellPosition(int absId)
+{ // 15-nov-04
+  if(absId<1 || absId>fXYZofCells->GetSize()) return 0;
+  return (TVector3*)fXYZofCells->At(absId-1);
+}
+
+void AliEMCALGeometryOfflineTrd1::PrintSuperModule()
+{ // 12-nov-04
+  printf(" ** Super module ** fSMMaxEta %i fSMMaxPHI %i\n ETA     eta(Z)          (X)\n",
+  fSMMaxEta,fSMPositionPhi.GetSize());
+  for(Int_t i=0; i<fSMMaxEta; i++) {
+    printf("%3i | %8.3f         %8.3f\n", i+1,fSMPositionEta[i].X(), fSMPositionEta[i].Y());
+  }
+  printf("\n PHI      (Y)\n");
+  for(Int_t i=0; i<fSMPositionPhi.GetSize(); i++) {
+    printf("%3i | %8.3f\n",i+1, fSMPositionPhi[i]);
+  }
+}
+
+void AliEMCALGeometryOfflineTrd1::PrintCell(Int_t absId)
+{
+  Int_t nSupMod, nTower, nIphi, nIeta, iphi, ieta;
+  if(fGeometry->GetCellIndex(absId, nSupMod,nTower,nIphi,nIeta)){
+     fGeometry->GetCellPhiEtaIndexInSModule(nTower,nIphi,nIeta, iphi,ieta);
+     TVector3 *v = CellPosition(absId);
+     printf("(%5i) X %8.3f Y %8.3f Z %8.3f | #sup.Mod %2i #tower %3i nIphi %1i nIeta %1i | iphi %2i ieta %2i\n",
+           absId, v->X(),v->Y(),v->Z(), nSupMod,nTower,nIphi,nIeta, iphi,ieta);
+  } else {
+    Warning("PrintCell","Wrong abs id %i\n",absId); 
+  }
+}
+
+void AliEMCALGeometryOfflineTrd1::Browse(TBrowser* b)
+{
+  if(fGeometry) b->Add(fGeometry);
+
+  for(Int_t i=0; i<fMaxInEta; i++)  if(fTrd1Modules[i]>0) b->Add(fTrd1Modules[i]);
+
+  for(Int_t i=0; i<fNPhiSuperModule; i++){ 
+    b->Add(&fSuperModuleRotationZ[i], fNameSuperModuleRotationZ[i].Data());
+    for(Int_t j=0; j<2; j++) {
+      TString name("rotationM_"); name += (2*i+j);
+      b->Add(&fSuperModuleRotation[2*i+j], name.Data());
+    }
+  }
+
+  b->Add(&fSuperModuleRotationX, "rotationX_180");
+
+  b->Add(fXYZofCells);
+}
+
+Bool_t AliEMCALGeometryOfflineTrd1::IsFolder() const
+{
+  return kTRUE;
+}
diff --git a/EMCAL/AliEMCALGeometryOfflineTrd1.h b/EMCAL/AliEMCALGeometryOfflineTrd1.h
new file mode 100644 (file)
index 0000000..5b83ba3
--- /dev/null
@@ -0,0 +1,78 @@
+#ifndef ALIEMCALGEOMETRYOFFLINETRD1_H
+#define ALIEMCALGEOMETRYOFFLINETRD1_H
+/* Copyright(c) 1998-2004, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+//_________________________________________________________________________
+// GeometryOfflineTrd1 class  for EMCAL : singleton
+//*-- Author: Aleksei Pavlinov (WSU)
+
+#include <TNamed.h>
+#include <TArrayD.h>
+#include <TRotation.h>
+
+class TBrowser;
+class TString;
+class TVector2;
+class TVector3;
+class TObjArray;
+class AliEMCALGeometry;
+class AliEMCALShishKebabTrd1Module; 
+// --- ROOT system ---
+// --- AliRoot header files ---
+
+class AliEMCALGeometryOfflineTrd1 : public TNamed {
+ public:
+  AliEMCALGeometryOfflineTrd1(const AliEMCALGeometryOfflineTrd1& geom):TNamed(geom.GetName(),geom.GetTitle()) {
+    // cpy ctor requested by Coding Convention but not yet needed
+    Fatal("Cpy ctor", "Not implemented");
+  }
+
+  virtual ~AliEMCALGeometryOfflineTrd1() { /* nothing */ };
+  static   AliEMCALGeometryOfflineTrd1* GetInstance();
+  // positon in SuperModule
+  TVector3&  PosInSuperModule(const int nTower, const int nIphi, const int nIeta); 
+
+ private:
+  AliEMCALGeometryOfflineTrd1();
+  void Init();
+
+  static  AliEMCALGeometryOfflineTrd1* fgGeomOfflineTrd1; //!
+  AliEMCALGeometry *fGeometry;                            //!  
+
+  Int_t fMaxInEta;                                        //! max module in ETA direction 
+  AliEMCALShishKebabTrd1Module* fTrd1Modules[26];         //! 
+  // positon in SuperModule
+  Int_t     fSMMaxEta;                                    // = 2*fMaxInEta(52 now)
+  TVector2* fSMPositionEta;                               //! array of TVector2 [fSMMaxEta]
+  TArrayD   fSMPositionPhi;                               //! [24] now; start from 0;
+  Double_t  fShiftOnPhi;                                  //! 
+
+  Int_t fNPhiSuperModule;                                //! number phi position for super modules  
+  TRotation fSuperModuleRotationZ[6];                    //!
+  TString   fNameSuperModuleRotationZ[6];                //!
+  TRotation fSuperModuleRotationX;                       //!
+  TRotation fSuperModuleRotation[12];                    //! 
+  // position of cells in global coordinate system
+  TObjArray *fXYZofCells;                                //! 
+ public:
+  // One Super Module
+  void PositionInSuperModule(const int iphi, const int ieta, double &lphi, double &leta);
+  void PositionInSuperModule(const int nTower, const int nIphi, const int nIeta, double &lphi, double &leta);
+  // Position towers(cells)
+  TVector3* CellPosition(int absId); // from 0 to fGeometry->GetNCells()
+  // Global System
+  TRotation* Rotation(Int_t module); // module chabge from 1 to 12;
+
+  // service methods
+  void    PrintSuperModule();       // *MENU*
+  void    PrintCell(Int_t absId=1); // *MENU*
+  virtual void Browse(TBrowser* b);
+  virtual Bool_t  IsFolder() const;
+
+  ClassDef(AliEMCALGeometryOfflineTrd1, 0) // EMCAL geometry for offline 
+};
+
+#endif // ALIEMCALGEOMETRYOFFLINETRD1_H
index 79328b6..496cc9f 100644 (file)
@@ -759,3 +759,61 @@ TList* AliEMCALJetMicroDst::MoveHistsToList(const char* name, Bool_t putToBrowse
   if(putToBrowser) gROOT->GetListOfBrowsables()->Add((TObject*)list);
   return list;
 }
+
+void AliEMCALJetMicroDst::FillH1(TList *l, Int_t ind, Double_t x, Double_t w)
+{
+  static TH1* hid=0;
+  if(l == 0) return;
+  if(ind < l->GetSize()){
+    hid = (TH1*)l->At(ind);
+    hid->Fill(x,w);
+  }
+}
+
+void AliEMCALJetMicroDst::FillH2(TList *l, Int_t ind, Double_t x, Double_t y, Double_t w)
+{
+  static TH2* hid=0;
+  if(l == 0) return;
+  if(ind < l->GetSize()){
+    hid = (TH2*)l->At(ind);
+    hid->Fill(x,y,w);
+  }
+}
+
+int AliEMCALJetMicroDst::SaveListOfHists(TList *list,const char* name,Bool_t kSingleKey,const char* opt)
+{
+  printf(" Name of out file |%s|\n", name); 
+  int save = 0;
+  if(list && list->GetSize() && strlen(name)){
+    TString nf(name); 
+    if(nf.Contains(".root") == kFALSE) nf += ".root";
+    TFile file(nf.Data(),opt);
+    TIter nextHist(list);
+    TObject* objHist=0;
+    int nh=0;
+    if(kSingleKey) {
+       file.cd();
+       list->Write(list->GetName(),TObject::kSingleKey);
+       list->ls();
+       save = 1;
+    } else {
+      while((objHist=nextHist())) { // loop over list 
+        if(objHist->InheritsFrom("TH1")) {
+          TH1* hid = (TH1*)objHist;
+          file.cd();
+          hid->Write();
+          nh++;
+          printf("Save hist. %s \n",hid ->GetName());
+        }
+      }
+      printf("%i hists. save to file -> %s\n", nh,file.GetName());
+      if(nh>0) save = 1;
+    }
+    file.Close();
+  } else {
+    printf("TAliasPAI::saveListOfHists : N O  S A V I N G \n");
+    if(list==0) printf("List of object 0 : %i \n", (int)list);
+    else printf("Size of list %i \n", list->GetSize());
+  }
+  return save;
+}
index fea92a0..8ddc6bf 100644 (file)
@@ -23,8 +23,6 @@ class TVector3;
 class TBrowser;
 
 class AliEMCALJetMicroDst: public TNamed {
-
-
   public:
   AliEMCALJetMicroDst(const char *name="jetMicroDst",
   const char *tit="jet Micro Dst for preparation of proposal");
@@ -80,7 +78,12 @@ class AliEMCALJetMicroDst: public TNamed {
   virtual Bool_t  IsFolder() const;
   virtual void Browse(TBrowser* b) const ;
 
-  static TList *MoveHistsToList(const char* name="List of Hist", Bool_t putToBrowser=kTRUE);
+  // service routine
+  static TList *MoveHistsToList(const char* name="ListOfHists", Bool_t putToBrowser=kTRUE);
+  static void FillH1(TList *l=0, Int_t ind=0, Double_t x=-99999., Double_t w=1.);
+  static void FillH2(TList *l=0, Int_t ind=0, Double_t x=-99999., Double_t w=-99999., Double_t w=1.);
+  static int  SaveListOfHists(TList *list=0, const char* name="test", Bool_t kSingleKey=kFALSE,
+  const char* opt="RECREATE");
 
   AliEMCALJetMicroDst & operator = (const AliEMCALJetMicroDst &) {
     Fatal("operator =", "not implemented") ; return *this ; }
@@ -142,6 +145,9 @@ class AliEMCALJetMicroDst: public TNamed {
 };
 
 #endif // AliEMCALJETMICRODST_H
+
+typedef AliEMCALJetMicroDst sv; // for convinience
+
 /*
 What to do
 1. Common info about event
index 98791d5..5750736 100644 (file)
@@ -38,7 +38,6 @@
 
 ClassImp(AliEMCALRecPoint)
 
-
 //____________________________________________________________________________
 AliEMCALRecPoint::AliEMCALRecPoint()
   : AliRecPoint()
@@ -135,10 +134,26 @@ Bool_t AliEMCALRecPoint::AreNeighbours(AliEMCALDigit * digit1, AliEMCALDigit * d
   AliEMCALGeometry * geom =  (AliEMCALGetter::Instance())->EMCALGeometry();
 
   Int_t relid1[2] ; 
-  geom->AbsToRelNumbering(digit1->GetId(), relid1) ; 
+    //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(nTower,nIphi,nIeta, iphi,ieta);
+       relid1[0]=ieta;
+       relid1[1]=iphi;
+//  geom->AbsToRelNumbering(digit1->GetId(), relid1) ; 
 
   Int_t relid2[2] ; 
-  geom->AbsToRelNumbering(digit2->GetId(), relid2) ; 
+    //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(nTower1,nIphi1,nIeta1, iphi1,ieta1);
+       relid2[0]=ieta1;
+       relid2[1]=iphi1;
+//  geom->AbsToRelNumbering(digit2->GetId(), relid2) ; 
   
   Int_t rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ;  
   Int_t coldiff = TMath::Abs( relid1[1] - relid2[1] ) ;  
@@ -293,13 +308,20 @@ void AliEMCALRecPoint::EvalAll(Float_t logWeight,TClonesArray * digits)
   // Evaluates all shower parameters
 
   EvalLocalPosition(logWeight, digits) ;
+  printf("eval position done\n");
   EvalElipsAxis(logWeight, digits) ;
+  printf("eval axis done\n");
   EvalDispersion(logWeight, digits) ;
-  EvalCoreEnergy(logWeight, digits);
+  printf("eval dispersion done\n");
+ // EvalCoreEnergy(logWeight, digits);
+ // printf("eval energy done\n");
   EvalTime(digits) ;
+  printf("eval time done\n");
 
   EvalPrimaries(digits) ;
+  printf("eval pri done\n");
   EvalParents(digits);
+  printf("eval parent done\n");
 }
 
 //____________________________________________________________________________
@@ -320,7 +342,7 @@ void  AliEMCALRecPoint::EvalDispersion(Float_t logWeight, TClonesArray * digits)
   if (!fLocPos.X() || !fLocPos.Y()) 
     EvalLocalPosition(logWeight, digits) ;
   
-  const Float_t kDeg2Rad = TMath::DegToRad() ; 
+  //Sub const Float_t kDeg2Rad = TMath::DegToRad() ; 
     
   Float_t cluEta =  fLocPos.X() ; 
   Float_t cluPhi =  fLocPos.Y() ; 
@@ -335,8 +357,24 @@ void  AliEMCALRecPoint::EvalDispersion(Float_t logWeight, TClonesArray * digits)
     digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit])  ;
     Float_t etai = 0.;
     Float_t phii = 0.;
-    geom->EtaPhiFromIndex(digit->GetId(), etai, 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(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) ;
 
@@ -351,7 +389,7 @@ void  AliEMCALRecPoint::EvalDispersion(Float_t logWeight, TClonesArray * digits)
     d = 0. ; 
 
   fDispersion = TMath::Sqrt(d) ;
+      printf("Dispersion: = %f", fDispersion) ;
 }
 
 //____________________________________________________________________________
@@ -367,15 +405,23 @@ void AliEMCALRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digit
   Int_t iDigit;
   Float_t cluEta = 0;
   Float_t cluPhi = 0;
-  const Float_t kDeg2Rad = TMath::DegToRad();
+ //Sub const Float_t kDeg2Rad = TMath::DegToRad();
 
   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
     digit = dynamic_cast<AliEMCALDigit *>(digits->At(fDigitsList[iDigit])) ;
 
     Float_t etai ;
     Float_t phii ;
-    geom->EtaPhiFromIndex(digit->GetId(), etai, phii);
-    phii = phii * kDeg2Rad;
+    //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(nTower,nIphi,nIeta, iphi,ieta);
+       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 );
@@ -394,7 +440,7 @@ void AliEMCALRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digit
   fLocPos.SetY(cluPhi);
   fLocPos.SetZ(geom->GetIP2ECASection());
     
-  if (gDebug==2)
+//  if (gDebug==2)
     printf("EvalLocalPosition: eta,phi,r = %f,%f,%f", fLocPos.X(), fLocPos.Y(), fLocPos.Z()) ; 
   fLocPosM = 0 ;
 }
@@ -445,6 +491,7 @@ void  AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
   AliEMCALDigit * digit ;
 
   AliEMCALGeometry * geom = (AliEMCALGetter::Instance())->EMCALGeometry();
+  TString gn(geom->GetName());
 
   Int_t iDigit;
   
@@ -452,15 +499,29 @@ void  AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
     digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit])  ;
     Float_t etai = 0. ;
     Float_t phii = 0. ; 
-    geom->EtaPhiFromIndex(digit->GetId(), etai, phii);
-    phii = phii * kDeg2Rad;
+    if(gn.Contains("SHISH")) {
+    //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(nTower,nIphi,nIeta, iphi,ieta);
+      etai=(Float_t)ieta;
+      phii=(Float_t)iphi;
+    } else {
+      geom->EtaPhiFromIndex(digit->GetId(), etai, phii);
+      phii = phii * kDeg2Rad;
+    }
+
     Double_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
 
     dxx  += w * etai * etai ;
     x    += w * etai ;
     dzz  += w * phii * phii ;
     z    += w * phii ; 
+
     dxz  += w * etai * phii ; 
+
     wtot += w ;
   }
 
@@ -491,6 +552,7 @@ void  AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
     fLambda[1]= 0. ;
   }
 
+  //  printf("Evalaxis: lambdas  = %f,%f", fLambda[0],fLambda[1]) ; 
 
 }
 
index b7421be..d3823cc 100644 (file)
@@ -238,5 +238,4 @@ void EMCALHits2Digits (Bool_t split=kFALSE, TString fileName = "galice.root") {
   }
 }
 
-
  
index 204a0dc..d5caade 100644 (file)
@@ -31,7 +31,6 @@
 #include "AliEMCALClusterizerv1.h"
 #include "AliEMCALPIDv1.h"
 #include "AliEMCALGetter.h"
-#include "AliEMCALTracker.h"
 #include "AliRawReaderFile.h"
 
 ClassImp(AliEMCALReconstructor)
@@ -100,6 +99,19 @@ void AliEMCALReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
   
   Int_t eventNumber = runLoader->GetEventNumber() ;
 
+  TString headerFile(runLoader->GetFileName()) ; 
+  TString branchName(runLoader->GetEventFolder()->GetName()) ;  
+
+  AliEMCALPIDv1 pid(headerFile, branchName);
+
+
+  // do current event; the loop over events is done by AliReconstruction::Run()
+  pid.SetEventRange(eventNumber, eventNumber) ; 
+  if ( Debug() ) 
+    pid.ExecuteTask("deb all") ;
+  else 
+    pid.ExecuteTask("") ;
   // Creates AliESDtrack from AliEMCALRecParticles 
   AliEMCALGetter::Instance()->Event(eventNumber, "P") ; 
   TClonesArray *recParticles = AliEMCALGetter::Instance()->RecParticles();
@@ -123,10 +135,3 @@ void AliEMCALReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
     delete et;
   }
 }
-
-AliTracker* AliEMCALReconstructor::CreateTracker(AliRunLoader* runLoader) const
-{
-// creates the EMCAL tracker
-  if (!runLoader) return NULL; 
-  return new AliEMCALTracker(runLoader);
-}
index dbe71eb..fad0518 100644 (file)
@@ -40,8 +40,7 @@ public:
   virtual ~AliEMCALReconstructor() ; //dtor
 
   Bool_t                     Debug() const { return fDebug ; }
-  AliTracker *CreateTracker(AliRunLoader* runLoader) const;  
-  virtual void FillESD(AliRunLoader* runLoader, AliESD* esd) const ;
+  virtual void               FillESD(AliRunLoader* runLoader, AliESD* esd) const ;
   virtual void               Reconstruct(AliRunLoader* runLoader) const ;
   virtual void               Reconstruct(AliRunLoader* runLoader, AliRawReaderFile * rawreader) const ;
   
index 89c786d..cf18702 100644 (file)
@@ -50,7 +50,9 @@
 
 // --- ROOT system ---
 #include "TBenchmark.h"
-              //#include "TObjectTable.h"
+#include "TH1.h"
+#include "TBrowser.h"
+//#include "TObjectTable.h"
 
 // --- Standard library ---
 #include "stdlib.h"
@@ -62,7 +64,8 @@
 #include "AliEMCALHit.h"
 #include "AliEMCALSDigitizer.h"
 #include "AliEMCALGeometry.h"
-              //#include "AliMemoryWatcher.h"
+#include "AliEMCALJetMicroDst.h"
+//#include "AliMemoryWatcher.h"
 
 ClassImp(AliEMCALSDigitizer)
            
@@ -70,8 +73,10 @@ ClassImp(AliEMCALSDigitizer)
   AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("","") 
 {
   // ctor
-  fFirstEvent = fLastEvent  = 0 ;  
+  Info("AliEMCALSDigitizer()", " CTOR # 1 #");  
+  fFirstEvent = fLastEvent  = fControlHists = 0;  
   fDefaultInit = kTRUE ; 
+  fHists = 0;
 }
 
 //____________________________________________________________________________ 
@@ -81,16 +86,19 @@ AliEMCALSDigitizer::AliEMCALSDigitizer(const char * alirunFileName,
   fEventFolderName(eventFolderName)
 {
   // ctor
-  fFirstEvent = fLastEvent  = 0 ; // runs one event by defaut  
+  Info("AliEMCALSDigitizer()", " CTOR # 2 #");  
+  fFirstEvent = fLastEvent  = fControlHists = 0 ; // runs one event by defaut  
   Init();
   InitParameters() ; 
   fDefaultInit = kFALSE ; 
+  fHists = 0;
 }
 
 
 //____________________________________________________________________________ 
 AliEMCALSDigitizer::AliEMCALSDigitizer(const AliEMCALSDigitizer & sd) : TTask(sd) {
   //cpy ctor 
+  Info("AliEMCALSDigitizer()", " CPY CTOR # 3 #");  
 
   fFirstEvent    = sd.fFirstEvent ; 
   fLastEvent     = sd.fLastEvent ;
@@ -144,26 +152,25 @@ void AliEMCALSDigitizer::Init(){
 //____________________________________________________________________________ 
 void AliEMCALSDigitizer::InitParameters()
 {
-  // Initializes parameters
-  fA                      = 0;
-  fB                      = 10000000.;
-
   AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
   const AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
   if (geom->GetSampling() == 0.) {
     Fatal("InitParameters", "Sampling factor not set !") ; 
   }
-//   else
-//     Info("InitParameters", "Sampling factor set to %f", geom->GetSampling()) ; 
-  
-  // this threshold corresponds approximately to 100 MeV
-  fECPrimThreshold     = 100E-3;
+  // Initializes parameters
+  fA         = 0;
+  fB         = 1.e+7;
+  fSampling  = geom->GetSampling();
+
+ // threshold for deposit energy of hit
+  fECPrimThreshold  = 0.; // 24-nov-04 - was 1.e-6;
 }
 
 //____________________________________________________________________________
 void AliEMCALSDigitizer::Exec(Option_t *option) 
 { 
   // Collects all hit of the same tower into digits
+  TString o(option); o.ToUpper();
   if (strstr(option, "print") ) {
     Print() ; 
     return ; 
@@ -183,16 +190,19 @@ void AliEMCALSDigitizer::Exec(Option_t *option)
 
   if (fLastEvent == -1) 
     fLastEvent = gime->MaxEvent() - 1 ;
-  else 
-    fLastEvent = TMath::Min(fFirstEvent, gime->MaxEvent()); // only ine event at the time
+  else {
+    if(o != "EXACT") fLastEvent = TMath::Min(fFirstEvent, gime->MaxEvent());
+    fLastEvent = TMath::Min(fLastEvent, gime->MaxEvent());
+    printf("AliEMCALSDigitizer::Exec : option: %s | %i -> %i events : Max events %i \n \n", 
+    option, fFirstEvent, fLastEvent, gime->MaxEvent());
+  }
   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
 
-  Int_t ievent ;   
-
-  //AliMemoryWatcher memwatcher;
-
-  for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
-    gime->Event(ievent,"H") ;  
+  Int_t ievent;
+  Float_t energy=0.; // de * fSampling - 23-nov-04
+  for (ievent = fFirstEvent; ievent < fLastEvent; ievent++) {
+    if(ievent%100==0 || ievent==fLastEvent-1) printf(" processed event %i \n", ievent);  
+    gime->Event(ievent,"XH"); // read primaries and hits onl  
     TTree * treeS = gime->TreeS(); 
     TClonesArray * hits = gime->Hits() ; 
     TClonesArray * sdigits = gime->SDigits() ;
@@ -207,28 +217,30 @@ void AliEMCALSDigitizer::Exec(Option_t *option)
       //=========== Get the EMCAL branch from Hits Tree for the Primary iprim
       gime->Track(iprim) ;
       Int_t i;
+      AliEMCALGeometry *geom = gime->EMCALGeometry(); 
       for ( i = 0 ; i < hits->GetEntries() ; i++ ) {
        AliEMCALHit * hit = dynamic_cast<AliEMCALHit*>(hits->At(i)) ;
        AliEMCALDigit * curSDigit = 0 ;
        AliEMCALDigit * sdigit = 0 ;
        Bool_t newsdigit = kTRUE; 
 
+        // hit->GetId() - Absolute Id number EMCAL segment
+       if(geom->CheckAbsCellId(hit->GetId())) { // was IsInECA(hit->GetId())
+          energy = hit->GetEnergy() * fSampling; // 23-nov-04
+         if(energy >  fECPrimThreshold )
        // Assign primary number only if deposited energy is significant
-       AliEMCALGeometry * geom = gime->EMCALGeometry() ; 
-       if( geom->IsInECA(hit->GetId()) ){
-         if( hit->GetEnergy() >  fECPrimThreshold )
            curSDigit =  new AliEMCALDigit( hit->GetPrimary(),
                                            hit->GetIparent(), hit->GetId(), 
-                                           Digitize(hit->GetEnergy()), hit->GetTime() ) ;
+                                           Digitize(energy), hit->GetTime() ) ;
          else
            curSDigit =  new AliEMCALDigit( -1               , 
                                            -1               ,
                                            hit->GetId(), 
-                                           Digitize(hit->GetEnergy()), hit->GetTime() ) ;
-       }
-       else{
-         newsdigit = kFALSE;
-       }
+                                           Digitize(energy), hit->GetTime() ) ;
+       } else {
+          Warning("Exec"," abs id %i is bad \n", hit->GetId());
+          newsdigit = kFALSE;
+        }
 
        Int_t check = 0 ;
 
@@ -256,12 +268,21 @@ void AliEMCALSDigitizer::Exec(Option_t *option)
 
     Int_t nPrimarymax = -1 ; 
     Int_t i ;
+    Double_t e=0.,esum=0.;
+    sv::FillH1(fHists, 0, double(sdigits->GetEntriesFast()));
     for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { 
       AliEMCALDigit * sdigit = dynamic_cast<AliEMCALDigit *>(sdigits->At(i)) ;
       sdigit->SetIndexInList(i) ;
+
+      sv::FillH1(fHists, 2, double(sdigit->GetAmp()));
+      e = double(Calibrate(sdigit->GetAmp()));
+      esum += e;
+      sv::FillH1(fHists, 3, e);
+      sv::FillH1(fHists, 4, double(sdigit->GetId()));
     }
+    if(esum>0.) sv::FillH1(fHists, 1, esum);
     
-    for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) {   
+    for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { // for what 
       if (((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) > nPrimarymax)
        nPrimarymax = ((dynamic_cast<AliEMCALDigit *>(sdigits->At(i)))->GetNprimary()) ;
     }
@@ -277,15 +298,12 @@ void AliEMCALSDigitizer::Exec(Option_t *option)
     gime->WriteSDigits("OVERWRITE");
     
     //NEXT - SDigitizer
-
-    gime->WriteSDigitizer("OVERWRITE");
+    if(ievent == fFirstEvent) gime->WriteSDigitizer("OVERWRITE");  // why in event cycle ?
     
     if(strstr(option,"deb"))
       PrintSDigits(option) ;  
-    
-    //gObjectTable->Print() ; 
-    //memwatcher.Watch(ievent); 
-  }// event loop
+  }
+  //  gime->WriteSDigitizer("OVERWRITE"); // 12-jan-04
 
   Unload();
   
@@ -293,7 +311,7 @@ void AliEMCALSDigitizer::Exec(Option_t *option)
   
   if(strstr(option,"tim")){
     gBenchmark->Stop("EMCALSDigitizer"); 
-    printf("Exec: took %f seconds for SDigitizing %f seconds per event", 
+    printf("\n Exec: took %f seconds for SDigitizing %f seconds per event\n", 
         gBenchmark->GetCpuTime("EMCALSDigitizer"), gBenchmark->GetCpuTime("EMCALSDigitizer")/nEvents ) ; 
   }
 
@@ -303,17 +321,26 @@ void AliEMCALSDigitizer::Exec(Option_t *option)
 }
 
 
+void AliEMCALSDigitizer::Print1(Option_t * option)
+{
+  Print(); 
+  PrintSDigits(option);
+}
+
 //__________________________________________________________________
-void AliEMCALSDigitizer::Print()const
+void AliEMCALSDigitizer::Print() const
 { 
   // Prints parameters of SDigitizer
-  printf("Print: \n------------------- %s -------------", GetName() ) ; 
+  printf("Print: \n------------------- %s -------------\n", GetName() ) ; 
+  printf("   fInit                                 %i\n", int(fInit));
+  printf("   fFirstEvent                           %i\n", fFirstEvent);
+  printf("   fLastEvent                            %i\n", fLastEvent);
   printf("   Writing SDigits to branch with title  %s\n", fEventFolderName.Data()) ;
-  printf("   with digitization parameters  A = %f\n", fA) ; 
-  printf("                                 B = %f\n", fB) ;
-  printf("   Threshold for EC Primary assignment= %f\n", fECPrimThreshold)  ; 
+  printf("   with digitization parameters       A = %f\n", fA) ; 
+  printf("                                      B = %f\n", fB) ;
+  printf("   Threshold for EC Primary assignment  = %f\n", fECPrimThreshold)  ;
+  printf("   Sampling                             = %f\n", fSampling);
   printf("---------------------------------------------------\n") ;
-
 }
 
 //__________________________________________________________________
@@ -339,19 +366,20 @@ void AliEMCALSDigitizer::PrintSDigits(Option_t * option)
   
   printf("\n") ;  
   printf("event %i", gAlice->GetEvNumber()) ;
-  printf("\n      Number of entries in SDigits list %i", sdigits->GetEntriesFast()); 
+  printf(" Number of entries in SDigits list %i", sdigits->GetEntriesFast()); 
   if(strstr(option,"all")||strstr(option,"EMC")){
     
     //loop over digits
     AliEMCALDigit * digit;
     printf("\n   Id  Amplitude    Time          Index Nprim: Primaries list \n") ;    
-    Int_t index ;
+    Int_t index, isum=0;
     char * tempo = new char[8192]; 
     for (index = 0 ; index < sdigits->GetEntries() ; index++) {
       digit = dynamic_cast<AliEMCALDigit *>( sdigits->At(index) ) ;
       sprintf(tempo, "\n%6d  %8d    %6.5e %4d      %2d :",
              digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;  
-      printf(tempo); 
+      printf(tempo);
+      isum += digit->GetAmp();
       
       Int_t iprimary;
       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
@@ -360,7 +388,8 @@ void AliEMCALSDigitizer::PrintSDigits(Option_t * option)
       }         
     }
     delete tempo ;
-  }
+    printf("\n** Sum %i : %10.3f GeV/c **\n ", isum, double(isum)*1.e-6);
+  } else printf("\n");
 }
 
 //____________________________________________________________________________ 
@@ -372,3 +401,31 @@ void AliEMCALSDigitizer::Unload() const
   loader->UnloadHits() ; 
   loader->UnloadSDigits() ; 
 }
+
+void AliEMCALSDigitizer::Browse(TBrowser* b)
+{
+  if(fHists) b->Add(fHists);
+  TTask::Browse(b);
+}
+
+TList *AliEMCALSDigitizer::BookControlHists(const int var)
+{ // 22-nov-04
+  gROOT->cd();
+  AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName);
+  const AliEMCALGeometry *geom = gime->EMCALGeometry() ;
+  if(var>=1){
+    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);
+    new TH1F("HSDigiEnergy","EMCAL cell energy", 1000, 0.0, 100.);
+    new TH1F("HSDigiAbsId","EMCAL absID for sdigits",
+    geom->GetNCells(), 0.5, Double_t(geom->GetNCells())+0.5);
+  }
+  fHists = sv::MoveHistsToList("EmcalSDigiControlHists", kFALSE);
+  return fHists;
+}
+
+void AliEMCALSDigitizer::SaveHists(const char* name, Bool_t kSingleKey, const char* opt)
+{
+  sv::SaveListOfHists(fHists, name, kSingleKey, opt); 
+}
index 4429ce5..ec5065b 100644 (file)
@@ -17,6 +17,9 @@
 // --- ROOT system ---
 #include "TTask.h"
 class TFile ;
+class TList;
+class TBrowser;
+//class TBrowser;
 
 // --- Standard library ---
 
@@ -32,16 +35,25 @@ public:
   virtual ~AliEMCALSDigitizer(); // dtor
 
   Float_t       Calibrate(Int_t amp)const {return (amp - fA)/fB ; }
-  Int_t         Digitize(Float_t Energy)const { return (Int_t ) ( fA + Energy*fB); }
+  Int_t         Digitize(Float_t energy)const { return (Int_t ) (fA + energy*fB); }
   virtual void  Exec(Option_t *option); 
   Int_t         GetSDigitsInRun() const {return fSDigitsInRun ;}  
-  virtual void  Print() const ;
+  virtual void  Print() const;
+  void          Print1(Option_t *option="all");  // *MENU*
   void          SetEventFolderName(TString name) { fEventFolderName = name ; }
   void          SetEventRange(Int_t first=0, Int_t last=-1) {fFirstEvent=first; fLastEvent=last; }
 
   Bool_t operator == (const AliEMCALSDigitizer & sd) const ;
   const AliEMCALSDigitizer & operator = (const AliEMCALSDigitizer & /*sd*/) {return *this ;}
 
+  virtual void Browse(TBrowser* b);
+  // hists
+  void   SetControlHists(const Int_t var=0) {fControlHists=var;}
+  Int_t  GetControlHist() const {return fControlHists;}
+  TList *GetListOfHists() {return fHists;}
+  TList* BookControlHists(const int var=0);
+  void   SaveHists(const char* name="RF/TRD1/Digitizations/SDigiVar?", 
+  Bool_t kSingleKey=kTRUE, const char* opt="RECREATE"); // *MENU*
 
 private:
   void     Init() ;
@@ -55,13 +67,16 @@ private:
   Float_t fECPrimThreshold ;       // To store primary if EC Shower Elos > threshold
   Bool_t  fDefaultInit;            //! Says if the task was created by defaut ctor (only parameters are initialized)
   TString fEventFolderName;        // event folder name
-  Bool_t  fInit ;                  //! tells if initialisation wennt OK, will revent exec if not
+  Bool_t  fInit ;                  //! tells if initialisation went OK, will revent exec if not
   Int_t   fSDigitsInRun ;          //! Total number of sdigits in one run
   Int_t   fFirstEvent;             // first event to process
   Int_t   fLastEvent;              // last  event to process
+  Float_t fSampling;               // See AliEMCALGeometry
+  // Control hists
+  Int_t   fControlHists;          //!
+  TList  *fHists;                 //!
 
   ClassDef(AliEMCALSDigitizer,5)  // description 
-
 };
 
 #endif // AliEMCALSDigitizer_H
diff --git a/EMCAL/AliEMCALShishKebabModule.cxx b/EMCAL/AliEMCALShishKebabModule.cxx
new file mode 100644 (file)
index 0000000..94e89cd
--- /dev/null
@@ -0,0 +1,183 @@
+/**************************************************************************
+ * Copyright(c) 1998-2004, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+//*-- Author: Aleksei Pavlinov(WSU)
+#include "AliEMCALShishKebabModule.h"
+#include "AliEMCALGeometry.h"
+#include <TMath.h>
+#include <TGraph.h>
+
+#include <assert.h>
+
+ClassImp(AliEMCALShishKebabModule)
+
+  AliEMCALGeometry *AliEMCALShishKebabModule::fgGeometry=0; 
+  Double_t AliEMCALShishKebabModule::fga=0.; 
+  Double_t AliEMCALShishKebabModule::fgb=0.; 
+  Double_t AliEMCALShishKebabModule::fgr=0.; 
+
+AliEMCALShishKebabModule::AliEMCALShishKebabModule(const double theta) : TNamed()
+{ // theta in radians ; first object shold be with theta=pi/2.
+  fTheta = theta;
+  if(fgGeometry==0) {
+    if(GetParameters()) {
+      DefineFirstModule();
+    }
+  } else Warning("AliEMCALShishKebabModule(theta)","You should call this constractor just once !!");
+  DefineName(fTheta);
+}
+
+AliEMCALShishKebabModule::AliEMCALShishKebabModule(AliEMCALShishKebabModule &leftNeighbor) : TNamed()
+{ // 22-sep-04
+  TObject::SetUniqueID(leftNeighbor.GetUniqueID()+1);
+  Init(leftNeighbor.GetA(),leftNeighbor.GetB());
+}
+
+void AliEMCALShishKebabModule::Init(const double A,const double B)
+{ 
+  Double_t thetaMin, thetaMax, par[4];
+  Int_t npar=0;
+  if(A<0){
+    //    DefineSecondModuleFirstAssumption();
+    thetaMax = TMath::ATan2(fgr, 1.4*fga);
+    thetaMin = TMath::ATan2(fgr, 1.6*fga);
+    fTheta   = Solve(AliEMCALShishKebabModule::Y2, thetaMin,thetaMax,npar,par);
+  } else{
+    npar = 4;
+    par[0] = fga;
+    par[1] = fgr;
+    par[2] = A;
+    par[3] = B;
+    Double_t x = fgr/A;
+    thetaMax = TMath::ATan2(fgr,x + 0.5*fga);
+    thetaMin = TMath::ATan2(fgr,x + 1.5*fga);
+    fTheta   = Solve(AliEMCALShishKebabModule::YALL, thetaMin,thetaMax,npar,par);
+  }
+
+  Double_t rOK = fgr / TMath::Sin(fTheta) + (fga/2.)/TMath::Tan(fTheta) + fgb/2.;
+  fOK.SetMagPhi(rOK, fTheta); 
+  // have to define A and B
+  fA = TMath::Tan(fTheta);
+  fB = -fga/(2.*TMath::Cos(fTheta));
+  DefineName(fTheta);
+}
+
+void AliEMCALShishKebabModule::DefineFirstModule()
+{
+  fOK.Set(fga/2., fgr + fgb/2.); // position the center of module vs o
+
+  fB = fga/2.;    // z=fB
+  fA = -999.;     // fA=infinite in this case
+  TObject::SetUniqueID(1); //
+}
+
+void AliEMCALShishKebabModule::DefineSecondModuleFirstAssumption()
+{ // Keep for testing and checking
+  // cos(theta) << 1, theta ~ pi/2.; a/r = 11.4/462.54 = 0.0246465 << 1; 
+  // theta=1.53382  from this case; theta=1.533869 from TGraph::Zero 
+  Double_t x = (3.*fga)/(2.*fgr);
+  fTheta = TMath::ACos(x);
+  /*
+  Double_t rOK = fgr / TMath::Sin(fTheta) + (fga/2.)/TMath::Tan(fTheta) + fgb/2.;
+  fOK.SetMagPhi(rOK, fTheta); 
+  // have to define A and B
+  fA = TMath::Tan(fTheta);
+  fB = -fga/(2.*TMath::Cos(fTheta));
+  DefineName(fTheta);
+  */
+}
+
+Double_t AliEMCALShishKebabModule::Solve(Double_t (*fcn)(Double_t*,Double_t*), 
+Double_t xmin, Double_t xmax, Int_t npar, Double_t *par, Double_t eps, Int_t maxIter)
+{
+  if(npar); // unused now
+  TGraph gr;
+  double X,Y;
+  Int_t k = 0;
+  gr.Zero(k, xmin,xmax, eps, X,Y, maxIter); // remember initial interval
+  while(k!=2) {
+    Y = fcn(&X, par); 
+    gr.Zero(k, xmin,xmax, eps, X,Y, maxIter);
+  }
+  return X;
+}
+
+Double_t AliEMCALShishKebabModule::Y2(double *x, double *par)
+{ // For position calulation of second module
+  if(par);
+  Double_t theta = x[0];
+  Double_t cos = TMath::Cos(theta);
+  Double_t sin = TMath::Sin(theta);
+  Double_t y1  = fgr*cos/sin + fga/(2.*sin) - fga*sin;       
+  Double_t y2  = fga, y = y1-y2;;
+  //  printf(" theta %f Y %12.5e \n", theta, y);
+  return y;
+}
+
+Double_t AliEMCALShishKebabModule::YALL(double *x, double *par)
+{ // For position calulation of 3th, 4th to 30th modules
+  Double_t a=par[0], r=par[1], A=par[2], B=par[3]; 
+  Double_t theta = x[0];
+  Double_t cos = TMath::Cos(theta);
+  Double_t sin = TMath::Sin(theta);
+
+  Double_t y1  = r + a*cos;       
+  Double_t y2  = A*(r*cos/sin + a/(2.*sin) - a*sin) + B;
+  Double_t y   = y1-y2;
+  //  printf(" theta %f Y %12.5e \n", theta, y);
+  return y;
+}
+
+void AliEMCALShishKebabModule::DefineName(const double theta)
+{
+  char name[100];
+  // sprintf(name,"theta_%5.2f",theta*180./TMath::Pi());
+  sprintf(name,"%2i(%5.2f)", TObject::GetUniqueID(), theta*180./TMath::Pi());
+  SetName(name);
+}
+
+Bool_t AliEMCALShishKebabModule::GetParameters()
+{
+  fgGeometry = AliEMCALGeometry::GetInstance();
+  //  if(!fgGeometry) assert(0);
+  if(!fgGeometry) {
+    Warning("GetParameters()"," No geometry ");
+    return kFALSE; 
+  }
+
+  fga = (Double_t)fgGeometry->GetPhiModuleSize();
+  fgb = (Double_t)fgGeometry->GetLongModuleSize();
+  fgr = (Double_t)(fgGeometry->GetIPDistance() + fgGeometry->GetSteelFrontThickness());
+  Print(0);
+  return kTRUE;
+}
+
+// service methods
+void AliEMCALShishKebabModule::Print(const int pri) const
+{
+  if(pri>=0) {
+    Info("Print()", " a %7.2f | b %7.2f | r %7.2f ", fga, fgb, fgr);
+    printf(" fTheta %f : %5.2f : cos(theta) %f\n", fTheta, GetThetaInDegree(),TMath::Cos(fTheta)); 
+    if(pri>0) {
+      printf("%i %s | theta %f -> %f\n", GetUniqueID(), GetName(), fTheta, fOK.Phi());
+      printf(" A %f B %f \n", fA, fB);
+
+      fOK.Dump();
+    }
+  }
+}
+
diff --git a/EMCAL/AliEMCALShishKebabModule.h b/EMCAL/AliEMCALShishKebabModule.h
new file mode 100644 (file)
index 0000000..59a9e36
--- /dev/null
@@ -0,0 +1,63 @@
+#ifndef ALIEMCALSHISHKEBABMODULE_H
+#define ALIEMCALSHISHKEBABMODULE_H
+
+/* Copyright(c) 1998-2004, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+#include "TNamed.h"
+#include "TMath.h"
+#include "TVector2.h"
+
+class AliEMCALGeometry;
+
+class AliEMCALShishKebabModule : public TNamed {
+ public:
+  AliEMCALShishKebabModule(const double theta=TMath::Pi()/2.);
+  AliEMCALShishKebabModule(AliEMCALShishKebabModule &leftNeighbor);
+  void Init(const double A,const double B);
+
+  virtual ~AliEMCALShishKebabModule(void) {}
+  Bool_t GetParameters();
+  void DefineName(const double theta);
+  void DefineFirstModule();
+  void DefineSecondModuleFirstAssumption(); // need for testing
+
+  Double_t Solve(Double_t (*fcn)(Double_t*, Double_t*), Double_t xmin=0., Double_t xmax=1.,
+  Int_t npar=0, Double_t *par=0, Double_t eps=1.0e-8, Int_t maxIter=1000);
+
+  static Double_t Y2(double *x, double *par);
+  static Double_t YALL(double *x, double *par);
+
+  Double_t GetTheta() const {return fTheta;}
+  Double_t GetThetaInDegree() const {return fTheta*180./TMath::Pi();}
+
+  Double_t GetPosX() {return fOK.Y();}
+  Double_t GetPosZ() {return fOK.X();}
+  Double_t GetPosXfromR() {return fOK.Y() - fgr;}
+  Double_t GetA() {return fA;}
+  Double_t GetB() {return fB;}
+
+  // geometry info
+  static AliEMCALGeometry *fgGeometry; //!
+  static Double_t fga; // default 11.2cm 
+  static Double_t fgb; // 
+  // radius to IP
+  static Double_t fgr;
+
+  TVector2 fOK; // position the module center x->y; z->x;
+  Double_t fA;  // parameters of line = y = A*z + B
+  Double_t fB;  // 
+  // service methods
+  void Print(const int pri=1) const;  // *MENU*
+ protected:
+  // size of SK module
+  Double_t fTheta; // theta for SK module
+
+  ClassDef(AliEMCALShishKebabModule,0) // Turned Shish-Kebab module 
+};
+
+#endif
+/* To do
+ 1. Insert position the center of towers - 2 additional TVector2
+ */
diff --git a/EMCAL/AliEMCALShishKebabTrd1Module.cxx b/EMCAL/AliEMCALShishKebabTrd1Module.cxx
new file mode 100644 (file)
index 0000000..f0d07d3
--- /dev/null
@@ -0,0 +1,155 @@
+/**************************************************************************
+ * Copyright(c) 1998-2004, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+//*-- Author: Aleksei Pavlinov(WSU)
+
+#include "AliEMCALShishKebabTrd1Module.h"
+#include <assert.h>
+#include <TMath.h>
+#include "AliEMCALGeometry.h"
+
+ClassImp(AliEMCALShishKebabTrd1Module)
+
+  AliEMCALGeometry *AliEMCALShishKebabTrd1Module::fgGeometry=0; 
+  Double_t AliEMCALShishKebabTrd1Module::fga=0.; 
+  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::fgtanBetta=0; //
+
+AliEMCALShishKebabTrd1Module::AliEMCALShishKebabTrd1Module(const double theta) : TNamed()
+{ // theta in radians ; first object shold be with theta=pi/2.
+  fTheta = theta;
+  if(fgGeometry==0) {
+    if(GetParameters()) {
+      DefineFirstModule();
+    }
+  } else Warning("AliEMCALShishKebabTrd1Module(theta)","You should call this constractor just once !!");
+  DefineName(fTheta);
+}
+
+AliEMCALShishKebabTrd1Module::AliEMCALShishKebabTrd1Module(AliEMCALShishKebabTrd1Module &leftNeighbor) : TNamed()
+{ // 22-sep-04
+  //  printf("** Left Neighbor : %s **\n", leftNeighbor.GetName());
+  TObject::SetUniqueID(leftNeighbor.GetUniqueID()+1);
+  fTheta  = leftNeighbor.GetTheta() - fgangle; 
+  Init(leftNeighbor.GetA(),leftNeighbor.GetB());
+}
+
+void AliEMCALShishKebabTrd1Module::Init(const double A,const double B)
+{ // Define parameter module from parameters A,B from previos.
+  Double_t yl = (fgb/2)*TMath::Sin(fTheta) + (fga/2)*TMath::Cos(fTheta) + fgr, y = yl;
+  Double_t xl = (yl - B) / A;     // y=A*x+B
+
+  //  Double_t xp1 = (fga/2. + fgb/2.*fgtanBetta)/(TMath::Sin(fTheta) + fgtanBetta*TMath::Cos(fTheta));
+  //  printf(" xp1 %9.3f \n ", xp1);
+  // xp1 == xp => both methods give the same results - 3-feb-05
+  Double_t alpha = TMath::Pi()/2. + fgangle/2;
+  Double_t xt = (fga+fga2)*TMath::Tan(fTheta)*TMath::Tan(alpha)/(4.*(1.-TMath::Tan(fTheta)*TMath::Tan(alpha)));
+  Double_t yt = xt / TMath::Tan(fTheta), xp = TMath::Sqrt(xt*xt + yt*yt);
+  Double_t x  = xl + xp;
+  fOK.Set(x, y);
+  //  printf(" yl %9.3f | xl %9.3f | xp %9.3f \n", yl, xl, xp);
+
+  // have to define A and B; 
+  Double_t yCprev = fgr + fga*TMath::Cos(fTheta);
+  Double_t xCprev = (yCprev - B) / A;
+  Double_t xA     = xCprev + fga*TMath::Sin(fTheta), yA = fgr;
+
+  fThetaA = fTheta - fgangle/2.;
+  fA = TMath::Tan(fThetaA); // !!
+  fB = yA - fA*xA;
+
+  DefineName(fTheta);
+  // Centers of modules
+  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);
+  fOK1.Set(xk1,yk1);
+
+  Double_t xk2 = fOK.X() + kk1*TMath::Sin(fTheta);
+  Double_t yk2 = fOK.Y() - kk1*TMath::Cos(fTheta);
+  fOK2.Set(xk2,yk2);
+}
+
+void AliEMCALShishKebabTrd1Module::DefineFirstModule()
+{
+  fOK.Set(fga2/2., fgr + fgb/2.); // position the center of module vs o
+
+  // parameters of right line : y = A*z + B in system where zero point is IP.
+  fThetaA = fTheta - fgangle/2.;
+  fA      = TMath::Tan(fThetaA);
+  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());
+  fOK2.Set(fOK.X() + kk1, fOK.Y());
+
+  TObject::SetUniqueID(1); //
+}
+
+void AliEMCALShishKebabTrd1Module::DefineName(const double theta)
+{
+  char name[100];
+  // sprintf(name,"theta_%5.2f",theta*180./TMath::Pi());
+  sprintf(name,"%2i(%5.2f)", TObject::GetUniqueID(), theta*180./TMath::Pi());
+  SetName(name);
+}
+
+Bool_t AliEMCALShishKebabTrd1Module::GetParameters()
+{
+  fgGeometry = AliEMCALGeometry::GetInstance();
+  TString sn(fgGeometry->GetName()); // 2-Feb-05
+  sn.ToUpper();
+  if(!fgGeometry) {
+    Warning("GetParameters()"," No geometry ");
+    return kFALSE; 
+  }
+
+  fga        = (Double_t)fgGeometry->GetEtaModuleSize();
+  fgb        = (Double_t)fgGeometry->GetLongModuleSize();
+  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());
+  Print(0);
+  return kTRUE;
+}
+
+// service methods
+void AliEMCALShishKebabTrd1Module::Print(const int pri) const
+{
+  if(pri>=0) {
+    Info("Print()", "\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) {
+      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());
+    }
+  }
+}
diff --git a/EMCAL/AliEMCALShishKebabTrd1Module.h b/EMCAL/AliEMCALShishKebabTrd1Module.h
new file mode 100644 (file)
index 0000000..e23240a
--- /dev/null
@@ -0,0 +1,75 @@
+#ifndef ALIEMCALSHISHKEBABTRD1MODULE_H
+#define ALIEMCALSHISHKEBABTRD1MODULE_H
+
+/* Copyright(c) 1998-2004, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+//*-- Author: Aleksei Pavlinov (WSU)
+// TO DO : create class for Super module Geometry - 4-nov-04 
+
+#include "TNamed.h"
+#include "TMath.h"
+#include "TVector2.h"
+
+class AliEMCALGeometry;
+
+class AliEMCALShishKebabTrd1Module : public TNamed {
+ public:
+  AliEMCALShishKebabTrd1Module(const double theta=TMath::Pi()/2.);
+  AliEMCALShishKebabTrd1Module(AliEMCALShishKebabTrd1Module &leftNeighbor);
+  void Init(const double A,const double B);
+
+  virtual ~AliEMCALShishKebabTrd1Module(void) {}
+  Bool_t GetParameters();
+  void DefineName(const double theta);
+  void DefineFirstModule();
+
+  Double_t GetTheta() const{return fTheta;}
+  Double_t GetThetaInDegree() const {return fTheta*180./TMath::Pi();}
+  TVector2& GetCenterOfModule() {return fOK;}
+  Double_t  GetEtaOfCenterOfModule(){return -TMath::Log(TMath::Tan(fOK.Phi()/2.));}
+
+  Double_t  GetPosX() {return fOK.Y();}
+  Double_t  GetPosZ() {return fOK.X();}
+  Double_t  GetPosXfromR() {return fOK.Y() - fgr;}
+  Double_t  GetA() {return fA;}
+  Double_t  GetB() {return fB;}
+  //  Additional offline staff 
+  TVector2& GetCenterOfCell(const Int_t ieta)
+  { if(ieta<=1) return fOK1;
+    else        return fOK2;}
+  // 
+  Double_t GetTanBetta() {return fgtanBetta;}
+  Double_t Getb()        {return fgb;}
+  // service methods
+  void Print(const int pri=1) const;  // *MENU*
+
+  // geometry info
+  static AliEMCALGeometry *fgGeometry; //!
+  static Double_t fga;        // 2*dx1=2*dy1
+  static Double_t fga2;       // 2*dx2
+  static Double_t fgb;        // 2*dz1
+  static Double_t fgangle;    // ~1 degree
+  static Double_t fgtanBetta; // tan(fgangle/2.)
+  // radius to IP
+  static Double_t fgr;
+
+ protected:
+  TVector2 fOK;     // position the module center x->y; z->x;
+  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
+  // position of towers with differents ieta (1 or 2) -  4-nov-04
+  TVector2 fOK1;
+  TVector2 fOK2;
+
+ public:
+  ClassDef(AliEMCALShishKebabTrd1Module,0) // Turned Shish-Kebab module 
+};
+
+#endif
+/* To do
+ 1. Insert position the center of towers - 2 additional TVector2
+ */
diff --git a/EMCAL/AliEMCALWsuCosmicRaySetUp.cxx b/EMCAL/AliEMCALWsuCosmicRaySetUp.cxx
new file mode 100644 (file)
index 0000000..47df67e
--- /dev/null
@@ -0,0 +1,132 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Wsu Cosmic Ray SetUp                                                     //
+//  This class contains the description of the  Wsu Cosmic Ray SetUp         //
+//  external volume                                                          //
+//                                                                           //
+//Begin_Html
+/*
+<img src="picts/AliEMCALWsuCosmicRaySetUpClass.gif">
+</pre>
+<br clear=left>
+<font size=+2 color=red>
+<p>The responsible person for this module is
+<a href="mailto:pavlinov@physics.wayne.edu">Aleksei Pavlino, WSU</a>.
+</font>
+<pre>
+*/
+//End_Html
+//                                                                           //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include <TVirtualMC.h>
+
+#include "AliEMCALWsuCosmicRaySetUp.h"
+//#include "AliMagF.h"
+#include "AliRun.h"
+
+ClassImp(AliEMCALWsuCosmicRaySetUp)
+//_____________________________________________________________________________
+AliEMCALWsuCosmicRaySetUp::AliEMCALWsuCosmicRaySetUp()
+{
+  //
+  // Default constructor
+  //
+}
+//_____________________________________________________________________________
+AliEMCALWsuCosmicRaySetUp::AliEMCALWsuCosmicRaySetUp(const char *name, const char *title)
+       : AliModule(name,title)
+{
+  //
+  // Standard constructor of the  Wsu Cosmic Ray SetUp external volume
+  //
+  SetMarkerColor(7);
+  SetMarkerStyle(2);
+  SetMarkerSize(0.4);
+}
+//_____________________________________________________________________________
+void AliEMCALWsuCosmicRaySetUp::CreateGeometry()
+{
+  //
+  // Create the geometry of the Alice external body
+  //
+  //Begin_Html
+  /*
+    <img src="picts/AliEMCALWsuCosmicRaySetUpTree.gif">
+  */
+  //End_Html
+
+  Float_t dASUC[3];
+  Int_t *idtmed = fIdtmed->GetArray()+1;
+  int idSC = idtmed[0];
+  //
+  dASUC[0]=50;
+  dASUC[1]=50;
+  dASUC[2]=50;
+  //  TString tmp(GetTitle());
+  gMC->Gsvolu(GetName(),"BOX",idSC, dASUC,3); // WSUC - Wsu Cosmic Ray SetUp
+}
+//_____________________________________________________________________________
+void AliEMCALWsuCosmicRaySetUp::CreateMaterials()
+{
+// Create materials and media
+  Int_t   isxfld = 0;
+  Float_t sxmgmx = 0.;
+  
+  // AIR
+  Float_t aAir[4]={12.0107,14.0067,15.9994,39.948};
+  Float_t zAir[4]={6.,7.,8.,18.};
+  Float_t wAir[4]={0.000124,0.755267,0.231781,0.012827};
+  Float_t dAir = 1.20479E-3;
+  //  Float_t dAir1 = 1.20479E-10;
+  //
+  AliMixture(1,"Air     $",aAir,zAir,dAir,4,wAir);
+  //
+  AliMedium(1,"Air     $",1,0,isxfld,sxmgmx,10,-1,-0.1,0.1 ,-10);
+}
+//_____________________________________________________________________________
+void AliEMCALWsuCosmicRaySetUp::DrawWSUC(float cxy) const
+{
+  //
+  // Draw a view of the Wsu Cosmic Ray SetUp 
+  //
+  // Set everything unseen
+  gMC->Gsatt("*", "seen", -1);
+  // 
+  // Set WSUC mother visible
+  gMC->Gsatt("WSUC","SEEN",1);
+  //
+  // Set the volumes visible
+  //
+  gMC->Gdopt("hide","off");
+
+  gMC->Gdraw("WSUC", 40, 30, 0, 10, 9, cxy, cxy);
+  gMC->Gdhead(1111, "WSU Cosmic Ray Setup ");
+
+  gMC->Gdman(18, 4, "MAN");
+}
+
diff --git a/EMCAL/AliEMCALWsuCosmicRaySetUp.h b/EMCAL/AliEMCALWsuCosmicRaySetUp.h
new file mode 100644 (file)
index 0000000..b3fb8ce
--- /dev/null
@@ -0,0 +1,29 @@
+#ifndef ALIBODY_H
+#define ALIBODY_H
+/* Copyright(c) 1998-2005, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////////////////
+//  Manager class for detector:  AliEMCALWsuCosmicRaySetUp        //
+//   This is the envelop for Alice                                //
+///////////////////////////////////////////////////////////////////
+#include "AliModule.h"
+
+class AliEMCALWsuCosmicRaySetUp : public AliModule {
+public:
+  AliEMCALWsuCosmicRaySetUp();
+  AliEMCALWsuCosmicRaySetUp(const char *name, const char *title);
+  virtual     ~AliEMCALWsuCosmicRaySetUp() {}
+  virtual void  CreateGeometry();
+  virtual void  CreateMaterials();
+  virtual Int_t IsVersion() const {return 0;}
+  void  DrawWSUC(float cxy=0.025) const; // *MENU*
+
+  ClassDef(AliEMCALWsuCosmicRaySetUp,1)  // Class manager for the Wsu Cosmic Ray SetUp
+};
+
+#endif
index 5e17ae6..9d4e2a9 100644 (file)
 
 // --- ROOT system ---
 
-//#include "TPGON.h"
-#include "TTUBS.h"
 #include "TNode.h"
+#include "TBRIK.h"
+#include "TTRD1.h"
+#include "TTRD2.h"
+#include "TTRAP.h"
+#include "TPGON.h"
+#include "TTUBS.h"
 #include "TGeometry.h"
 #include "TVirtualMC.h"
 #include "TArrayI.h"
+#include "TROOT.h"
+#include "TArrayF.h"
+#include "TList.h"
+#include "TVector2.h"
+
+#include "AliEMCALShishKebabModule.h"
+#include "AliEMCALShishKebabTrd1Module.h"
+//#include <TGeant3.h> // can not include - I don't know why
 
 // --- Standard library ---
 
 
 ClassImp(AliEMCALv0)
 
+TArrayF ENVELOP1; // now global for BuildGeometry() - 30-aug-2004
+int idAIR=1599, idPB = 1600, idSC = 1601, idSTEEL = 1603; // global
+Int_t *idtmed=0, idrotm=0;
+double sampleWidth=0.;
+double parEMOD[5], smodPar0=0., smodPar1=0., smodPar2=0.;
+
+
 //______________________________________________________________________
 AliEMCALv0::AliEMCALv0(const char *name, const char *title):
   AliEMCAL(name,title)
 {
   // ctor : title is used to identify the layout
   GetGeometry() ; 
-
+  fShishKebabModules = 0;
 }
 
 //______________________________________________________________________
@@ -67,23 +86,153 @@ void AliEMCALv0::BuildGeometry()
     const Int_t kColorArm1   = kBlue ;
 
     AliEMCALGeometry * geom = GetGeometry() ; 
+    TString gn(geom->GetName());
+    gn.ToUpper(); 
 
     // Define the shape of the Calorimeter 
-    TNode * top = gAlice->GetGeometry()->GetNode("alice") ;
-    new TTUBS("Envelop1", "Tubs that contains arm 1", "void", 
-             geom->GetEnvelop(0),     // rmin 
-             geom->GetEnvelop(1) +30 ,// rmax
+    TNode * top = gAlice->GetGeometry()->GetNode("alice") ; // See AliceGeom/Nodes
+    TNode * envelopNode = 0;
+    char *envn = "Envelop1";
+    if(!gn.Contains("SHISH") || gn.Contains("TRD2")){
+      new TTUBS(envn, "Tubs that contains arm 1", "void", 
+             geom->GetEnvelop(0) -10, // rmin 
+             geom->GetEnvelop(1) +40 ,// rmax
              geom->GetEnvelop(2)/2.0, // half length in Z
              geom->GetArm1PhiMin(),   // minimum phi angle
              geom->GetArm1PhiMax()    // maximum phi angle
        );
+      top->cd();
+      envelopNode = new TNode(envn, "Arm1 Envelop", "Envelop1", 0., 0., 0., "") ;
+    } else {
+      if(gn.Contains("WSUC")) {
+        envelopNode = BuildGeometryOfWSUC();
+      } else { // Shish-kebab now for compact, twist and TRD1 cases (ALIC)
+        envn="Envelop2";
+        TPGON *pgon = new TPGON(envn, "PGON that contains arm 1", "void", 
+        geom->GetArm1PhiMin(),geom->GetArm1PhiMax()-geom->GetArm1PhiMin(),geom->GetNPhiSuperModule(), 2);
+      // define section
+        pgon->DefineSection(0, ENVELOP1[4],  ENVELOP1[5], ENVELOP1[6]);
+        pgon->DefineSection(1, ENVELOP1[7],  ENVELOP1[5], ENVELOP1[6]);
+        top->cd();
+        envelopNode = new TNode(envn, "Arm1 Envelop2", envn, 0., 0., 0., "") ;
+      }
+    }
+                               
+    envelopNode->SetLineColor(kColorArm1) ;
+    fNodes->Add(envelopNode);
+}
+
+TNode *AliEMCALv0::BuildGeometryOfWSUC()
+{ // June 8, 2005; see directory geant3/TGeant3/G3toRoot.cxx
+  // enum EColor { kWhite, kBlack, kRed, kGreen, kBlue, kYellow, kMagenta, kCyan } - see $ROOTSYS/include/Gtypes.h
+   AliEMCALGeometry * g = GetGeometry(); 
+   TNode * top = gAlice->GetGeometry()->GetNode("alice") ; // See AliceGeom/Nodes
+   top->cd();
+
+   TNode *envelopNode = 0;
+   char *name = "";
+   /*
+    name = "WSUC";
+   new TBRIK(name, "WSUC(XEN1 in Geant)","void",ENVELOP1[0],ENVELOP1[1],ENVELOP1[2]);
+    envelopNode = new TNode(name, "envelope for WSUC", name, 0., 0., 0., "");
+   envelopNode->SetVisibility(0);
+   */
+
+   TNode *emod=0, *scmx=0;
+   name = "SMOD"; // super module
+   new TBRIK(name, "SMOD(SMOD in Geant)","void", smodPar0,smodPar1,smodPar2);
+   if(envelopNode) envelopNode->cd();
+   TNode *smod = new TNode(name, "SMOD", name, 0., 0., 0., "");   
+   smod->SetLineColor(kBlue) ;
+   if(envelopNode==0) envelopNode = smod;
+
+   name = "EMOD"; // see CreateEMOD 
+   TTRD1 *EMOD = new TTRD1(name, "EMOD(EMOD in Geant)","void", float(parEMOD[0]),
+   float(parEMOD[1]),float(parEMOD[2]),float(parEMOD[3]));
+
+   // SCMX = EMOD/4 for simplicity of drawing
+   name = "SCMX";
+   Float_t dz=0.,theta=0.,phi=0.,h1=0.,bl1=0.,tl1=0.,alpha1=0.,h2=0.,bl2=0.,tl2=0.,alpha2=0.;
+   h1     = EMOD->GetDy()/2.;
+   bl1    = EMOD->GetDx()/2.;
+   tl1    = bl1;
+   alpha1 = 0.;
+   h2     = EMOD->GetDy()/2.;
+   bl2    = EMOD->GetDx2()/2.;
+   tl2    = bl2;
+   alpha2 = 0.;
+
+   dz       = EMOD->GetDz();
+   double dr = TMath::Sqrt((h2-h1)*(h2-h1)+(bl2-bl1)*(bl2-bl1));
+   theta    = TMath::ATan2(dr,2.*dz) * TMath::RadToDeg(); 
+   phi      = 180.;
+
+   TTRAP *SCMX = new TTRAP(name, "SCMX(SCMX as in Geant)","void",
+            dz,theta,phi, h1,bl1,tl1,alpha1, h2,bl2,tl2,alpha2);
+   //   SCMX->Dump();
+   Float_t xShiftSCMX = (EMOD->GetDx() + EMOD->GetDx2())/4.;
+   Float_t yShiftSCMX = EMOD->GetDy()/2.;
+   printf(" xShiftSCMX %7.4f yShiftSCMX %7.4f \n",xShiftSCMX,yShiftSCMX);
+
+   name = "EMOD"; // see CreateEMOD 
+   smod->cd();
+   
+   AliEMCALShishKebabTrd1Module *mod=0;
+   Double_t  angle=90., xpos=0.,ypos=0.,zpos=0., xposSCMX=0.,yposSCMX=0.,zposSCMX=0.;
+   char rtmn[100], rtmt[100];
+   TRotMatrix *rtm=0, *rtmSCMX=0;
+   int numEmod=0;
+   for(int iz=0; iz<g->GetNZ(); iz++) {
+     mod   = (AliEMCALShishKebabTrd1Module*)fShishKebabModules->At(iz);
+     zpos = mod->GetPosZ()      - smodPar2;
+     ypos = mod->GetPosXfromR() - smodPar1;
+
+     angle = mod->GetThetaInDegree();
+     sprintf(rtmn,"rmEmod%5.1f",angle);
+     sprintf(rtmt,"rotation matrix for EMOD, iz=%i, angle = %6.3f",iz, angle);
+     if(iz==0) rtm = new TRotMatrix(rtmn, rtmt,0.,0., 90.,0., 90.,90.); // z'(x); y'(y); x'(z)
+     else      rtm = new TRotMatrix(rtmn, rtmt,90.-angle,270., 90.0,0.0, angle,90.);
 
-    // Place the Node
-    top->cd();
-    TNode * envelop1node = new TNode("Envelop1", "Arm1 Envelop", "Envelop1"
-                                    ,0., 0., 0., "") ;
-    envelop1node->SetLineColor(kColorArm1) ;
-    fNodes->Add(envelop1node) ;
+     TGeometry *tg = gAlice->GetGeometry();
+     for(int ix=0; ix<g->GetNPhi(); ix++) { // flat in phi
+       xpos = g->GetPhiModuleSize()*(2*ix+1 - g->GetNPhi())/2.;
+       sprintf(rtmt,"EMOD, iz %i, ix %i, angle %6.3f",iz,ix, angle);
+       TString namNode=name;
+       namNode += numEmod++;
+       smod->cd();
+       emod = new TNode(namNode.Data(), rtmt, (TShape*)EMOD, xpos,ypos,zpos,rtm);   
+       //       emod->SetLineColor(kGreen) ;
+       emod->SetVisibility(0); // SCMX will bi visible 
+       if(SCMX) { // 4(2x2) sensetive volume inside EMOD
+         emod->cd();
+         zposSCMX = 0.;
+         for(int jy=0; jy<2; jy++){ // division on y
+           yposSCMX = yShiftSCMX *(2*jy - 1); 
+           for(int jx=0; jx<2; jx++){ // division on x
+             Double_t theta1=90.,phi1=0., theta2=90.,phi2=90., theta3=0.,phi3=0 ; 
+             xposSCMX = xShiftSCMX *(2*jx - 1);              
+             namNode = "SCMX";
+             namNode += jy;
+             namNode += jx;
+             sprintf(rtmn,"rm%s",namNode.Data());
+             sprintf(rtmt,"rotation matrix for %s inside EMOD",namNode.Data());
+             rtmSCMX = tg->GetRotMatrix(rtmn);
+             if(jx == 1) {
+               phi1 = 180.; // x' = -x
+               phi2 = 270.; // y' = -y
+             }
+             if(rtmSCMX == 0) rtmSCMX = new TRotMatrix(rtmn,rtmt, theta1,phi1, theta2,phi2, theta3,phi3);
+             sprintf(rtmt,"%s inside %s", namNode.Data(), emod->GetName());
+             scmx = new TNode(namNode.Data(), rtmt, (TShape*)SCMX, xposSCMX,yposSCMX,zposSCMX,rtmSCMX);   
+             scmx->SetLineColor(kMagenta);
+           }
+         }
+       }
+     }
+   }
+   //   emod->Draw(); // for testing
+
+   return envelopNode;
 }
 
 //______________________________________________________________________
@@ -118,44 +267,82 @@ void AliEMCALv0::CreateGeometry()
     Float_t *dum=0;
 
     AliEMCALGeometry * geom = GetGeometry() ; 
+    TString gn(geom->GetName());
+    gn.ToUpper(); 
 
     if(!(geom->IsInitialized())){
        Error("CreateGeometry","EMCAL Geometry class has not been set up.");
     } // end if
 
     // Get pointer to the array containing media indices
-    Int_t *idtmed = fIdtmed->GetArray() - 1599 ;
+    idtmed = fIdtmed->GetArray() - 1599 ;
 
-    Int_t idrotm = 1;
-    AliMatrix(idrotm, 90.0, 0., 90.0, 90.0, 0.0, 0.0) ;
+    idrotm = 1;
+    //  gMC->Matrix(nmat, theta1, phi1, theta2, phi2, theta3, phi3) - see AliModule
+    AliMatrix(idrotm, 90.0, 0., 90.0, 90.0, 0.0, 0.0) ; 
 
     // Create the EMCAL Mother Volume (a polygone) within which to place the Detector and named XEN1 
 
     Float_t envelopA[10];
-    envelopA[0] = geom->GetArm1PhiMin();                         // minimum phi angle
-    envelopA[1] = geom->GetArm1PhiMax() - geom->GetArm1PhiMin(); // angular range in phi
-    envelopA[2] = geom->GetNPhi();                               // number of sections in phi
-    envelopA[3] = 2;                                             // 2 z coordinates
-    envelopA[4] = geom->ZFromEtaR(geom->GetEnvelop(1),
+    if(gn.Contains("TRD2")) { // TUBS
+       envelopA[0] = geom->GetEnvelop(0) - 10.; // rmin 
+       envelopA[1] = geom->GetEnvelop(1) + 12.; // rmax
+       //       envelopA[2] = geom->ZFromEtaR(geom->GetEnvelop(1), geom->GetArm1EtaMin());
+       envelopA[2] = 390.; // 6-feb-05
+       envelopA[3] = geom->GetArm1PhiMin();
+       envelopA[4] = geom->GetArm1PhiMax();
+       gMC->Gsvolu("XEN1", "TUBS", idtmed[1599], envelopA, 5) ;   // Tubs filled with air 
+       ENVELOP1.Set(5, envelopA);
+    // Position the EMCAL Mother Volume (XEN1) in Alice (ALIC)  
+       gMC->Gspos("XEN1", 1, "ALIC", 0.0, 0.0, 0.0, idrotm, "ONLY") ;
+    } else if(gn.Contains("TRD1") && gn.Contains("WSUC") ) { // TRD1 for WSUC facility
+      // 17-may-05 - just BOX
+      envelopA[0] = 26;
+      envelopA[1] = 15;
+      envelopA[2] = 30;
+      gMC->Gsvolu("XEN1", "BOX", idtmed[idSC], envelopA, 3) ;
+      ENVELOP1.Set(3);
+      for(int i=0; i<3; i++) ENVELOP1[i] = envelopA[i]; // 23-may-05  
+      // Position the EMCAL Mother Volume (XEN1) in WSUC  
+      gMC->Gspos("XEN1", 1, "WSUC", 0.0, 0.0, 0.0, idrotm, "ONLY") ;
+    } else { 
+      envelopA[0] = geom->GetArm1PhiMin();                         // minimum phi angle
+      envelopA[1] = geom->GetArm1PhiMax() - geom->GetArm1PhiMin(); // angular range in phi
+      envelopA[2] = geom->GetNPhi();                               // number of sections in phi
+      envelopA[3] = 2;                                             // 2 z coordinates
+      envelopA[4] = geom->ZFromEtaR(geom->GetEnvelop(1),
                                   geom->GetArm1EtaMin());       // z coordinate 1
     //add some padding for mother volume
-    envelopA[5] = geom->GetEnvelop(0) ;                          // rmin at z1
-    envelopA[6] = geom->GetEnvelop(1) ;                          // rmax at z1
-    envelopA[7] = geom->ZFromEtaR(geom->GetEnvelop(1),
+      envelopA[5] = geom->GetEnvelop(0) ;                          // rmin at z1
+      envelopA[6] = geom->GetEnvelop(1) ;                          // rmax at z1
+      envelopA[7] = geom->ZFromEtaR(geom->GetEnvelop(1),
                                  geom->GetArm1EtaMax());        // z coordinate 2
-    envelopA[8] = envelopA[5] ;                                  // radii are the same.
-    envelopA[9] = envelopA[6] ;                                  // radii are the same.
+      envelopA[8] = envelopA[5] ;                                  // radii are the same.
+      envelopA[9] = envelopA[6] ;                                  // radii are the same.
 
-    gMC->Gsvolu("XEN1", "PGON ", idtmed[1599], envelopA, 10) ;   // Polygone filled with air 
+      if(gn.Contains("SHISH")) envelopA[2] = geom->GetNPhiSuperModule();
 
+      gMC->Gsvolu("XEN1", "PGON ", idtmed[1599], envelopA, 10) ;   // Polygone filled with air 
+      ENVELOP1.Set(10, envelopA);
+      if (gDebug==2) {
+        printf("CreateGeometry: XEN1 = %f, %f\n", envelopA[5], envelopA[6]); 
+        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") ;
+    }
+
+    if(gn.Contains("SHISH")){
+      // Compact and twist
+      CreateShishKebabGeometry();
+      return;
+    }
 
-    gMC->Gspos("XEN1", 1, "ALIC", 0.0, 0.0, 0.0, idrotm, "ONLY") ;
-    
     if (AliLog::GetGlobalDebugLevel()>=2) {
       printf("CreateGeometry: XEN1 = %f, %f\n", envelopA[5], envelopA[6]); 
       printf("CreateGeometry: XU0 = %f, %f\n", envelopA[5], envelopA[6]); 
     }
+
     // Create mini-envelopes which will contain the Tower scintillator-radiator
     
     TString label ;
@@ -343,3 +530,788 @@ void AliEMCALv0::Init(void)
     printf(message.Data() ) ; 
   }
 }
+
+// 24-aug-04 by PAI
+void AliEMCALv0::CreateShishKebabGeometry()
+{  // TWIST, TRD1 and TRD2 
+  AliEMCALGeometry * g = GetGeometry(); 
+  TString gn(g->GetName()); gn.ToUpper(); 
+  // see AliModule::fIdtmed
+  //  idtmed = fIdtmed->GetArray() - 1599 ; // see AliEMCAL::::CreateMaterials()
+  //  int idAIR=1599, idPB = 1600, idSC = 1601, idSTEEL = 1603;
+  //  idAL = 1602;
+  Double_t par[10], xpos=0., ypos=0., zpos=0.;
+
+  CreateSmod("XEN1"); // 18-may-05 
+
+  CreateEmod("SMOD","EMOD"); // 18-may-95
+
+  // Sensitive SC  (2x2 tiles)
+  double parSCM0[5], *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.;
+    par[1] = par[2] = g->GetPhiTileSize();   // Symetric case
+    gMC->Gsvolu("SCM0", "BOX", idtmed[idSC], par, 3); // 2x2 tiles
+    gMC->Gspos("SCM0", 1, "EMOD", 0., 0., 0., 0, "ONLY") ;
+    // Division to tile size
+    gMC->Gsdvn("SCM1","SCM0", g->GetNPHIdiv(), 2); // y-axis
+    gMC->Gsdvn("SCM2","SCM1", g->GetNETAdiv(), 3); // z-axis
+  // put LED to the SCM2 
+    par[0] = g->GetECPbRadThick()/2.;
+    par[1] = par[2] = g->GetPhiTileSize()/2.; // Symetric case
+    gMC->Gsvolu("PBTI", "BOX", idtmed[idPB], par, 3);
+
+    printf(" Pb tiles \n");
+    int nr=0;
+    ypos = zpos = 0.0;
+    xpos = -sampleWidth*g->GetNECLayers()/2. + g->GetECPbRadThick()/2.;
+    for(int ix=0; ix<g->GetNECLayers(); ix++){
+      gMC->Gspos("PBTI", ++nr, "SCM2", xpos, ypos, zpos, 0, "ONLY") ;
+    //    printf(" %i xpos %f \n", ix+1, xpos);
+      xpos += sampleWidth;
+    } 
+    printf(" Number of Pb tiles in SCM2 %i \n", nr);
+  } else if(gn.Contains("TRD1")) { // TRD1 - 30-sep-04
+    if(gn.Contains("MAY05")){
+      Double_t dzTmp = g->GetFrontSteelStrip()+g->GetPassiveScintThick();
+      parSCM0[0] = parEMOD[0] + tanTmp*dzTmp; // dx1
+      parSCM0[1] = parEMOD[1];                // dx2
+      parSCM0[2] = parEMOD[2];                // dy
+      for(int i=0; i<3; i++) parSCM0[i] -= g->GetLateralSteelStrip();
+      parSCM0[3] = parEMOD[3] - dzTmp/2.; // dz
+
+      gMC->Gsvolu("SCM0", "TRD1", idtmed[idAIR], 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
+      for(int i=0; i<3; i++) parSCM0[i] = parEMOD[i] - wallThickness;
+      parSCM0[3] = parEMOD[3];
+      gMC->Gsvolu("SCM0", "TRD1", idtmed[idAIR], parSCM0, 4);
+      gMC->Gspos("SCM0", 1, "EMOD", 0., 0., 0., 0, "ONLY") ;
+    }
+
+    if(g->GetNPHIdiv()==2 && g->GetNETAdiv()==2) {
+    // Division to tile size - 1-oct-04
+      printf("<I> Divide SCM0 on y-axis %i\n", g->GetNETAdiv());
+      gMC->Gsdvn("SCMY","SCM0", g->GetNETAdiv(), 2); // y-axis
+    // Trapesoid 2x2
+      parTRAP[0] = parSCM0[3];    // dz
+      parTRAP[1] = TMath::ATan2((parSCM0[1]-parSCM0[0])/2.,2.*parSCM0[3])*180./TMath::Pi(); // theta
+      parTRAP[2] = 0.;           // phi
+      // bottom
+      parTRAP[3] = parSCM0[2]/2.; // H1
+      parTRAP[4] = parSCM0[0]/2.; // BL1
+      parTRAP[5] = parTRAP[4];    // TL1
+      parTRAP[6] = 0.0;           // ALP1
+      // top
+      parTRAP[7] = parSCM0[2]/2.; // H2
+      parTRAP[8] = parSCM0[1]/2.; // BL2
+      parTRAP[9] = parTRAP[8];    // TL2
+      parTRAP[10]= 0.0;           // ALP2
+      printf(" ** TRAP ** \n");
+      for(int i=0; i<11; i++) printf(" par[%2.2i] %9.4f\n", i, parTRAP[i]);
+
+      gMC->Gsvolu("SCMX", "TRAP", idtmed[idSC], parTRAP, 11);
+      xpos = +(parSCM0[1]+parSCM0[0])/4.;
+      gMC->Gspos("SCMX", 1, "SCMY", xpos, 0.0, 0.0, 0, "ONLY") ;
+
+      // Using rotation because SCMX should be the same due to Pb tiles
+      xpos = -xpos; 
+      AliMatrix(idrotm, 90.0,180., 90.0, 270.0, 0.0,0.0) ;
+      gMC->Gspos("SCMX", 2, "SCMY", xpos, 0.0, 0.0, idrotm, "ONLY");
+    // put LED to the SCM0 
+      AliEMCALShishKebabTrd1Module *mod = (AliEMCALShishKebabTrd1Module*)fShishKebabModules->At(0);
+      gMC->Gsvolu("PBTI", "BOX", idtmed[idPB], dummy, 0);
+
+      par[1] = parSCM0[2]/2;            // y 
+      par[2] = g->GetECPbRadThick()/2.; // z
+
+      int nr=0;
+      ypos = 0.0;
+      zpos = -sampleWidth*g->GetNECLayers()/2. + g->GetECPbRadThick()/2.;
+      double xCenterSCMX =  (parTRAP[4] +  parTRAP[8])/2.;
+      if(!gn.Contains("NOPB")) { // for testing - 11-jul-05
+        printf(" Pb tiles \n");
+        for(int iz=0; iz<g->GetNECLayers(); iz++){
+          par[0] = (parSCM0[0] + mod->GetTanBetta()*sampleWidth*iz)/2.;
+          xpos   = par[0] - xCenterSCMX;
+          gMC->Gsposp("PBTI", ++nr, "SCMX", xpos, ypos, zpos, 0, "ONLY", par, 3) ;
+          printf(" %i xpos %f zpos %f par[0] %f \n", iz+1, xpos, zpos, par[0]);
+          zpos += sampleWidth;
+        } 
+        printf(" Number of Pb tiles in SCMX %i \n", nr);
+      }
+    } else if(g->GetNPHIdiv()==3 && g->GetNETAdiv()==3) {
+      Trd1Tower3X3(parSCM0);
+    } else if(g->GetNPHIdiv()==4 && g->GetNETAdiv()==4) {
+      Trd1Tower4X4();
+    }
+  } else if(gn.Contains("TRD2")) {    // TRD2 - 14-jan-05
+    //    Scm0InTrd2(g, parEMOD, parSCM0); // First dessin 
+    PbmoInTrd2(g, parEMOD, parSCM0); // Second dessin 
+  }
+}
+
+void AliEMCALv0::CreateSmod(const char* mother)
+{ // 18-may-05; mother="XEN1"; child="SMOD" or "SMON" and "SMOP"("TRD2" case)
+  AliEMCALGeometry * g = GetGeometry(); 
+  TString gn(g->GetName()); gn.ToUpper();
+
+  Double_t par[1], parTubs[5], xpos=0., ypos=0., zpos=0., rpos=0., dphi=0., phi=0.0, phiRad=0.;
+  //  ===== define Super Module from air - 14x30 module ==== ;
+  sampleWidth = double(g->GetECPbRadThick()+g->GetECScintThick());
+  printf("\n ## Super Module | sampleWidth %5.3f ## \n", sampleWidth);
+  par[0] = g->GetShellThickness()/2.;
+  par[1] = g->GetPhiModuleSize()*g->GetNPhi()/2.; 
+  par[2] = g->GetEtaModuleSize()*15.; 
+  idrotm=0;
+  int nphism = g->GetNumberOfSuperModules()/2; // 20-may-05
+  if(nphism>0) {
+    dphi = (g->GetArm1PhiMax() - g->GetArm1PhiMin())/nphism;
+    rpos = (g->GetEnvelop(0) + g->GetEnvelop(1))/2.;
+    printf(" rpos %8.2f \n", rpos);
+  }
+
+  if (gn.Contains("TRD2")) { // tubs - 27-jan-05
+    parTubs[0] = g->GetTubsR();                       // rmin
+    parTubs[1] = parTubs[0] + g->GetShellThickness(); // rmax ?? 
+    parTubs[2] = 380./2.;                             // DZ half length in z; 11-oct-04 - for 26 division
+    parTubs[3] = -dphi/2.;                            // PHI1 starting angle of the segment;
+    parTubs[4] = +dphi/2.;                            // PHI2 ending angle of the segment;
+
+    gMC->Gsvolu("SMOP", "TUBS", idtmed[idAIR], parTubs, 5); // pozitive Z
+    gMC->Gsvolu("SMON", "TUBS", idtmed[idAIR], parTubs, 5); // negative Z
+
+    printf(" SMOP,N ** TUBS **\n"); 
+    printf("tmed %i | Rmin %7.2f Rmax %7.2f dz %7.2f phi1,2 (%7.2f,%7.2f)\n", 
+    idtmed[idAIR], parTubs[0],parTubs[1],parTubs[2], parTubs[3],parTubs[4]);
+    // have to add 1 cm plastic before EMOD - time solution 
+  } else if(gn.Contains("WSUC")) {
+    par[0] = g->GetPhiModuleSize()*g->GetNPhi()/2.; 
+    par[1] = g->GetShellThickness()/2.;
+    par[2] = g->GetEtaModuleSize()*g->GetNZ()/2. + 5; 
+
+    gMC->Gsvolu("SMOD", "BOX", idtmed[idAIR], par, 3);
+
+    printf("SMOD in WSUC : tmed %i | dx %7.2f dy %7.2f dz %7.2f (SMOD, BOX)\n", 
+    idtmed[idAIR], par[0],par[1],par[2]);
+    smodPar0 = par[0]; 
+    smodPar1 = par[1];
+    smodPar2 = par[2];
+    nphism   =  g->GetNumberOfSuperModules();
+  } else {
+    if     (gn.Contains("TWIST")) {
+      par[2] += 0.4;      // for 27 division
+    } else if(gn.Contains("TRD")) {
+      par[2]  = 350./2.; // 11-oct-04 - for 26 division
+    }
+  //  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]);
+    smodPar0 = par[0]; 
+    smodPar2 = par[2];
+  // Steel plate
+    if(g->GetSteelFrontThickness() > 0.0) { // 28-mar-05
+      par[0] = g->GetSteelFrontThickness()/2.;
+      gMC->Gsvolu("STPL", "BOX", idtmed[idSTEEL], par, 3);
+      printf("tmed %i | dx %7.2f dy %7.2f dz %7.2f (STPL) \n", idtmed[idSTEEL], par[0],par[1],par[2]);
+      xpos = -(g->GetShellThickness() - g->GetSteelFrontThickness())/2.;
+      gMC->Gspos("STPL", 1, "SMOD", xpos, 0.0, 0.0, 0, "ONLY") ;
+    }
+  }
+
+  int nr=0, i0=0;
+  if(gn.Contains("TEST")) {nphism = 1;} // just only 2 super modules;
+
+  // Turn whole super module
+  int TurnSupMod = 1; // should be ONE; for testing =0
+  for(int i=i0; i<nphism; i++) {
+    if (gn.Contains("TRD2")) {      // tubs - 27-jan-05
+      if(i==i0) {
+        printf("** TRD2 ** ");
+        if(TurnSupMod==1) printf(" No 3 degree rotation !!! ");
+        printf("\n");
+      }
+      Double_t phic=0., phicRad=0.; // phi angle of arc center
+      phic    = g->GetArm1PhiMin() + dphi*(2*i+1)/2.; //
+      phicRad = phic*TMath::DegToRad();
+      phi     = phic - g->GetTubsTurnAngle();
+      phiRad  = phi*TMath::DegToRad();
+      if(TurnSupMod==1) {
+        TVector2  vc;     // position of super module center
+        vc.SetMagPhi(parTubs[0], phicRad);
+        TVector2  vcTurn; // position of super module center with turn
+        vcTurn.SetMagPhi(parTubs[0], phiRad);
+        TVector2 vcShift = vc - vcTurn;
+        phic = phi;
+
+        xpos = vcShift.X();
+        ypos = vcShift.Y();
+      } else { // 1-mar-05 ; just for testing - no turn od SMOD; looks good
+        xpos = ypos = 0.0;
+      }
+      zpos = parTubs[2];
+      AliMatrix(idrotm, 90.0, phic, 90.0, 90.0+phic, 0.0, 0.0);
+
+      gMC->Gspos("SMOP", ++nr, mother, xpos, ypos, zpos, idrotm, "ONLY") ;
+      printf("SMOP %2i | %2i idrotm %3i phi %6.1f(%5.3f) xpos %7.2f ypos %7.2f zpos %7.2f \n", 
+      i, nr, idrotm, phic, phicRad, xpos, ypos, zpos);
+      printf(" phiy(90+phic)  %6.1f \n", 90. + phic);
+
+      if(!gn.Contains("TEST1") && g->GetNumberOfSuperModules() > 1){
+       //        double  phiy = 90. + phic + 180.;
+       //        if(phiy>=360.) phiy -= 360.;
+       //        printf(" phiy  %6.1f \n", phiy);
+       //        AliMatrix(idrotm, 90.0, phic, 90.0, phiy, 180.0, 0.0);
+        gMC->Gspos("SMON", nr, mother, xpos, ypos, -zpos, idrotm, "ONLY") ;
+        printf("SMON %2i | %2i idrotm %3i phi %6.1f(%5.3f) xpos %7.2f ypos %7.2f zpos %7.2f \n", 
+        i, nr, idrotm, phic, phicRad, xpos, ypos, -zpos);
+      }
+   } else if(gn.Contains("WSUC")) {
+     xpos = ypos = zpos = 0.0;
+     idrotm = 0;
+     gMC->Gspos("SMOD", 1, mother, xpos, ypos, zpos, idrotm, "ONLY") ;
+     printf(" idrotm %3i phi %6.1f(%5.3f) xpos %7.2f ypos %7.2f zpos %7.2f \n", 
+     idrotm, phi, phiRad, xpos, ypos, zpos);
+     nr++;
+   } else {
+      phi    = g->GetArm1PhiMin() + dphi*(2*i+1)/2.; // phi= 70, 90, 110, 130, 150, 170
+      phiRad = phi*TMath::Pi()/180.;
+
+      AliMatrix(idrotm, 90.0, phi, 90.0, 90.0+phi, 0.0, 0.0);
+
+      xpos = rpos * TMath::Cos(phiRad);
+      ypos = rpos * TMath::Sin(phiRad);
+      zpos = smodPar2; // 21-sep-04
+
+      // 1th module in z-direction;
+      gMC->Gspos("SMOD", ++nr, mother, xpos, ypos, zpos, idrotm, "ONLY") ;
+      printf(" %2i idrotm %3i phi %6.1f(%5.3f) xpos %7.2f ypos %7.2f zpos %7.2f \n", 
+      nr, idrotm, phi, phiRad, xpos, ypos, zpos);
+      // 2th module in z-direction;
+      if(gn.Contains("TWIST") || gn.Contains("TRD")) {
+      // turn arround X axis; 0<phi<360
+        double phiy = 90. + phi + 180.;
+        if(phiy>=360.) phiy -= 360.;
+        AliMatrix(idrotm, 90.0, phi, 90.0, phiy, 180.0, 0.0);
+        gMC->Gspos("SMOD", ++nr, mother, xpos, ypos, -zpos, idrotm, "ONLY");
+        printf(" %2i idrotm %3i phiy %6.1f  xpos %7.2f ypos %7.2f zpos %7.2f \n", 
+        nr, idrotm, phiy, xpos, ypos, -zpos);
+      } else {
+        gMC->Gspos("SMOD", ++nr, mother, xpos, ypos, -zpos, idrotm, "ONLY");
+      }
+    }
+  }
+  printf(" Number of Super Modules %i \n", nr);
+}
+
+void AliEMCALv0::CreateEmod(const char* mother, const char* child)
+{ // 17-may-05; mother="SMOD"; child="EMOD"
+  AliEMCALGeometry * g = GetGeometry(); 
+  TString gn(g->GetName()); gn.ToUpper(); 
+  // Module definition
+  Double_t par[10], parTubs[5], xpos=0., ypos=0., zpos=0., rpos=0.;
+  Double_t parSCPA[5], zposSCPA=0.; // passive SC - 13-MAY-05, TRD1 case
+  Double_t trd1Angle = g->GetTrd1Angle()*TMath::DegToRad(), tanTrd1 = TMath::Tan(trd1Angle/2.);
+  Double_t tanTrd2y  = TMath::Tan(g->GetTrd2AngleY()*TMath::DegToRad()/2.);
+  int nr=0;
+  idrotm=0;
+  if(!gn.Contains("TRD")) { // standard module
+    par[0] = (sampleWidth*g->GetNECLayers())/2.; 
+    par[1] = par[2] = g->GetPhiModuleSize()/2.;
+    gMC->Gsvolu(child, "BOX", idtmed[idSTEEL], par, 3);
+
+  } else if (gn.Contains("TRD1")){ // TRD1 system coordinate iz differnet
+    parEMOD[0] = g->GetEtaModuleSize()/2.;   // dx1
+    parEMOD[1] = g->Get2Trd1Dx2()/2.;        // dx2
+    parEMOD[2] = g->GetPhiModuleSize()/2.;;  // dy
+    parEMOD[3] = g->GetLongModuleSize()/2.;  // dz
+    gMC->Gsvolu(child, "TRD1", idtmed[idSTEEL], parEMOD, 4);
+    if(gn.Contains("WSUC") || gn.Contains("MAY05")){
+      parSCPA[0] = g->GetEtaModuleSize()/2. + tanTrd1*g->GetFrontSteelStrip();   // dx1
+      parSCPA[1] = parSCPA[0]               + tanTrd1*g->GetPassiveScintThick(); // dx2
+      parSCPA[2] = g->GetPhiModuleSize()/2.;     // dy
+      parSCPA[3] = g->GetPassiveScintThick()/2.; // dz
+      gMC->Gsvolu("SCPA", "TRD1", idtmed[idSC], parSCPA, 4);
+      zposSCPA   = -parEMOD[3] + g->GetFrontSteelStrip() + g->GetPassiveScintThick()/2.;
+      gMC->Gspos ("SCPA", ++nr, child, 0.0, 0.0, zposSCPA, 0, "ONLY");
+    }
+  } else if (gn.Contains("TRD2")){ // TRD2 as for TRD1 - 27-jan-05
+    parEMOD[0] = g->GetEtaModuleSize()/2.;   // dx1
+    parEMOD[1] = g->Get2Trd1Dx2()/2.;        // dx2
+    parEMOD[2] = g->GetPhiModuleSize()/2.;   // dy1
+    parEMOD[3] = parEMOD[2] + tanTrd2y*g->GetLongModuleSize();// dy2
+    parEMOD[4] = g->GetLongModuleSize()/2.;  // dz
+    gMC->Gsvolu(child, "TRD2", idtmed[idSTEEL], parEMOD, 5);
+  }
+
+  nr   = 0;
+  if(gn.Contains("TWIST")) { // 13-sep-04
+    fShishKebabModules = new TList;
+    AliEMCALShishKebabModule *mod=0, *mTmp; // current module
+    for(int iz=0; iz<g->GetNZ(); iz++) {
+      //for(int iz=0; iz<4; iz++) {
+      if(iz==0) {
+        mod    = new AliEMCALShishKebabModule();
+        idrotm = 0;
+      } else {
+        mTmp = new AliEMCALShishKebabModule(*mod);
+        mod  = mTmp;
+        Double_t  angle = mod->GetThetaInDegree();
+        AliMatrix(idrotm, angle,0., 90.0,90.0, 90.-angle, 180.);
+      }
+      fShishKebabModules->Add(mod);
+
+      xpos = mod->GetPosXfromR() + g->GetSteelFrontThickness() - smodPar0;
+      zpos = mod->GetPosZ() - smodPar2;
+      for(int iy=0; iy<g->GetNPhi(); iy++) {
+        ypos = g->GetPhiModuleSize()*(2*iy+1 - g->GetNPhi())/2.;
+        gMC->Gspos(child, ++nr, mother, xpos, ypos, zpos, idrotm, "ONLY") ;
+      //printf(" %2i xpos %7.2f ypos %7.2f zpos %7.2f \n", nr, xpos, ypos, zpos);
+      }
+    }    
+  } 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)
+    fShishKebabModules = new TList;
+    AliEMCALShishKebabTrd1Module *mod=0, *mTmp; // current module
+    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")) { //27-jan-05 - as for TRD1 as for TRD2
+        mod = (AliEMCALShishKebabTrd1Module*)fShishKebabModules->At(iz);
+        angle = mod->GetThetaInDegree();
+        if(!gn.Contains("WSUC")) { // ALICE 
+          if(iz==0) AliMatrix(idrotm, 0.,0., 90.,90., 90.,0.); // z'(x); y'(y); x'(z)
+          else      AliMatrix(idrotm, 90.-angle,180., 90.0,90.0, angle, 0.);
+          phiOK = mod->GetCenterOfModule().Phi()*180./TMath::Pi(); 
+          printf(" %2i | angle | %6.3f - %6.3f = %6.3f(eta %5.3f)\n", 
+          iz+1, angle, phiOK, angle-phiOK, mod->GetEtaOfCenterOfModule());
+          xpos = mod->GetPosXfromR() + g->GetSteelFrontThickness() - smodPar0;
+          zpos = mod->GetPosZ() - smodPar2;
+          for(int iy=0; iy<g->GetNPhi(); iy++) { // flat in phi
+            ypos = g->GetPhiModuleSize()*(2*iy+1 - g->GetNPhi())/2.;
+            gMC->Gspos(child, ++nr, mother, xpos, ypos, zpos, idrotm, "ONLY") ;
+        //printf(" %2i xpos %7.2f ypos %7.2f zpos %7.2f idrotm %i\n", nr, xpos, ypos, zpos, idrotm);
+          }
+       } else {
+          if(iz==0) AliMatrix(idrotm, 0.,0., 90.,0., 90.,90.); // (x')z; y'(x); z'(y)
+          else      AliMatrix(idrotm, 90-angle,270., 90.0,0.0, angle,90.);
+          phiOK = mod->GetCenterOfModule().Phi()*180./TMath::Pi(); 
+          printf(" %2i | angle -phiOK | %6.3f - %6.3f = %6.3f(eta %5.3f)\n", 
+          iz+1, angle, phiOK, angle-phiOK, mod->GetEtaOfCenterOfModule());
+          zpos = mod->GetPosZ()      - smodPar2;
+          ypos = mod->GetPosXfromR() - smodPar1;
+          printf(" zpos %7.2f ypos %7.2f idrotm %i\n xpos ", zpos, xpos, idrotm);
+          for(int ix=0; ix<g->GetNPhi(); ix++) { // flat in phi
+            xpos = g->GetPhiModuleSize()*(2*ix+1 - g->GetNPhi())/2.;
+            gMC->Gspos(child, ++nr, mother, xpos, ypos, zpos, idrotm, "ONLY") ;
+            printf(" %7.2f ", xpos);
+          }
+          printf("\n");
+        }
+      } else if(gn.Contains("TRD2")){ // 1-feb-05 - TRD2;  curve in phi
+        double angEtaRow = 0.;
+       double theta1=0.,phi1=0., theta2=0.,phi2=0., theta3=0.,phi3=0.;
+        angle=90.;
+        if(iz==0) {
+          mod   = new AliEMCALShishKebabTrd1Module();
+        } else {
+          mTmp  = new AliEMCALShishKebabTrd1Module(*mod);
+          mod   = mTmp;
+          angle = mod->GetThetaInDegree();
+        }
+
+        fShishKebabModules->Add(mod);
+        phiOK = mod->GetCenterOfModule().Phi()*180./TMath::Pi(); 
+        printf(" %i | theta | %6.3f - %6.3f = %6.3f\n", iz+1, angle, phiOK, angle-phiOK);
+
+        zpos = mod->GetPosZ() - parTubs[2];
+        rpos = parTubs[0] + mod->GetPosXfromR();
+
+        angle     = mod->GetThetaInDegree();
+        Double_t stepAngle = (parTubs[4] -  parTubs[3])/g->GetNPhi(); // 11-mar-04
+       for(int iy=0; iy<g->GetNPhi(); iy++) {
+          angEtaRow = parTubs[3] + stepAngle*(0.5+double(iy));
+         //          angEtaRow = 0;
+          theta1  = 90. +  angle; phi1 = angEtaRow;      // x' ~-z;
+          theta2  = 90.;          phi2 = 90. + angEtaRow;// y' ~ y;
+          theta3  = angle;        phi3 = angEtaRow;      // z' ~ x;
+          if(phi3 < 0.0) phi3 += 360.; 
+          AliMatrix(idrotm, theta1,phi1, theta2,phi2, theta3,phi3);
+
+          xpos = rpos * TMath::Cos(angEtaRow*TMath::DegToRad());
+          ypos = rpos * TMath::Sin(angEtaRow*TMath::DegToRad());
+          gMC->Gspos(child, ++nr, "SMOP", xpos, ypos, zpos, idrotm, "ONLY") ;
+         // SMON; 
+         phi1    = 180 + angEtaRow;
+         theta3  = 180.-theta3;  phi3 = angEtaRow;
+          AliMatrix(idrotm, theta1,phi1, theta2,phi2, theta3,phi3);
+          gMC->Gspos(child,  nr, "SMON", xpos, ypos, -zpos, idrotm, "ONLY") ;
+          if(0) {
+           printf(" angEtaRow(phi) %7.2f |  angle(eta) %7.2f \n",  angEtaRow, angle);
+           printf("iy=%2i xpos %7.2f ypos %7.2f zpos %7.2f idrotm %i\n", iy, xpos, ypos, zpos, idrotm);
+          }
+        } // for(int iy=0; iy<g->GetNPhi(); iy++)
+      }
+    } 
+  } else {
+    xpos = g->GetSteelFrontThickness()/2.;
+    for(int iz=0; iz<g->GetNZ(); iz++) {
+      zpos = -smodPar2 + g->GetEtaModuleSize()*(2*iz+1)/2.;
+      for(int iy=0; iy<g->GetNPhi(); iy++) {
+        ypos = g->GetPhiModuleSize()*(2*iy+1 - g->GetNPhi())/2.;
+        gMC->Gspos(child, ++nr, mother, xpos, ypos, zpos, 0, "ONLY") ;
+      //printf(" %2i xpos %7.2f ypos %7.2f zpos %7.2f \n", nr, xpos, ypos, zpos);
+      }
+    }
+  }
+  printf(" Number of modules in Super Module %i \n", nr);
+}
+
+// 8-dec-04 by PAI
+void AliEMCALv0::Trd1Tower3X3(const double parSCM0[4])
+{// PB should be for whole SCM0 - ?
+  double parTRAP[11], *dummy=0;
+  AliEMCALGeometry * g = GetGeometry(); 
+  TString gn(g->GetName()), scmx; 
+  gn.ToUpper(); 
+ // Division to tile size 
+  Info("Trd1Tower3X3()","<I> Divide SCM0 on y-axis %i", g->GetNETAdiv());
+  gMC->Gsdvn("SCMY","SCM0", g->GetNETAdiv(), 2); // y-axis
+  double dx1=parSCM0[0], dx2=parSCM0[1], dy=parSCM0[2], dz=parSCM0[3];
+  double ndiv=3., xpos=0.0;
+  // should be defined once
+  gMC->Gsvolu("PBTI", "BOX", idtmed[idPB], dummy, 0);
+  if(gn.Contains("TEST")==0) { // one name for all trapesoid
+    scmx = "SCMX"; 
+    gMC->Gsvolu(scmx.Data(), "TRAP", idtmed[idSC], dummy, 0);
+  }
+
+  
+  for(int ix=1; ix<=3; ix++) { // 3X3
+    // ix=1
+    parTRAP[0] = dz;
+    parTRAP[1] = TMath::ATan2((dx2-dx1)/2.,2.*dz)*TMath::RadToDeg(); // theta
+    parTRAP[2] = 0.;           // phi
+    // bottom
+    parTRAP[3] = dy/ndiv;      // H1
+    parTRAP[4] = dx1/ndiv;     // BL1
+    parTRAP[5] = parTRAP[4];   // TL1
+    parTRAP[6] = 0.0;          // ALP1
+    // top
+    parTRAP[7] = dy/ndiv;      // H2
+    parTRAP[8] = dx2/ndiv;     // BL2
+    parTRAP[9] = parTRAP[8];   // TL2
+    parTRAP[10]= 0.0;          // ALP2
+    xpos = +(dx1+dx2)/3.;      // 6 or 3
+
+    if      (ix==3) {
+      parTRAP[1] = -parTRAP[1];
+      xpos = -xpos;
+    } else if(ix==2) { // central part is box but we treat as trapesoid due to numbering
+      parTRAP[1] = 0.;
+      parTRAP[8] = dx1/ndiv;     // BL2
+      parTRAP[9] = parTRAP[8];   // TL2
+      xpos = 0.0;
+    }
+    printf(" ** TRAP ** xpos %9.3f\n", xpos);
+    for(int i=0; i<11; i++) printf(" par[%2.2i] %9.4f\n", i, parTRAP[i]);
+    if(gn.Contains("TEST")){
+      scmx = "SCX"; scmx += ix;
+      gMC->Gsvolu(scmx.Data(), "TRAP", idtmed[idSC], 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) ;
+    }
+    PbInTrap(parTRAP, scmx);
+  }
+
+  Info("Trd1Tower3X3()", "Ver. 1.0 : was tested.");
+}
+
+// 8-dec-04 by PAI
+void AliEMCALv0::PbInTrap(const double parTRAP[11], TString n)
+{// see for example CreateShishKebabGeometry(); just for case TRD1
+  static int nr=0;
+  printf(" Pb tiles : nrstart %i\n", nr);
+  AliEMCALGeometry * g = GetGeometry(); 
+
+  double par[3];
+  double sampleWidth = double(g->GetECPbRadThick()+g->GetECScintThick());
+  double xpos = 0.0, ypos = 0.0;
+  double zpos = -sampleWidth*g->GetNECLayers()/2. + g->GetECPbRadThick()/2.;
+
+  double coef = (parTRAP[8] -  parTRAP[4]) / (2.*parTRAP[0]);
+  double xCenterSCMX =  (parTRAP[4] +  parTRAP[8])/2.; // ??
+  //  double tan = TMath::Tan(parTRAP[1]*TMath::DegToRad());
+
+  par[1] = parTRAP[3];              // y 
+  par[2] = g->GetECPbRadThick()/2.; // z
+  for(int iz=0; iz<g->GetNECLayers(); iz++){
+    par[0] = parTRAP[4] + coef*sampleWidth*iz;
+    xpos   = par[0] - xCenterSCMX;
+    if(parTRAP[1] < 0.) xpos = -xpos;
+    gMC->Gsposp("PBTI", ++nr, n.Data(), xpos, ypos, zpos, 0, "ONLY", par, 3) ;
+    printf(" %i xpos %9.3f zpos %9.3f par[0] %9.3f |", iz+1, xpos, zpos, par[0]);
+    zpos += sampleWidth;
+    if(iz%2>0) printf("\n");
+  } 
+  printf(" Number of Pb tiles in SCMX %i coef %9.7f \n", nr, coef);
+  printf(" par[1] %9.3f  par[2] %9.3f ypos %9.3f \n", par[1], par[2], ypos); 
+  Info("PbInTrap", "Ver. 1.0 : was tested.");
+}
+
+// 8-dec-04 by PAI
+void AliEMCALv0::Trd1Tower4X4()
+{// see for example CreateShishKebabGeometry()
+}
+// 3-feb-05
+void AliEMCALv0::Scm0InTrd2(const AliEMCALGeometry * g, const Double_t parEMOD[5], Double_t parSCM0[5])
+{
+  // Passive material inside the detector
+  double wallThickness = g->GetPhiModuleSize()/2. -  g->GetPhiTileSize(); //Need check
+  printf(" wall thickness %7.5f \n", wallThickness);
+  for(int i=0; i<4; i++) { // on pictures sometimes I can not see 0 -> be carefull!!
+    parSCM0[i] = parEMOD[i] - wallThickness;
+    printf(" %i parSCMO %7.3f parEMOD %7.3f : dif %7.3f \n", i, parSCM0[i],parEMOD[i], parSCM0[i]-parEMOD[i]);
+  }
+  parSCM0[4] = parEMOD[4];
+  gMC->Gsvolu("SCM0", "TRD2", idtmed[idSC], parSCM0, 5); // idAIR -> idSC
+  gMC->Gspos("SCM0", 1, "EMOD", 0., 0., 0., 0, "ONLY") ;
+  // Division 
+  if(g->GetNPHIdiv()==2 && g->GetNETAdiv()==2) {
+    Division2X2InScm0(g, parSCM0);
+  } else {
+    Info("Scm0InTrd2"," no division SCM0 in this geometry |%s|\n", g->GetName());
+    assert(0);
+  }
+}
+
+void AliEMCALv0::Division2X2InScm0(const AliEMCALGeometry * g, const Double_t parSCM0[5])
+{
+  Double_t parTRAP[11], xpos=0.,ypos=0., dx1=0.0,dx2=0.,dy1=0.0,dy2=0.,dz=0.
+  ,dr1=0.0,dr2=0.;
+  idrotm=0;
+
+  Info("Division2X2InScm0","Divide SCM0 on y-axis %i\n", g->GetNETAdiv());
+  TString n("SCMX"), overLapFlagSCMY("ONLY"), overLapFlagSCMX("ONLY");
+  n = "SCM0"; // for testing - 14-mar-05
+  if(n=="SCM0"){
+    PbInTrapForTrd2(parSCM0, n);
+    // overLapFlagSCMY=overLapFlagSCMX="MANY"; // do not work
+    return;
+  }
+
+  dy1 = parSCM0[2] , dy2 = parSCM0[3], dz = parSCM0[4];
+
+  parTRAP[0] = parSCM0[4];    // dz
+  parTRAP[1] = TMath::ATan2((dy2-dy1)/2.,2.*dz)*TMath::RadToDeg();
+  parTRAP[2] = 90.;           // phi
+  // bottom
+  parTRAP[3] = parSCM0[2]/2.; // H1
+  parTRAP[4] = parSCM0[0];    // BL1
+  parTRAP[5] = parTRAP[4];    // TL1
+  parTRAP[6] = 0.0;           // ALP1
+  // top
+  parTRAP[7] = parSCM0[3]/2.; // H2
+  parTRAP[8] = parSCM0[1];    // BL2
+  parTRAP[9] = parTRAP[8];    // TL2
+  parTRAP[10]= 0.0;           // ALP2
+  printf(" ** SCMY ** \n");
+  for(int i=0; i<11; i++) printf(" par[%2.2i] %9.4f\n", i, parTRAP[i]);
+
+  idrotm=0;
+  gMC->Gsvolu("SCMY", "TRAP", idtmed[idSC], parTRAP, 11); // idAIR -> idSC
+  ypos = +(parTRAP[3]+parTRAP[7])/2.; //
+  printf(" Y shift SCMY inside SCM0 : %7.3f : opt %s\n", ypos, overLapFlagSCMY.Data()); 
+  gMC->Gspos("SCMY", 1, "SCM0", 0.0, ypos, 0.0, idrotm,  overLapFlagSCMY.Data()) ;
+  // Rotation SCMY around z-axis on 180 degree; x'=-x; y'=-y and z=z
+  AliMatrix(idrotm, 90.0,180., 90.0, 270.0, 0.0,0.0) ;
+  // We may have problem with numeration due to rotation - 4-feb-05
+  gMC->Gspos("SCMY", 2, "SCM0", 0.0, -ypos, 0.0, idrotm,  overLapFlagSCMY.Data()); 
+
+  Info("Division2X2InScm0","Divide SCMY on x-axis %i\n", g->GetNPHIdiv());
+  dx1 = parSCM0[0]; 
+  dx2 = parSCM0[1]; 
+  dr1=TMath::Sqrt(dx1*dx1+dy1*dy1);
+  dr2=TMath::Sqrt(dx2*dx2+dy2*dy2);
+
+  parTRAP[0] = parSCM0[4];    // dz
+  parTRAP[1] = TMath::ATan2((dr2-dr1)/2.,2.*dz)*TMath::RadToDeg(); // 
+  parTRAP[2] = 45.;           // phi
+  // bottom
+  parTRAP[3] = parSCM0[2]/2.; // H1
+  parTRAP[4] = parSCM0[0]/2.; // BL1
+  parTRAP[5] = parTRAP[4];    // TL1
+  parTRAP[6] = 0.0;           // ALP1
+  // top
+  parTRAP[7] = parSCM0[3]/2.; // H2
+  parTRAP[8] = parSCM0[1]/2;  // BL2
+  parTRAP[9] = parTRAP[8];    // TL2
+  parTRAP[10]= 0.0;           // ALP2
+  printf(" ** SCMX ** \n");
+  for(int i=0; i<11; i++) printf(" par[%2.2i] %9.4f\n", i, parTRAP[i]);
+
+  idrotm=0;
+  gMC->Gsvolu("SCMX", "TRAP", idtmed[idSC], parTRAP, 11);
+  xpos = (parTRAP[4]+parTRAP[8])/2.;
+  printf(" X shift SCMX inside SCMX : %7.3f : opt %s\n", xpos, overLapFlagSCMX.Data()); 
+  gMC->Gspos("SCMX", 1, "SCMY", xpos, 0.0, 0.0, idrotm,  overLapFlagSCMX.Data()) ;
+  //  AliMatrix(idrotm, 90.0,270., 90.0, 0.0, 0.0,0.0); // x'=-y; y'=x; z'=z
+  AliMatrix(idrotm, 90.0,90., 90.0, -180.0, 0.0,0.0);     // x'=y;  y'=-x; z'=z
+  gMC->Gspos("SCMX", 2, "SCMY", -xpos, 0.0, 0.0, idrotm,  overLapFlagSCMX.Data()) ;
+  // PB:
+  if(n=="SCMX" && overLapFlagSCMY == "ONLY") {
+    PbInTrapForTrd2(parTRAP, n);
+  }
+} 
+
+// 4-feb-05 by PAI
+void AliEMCALv0::PbInTrapForTrd2(const double *parTRAP, TString name)
+{// TRD2 cases
+  Double_t *dummy=0;
+  TString pbShape("BOX"), pbtiChonly("ONLY");
+  if(name=="SCM0") {
+    pbShape    = "TRD2";
+    //    pbtiChonly = "MANY";
+  }
+  gMC->Gsvolu("PBTI", pbShape.Data(), idtmed[idPB], dummy, 0);
+
+  int nr=0;
+  Info("PbInTrapForTrd2"," Pb tiles inside %s: shape %s :pbtiChonly %s\n nrstart %i\n", 
+  name.Data(), pbShape.Data(), pbtiChonly.Data(), nr);
+  AliEMCALGeometry * g = GetGeometry(); 
+
+  double par[5], parPB[5], parSC[5];
+  double sampleWidth = double(g->GetECPbRadThick()+g->GetECScintThick());
+  double xpos = 0.0, ypos = 0.0;
+  double zpos = -sampleWidth*g->GetNECLayers()/2. + g->GetECPbRadThick()/2.;
+  if(name == "SCMX") { // common trapezoid - 11 parameters
+    double coef = (parTRAP[8] -  parTRAP[4]) / (2.*parTRAP[0]);
+    double xCenterSCMX =  (parTRAP[4] +  parTRAP[8])/2.; // the same for y
+    printf(" xCenterSCMX %8.5f : coef %8.7f \n", xCenterSCMX, coef);
+
+    par[2] = g->GetECPbRadThick()/2.; // z
+    for(int iz=0; iz<g->GetNECLayers(); iz++){
+      par[0] = parTRAP[4] + coef*sampleWidth*iz;
+      par[1] = par[0];
+      xpos   = ypos = par[0] - xCenterSCMX;
+    //if(parTRAP[1] < 0.) xpos = -xpos;
+      gMC->Gsposp("PBTI", ++nr, name.Data(), xpos, ypos, zpos, 0, "ONLY", par, 3) ;
+      printf(" %2.2i xpos %8.5f zpos %6.3f par[0,1] %6.3f |", iz+1, xpos, zpos, par[0]);
+      if(iz%2>0) printf("\n");
+      zpos += sampleWidth;
+    } 
+    printf(" Number of Pb tiles in SCMX %i coef %9.7f \n", nr, coef);
+    printf(" par[1] %9.5f  par[2] %9.5f ypos %9.5f \n", par[1], par[2], ypos); 
+  } else if(name == "SCM0") { // 1-mar-05 ; TRD2 - 5 parameters
+    printf(" SCM0 par = ");
+    for(int i=0; i<5; i++) printf(" %9.5f ", parTRAP[i]);
+    printf("\n zpos %f \n",zpos);
+
+    double tanx = (parTRAP[1] -  parTRAP[0]) / (2.*parTRAP[4]); //  tanx =  tany now
+    double tany = (parTRAP[3] -  parTRAP[2]) / (2.*parTRAP[4]), ztmp=0.;
+    parPB[4] = g->GetECPbRadThick()/2.;
+    parSC[2] = g->GetECScintThick()/2.;
+    for(int iz=0; iz<g->GetNECLayers(); iz++){
+      ztmp     = sampleWidth*double(iz);
+      parPB[0] = parTRAP[0] + tanx*ztmp;
+      parPB[1] = parPB[0]   + tanx*g->GetECPbRadThick();
+      parPB[2] = parTRAP[2] + tany*ztmp;
+      parPB[3] = parPB[2]   + tany*g->GetECPbRadThick();
+      gMC->Gsposp("PBTI", ++nr, name.Data(), xpos, ypos, zpos, 0, pbtiChonly.Data(), parPB, 5) ;
+      printf("\n PBTI %2i | zpos %6.3f | par = ", nr, zpos);
+      /*
+      for(int i=0; i<5; i++) printf(" %9.5f ", parPB[i]);
+      // individual SC tile
+      parSC[0] = parPB[0];
+      parSC[1] = parPB[1];
+      gMC->Gsposp("SCTI", nr, name.Data(), xpos, ypos, zpos+g->GetECScintThick(), 
+      0, pbtiChonly.Data(), parSC, 3) ;
+      printf("\n SCTI     zpos %6.3f | par = ", zpos+g->GetECScintThick());
+      for(int i=0; i<3; i++) printf(" %9.5f ", parPB[i]);
+      */
+      zpos  += sampleWidth;
+    }
+    printf("\n");
+  }
+  Info("PbInTrapForTrd2", "Ver. 0.03 : was tested.");
+}
+
+// 15-mar-05
+void AliEMCALv0::PbmoInTrd2(const AliEMCALGeometry * g, const Double_t parEMOD[5], Double_t parPBMO[5])
+{
+  Info("PbmoInTrd2"," started : geometry %s ", g->GetName());
+  double wallThickness = g->GetPhiModuleSize()/2. -  g->GetPhiTileSize();
+  printf(" wall thickness %7.5f \n", wallThickness);
+  for(int i=0; i<4; i++) {
+    parPBMO[i] = parEMOD[i] - wallThickness;
+    printf(" %i parPBMO %7.3f parEMOD %7.3f : dif %7.3f \n", 
+    i, parPBMO[i],parEMOD[i], parPBMO[i]-parEMOD[i]);
+  }
+  parPBMO[4] = parEMOD[4];
+  gMC->Gsvolu("PBMO", "TRD2", idtmed[idPB], parPBMO, 5);
+  gMC->Gspos("PBMO", 1, "EMOD", 0., 0., 0., 0, "ONLY") ;
+  // Division 
+  if(g->GetNPHIdiv()==2 && g->GetNETAdiv()==2) {
+    Division2X2InPbmo(g, parPBMO);
+    printf(" PBMO, division 2X2 | geometry |%s|\n", g->GetName());
+  } else {
+    printf(" no division PBMO in this geometry |%s|\n", g->GetName());
+    assert(0);
+  }
+}
+
+void AliEMCALv0::Division2X2InPbmo(const AliEMCALGeometry * g, const Double_t parPBMO[5])
+{
+  Info("Division2X2InPbmo"," started : geometry %s ", g->GetName());
+  //Double_t *dummy=0;
+  //  gMC->Gsvolu("SCTI", "BOX", idtmed[idSC], dummy, 0);
+
+  double parSC[3];
+  double sampleWidth = double(g->GetECPbRadThick()+g->GetECScintThick());
+  double xpos = 0.0, ypos = 0.0, zpos = 0.0, ztmp=0;;
+  double tanx = (parPBMO[1] -  parPBMO[0]) / (2.*parPBMO[4]); //  tanx =  tany now
+  double tany = (parPBMO[3] -  parPBMO[2]) / (2.*parPBMO[4]);
+  char name[10], named[10], named2[10];
+
+  printf(" PBMO par = ");
+  for(int i=0; i<5; i++) printf(" %9.5f ", parPBMO[i]);
+  printf("\n");
+
+  parSC[2] = g->GetECScintThick()/2.;
+  zpos = -sampleWidth*g->GetNECLayers()/2. + g->GetECPbRadThick() + g->GetECScintThick()/2.;
+  printf(" parSC[2] %9.5f \n", parSC[2]);
+  for(int iz=0; iz<g->GetNECLayers(); iz++){
+    ztmp     = g->GetECPbRadThick() + sampleWidth*double(iz); // Z for previous PB
+    parSC[0] =  parPBMO[0] + tanx*ztmp;
+    parSC[1] =  parPBMO[2] + tany*ztmp;
+
+    sprintf(name,"SC%2.2i", iz+1);
+    gMC->Gsvolu(name, "BOX", idtmed[idSC], parSC, 3);
+    gMC->Gspos(name, 1, "PBMO", xpos, ypos, zpos, 0, "ONLY") ;
+    printf("%s | zpos %6.3f | parSC[0,1]=(%7.5f,%7.5f) -> ", 
+    name, zpos, parSC[0], parSC[1]);
+    
+    sprintf(named,"SY%2.2i", iz+1);
+    printf(" %s -> ", named);
+    gMC->Gsdvn(named,name, 2, 2);
+
+    sprintf(named2,"SX%2.2i", iz+1);
+    printf(" %s \n", named2);
+    gMC->Gsdvn(named2,named, 2, 1);
+
+    zpos    += sampleWidth;
+  }
+}
index a722709..7be4820 100644 (file)
@@ -15,6 +15,8 @@
 // --- ROOT system ---
 
 class TFile;
+class TList;
+class TNode;
 
 // --- AliRoot header files ---
 #include "AliEMCAL.h"
@@ -35,6 +37,7 @@ class AliEMCALv0 : public AliEMCAL {
   virtual ~AliEMCALv0(){} 
 
   virtual void BuildGeometry();// creates the geometry for the ROOT display
+  TNode *BuildGeometryOfWSUC();  // WSUC - test environment
   virtual void CreateGeometry() ;// creates the geometry for GEANT
   virtual void   Init(void) ;                                       // does nothing
   virtual Int_t  IsVersion(void) const { 
@@ -51,6 +54,23 @@ class AliEMCALv0 : public AliEMCAL {
     Fatal("operator =", "not implemented") ;  
     return *this ; 
   }
+  // ShishKebab 
+  void CreateShishKebabGeometry();
+  void CreateSmod(const char* mother="XEN1");
+  void CreateEmod(const char* mother="SMOD", const char* child="EMOD");
+  // TRD1
+  void Trd1Tower3X3(const double parSCM0[5]);
+  void Trd1Tower4X4();
+  void PbInTrap(const double parTRAP[11], TString n);
+  // TRD2 - 1th design
+  void Scm0InTrd2(const AliEMCALGeometry * g, const Double_t parEMOD[5], Double_t parSCM0[5]);
+  void Division2X2InScm0(const AliEMCALGeometry * g, const Double_t parSCM0[5]);
+  void PbInTrapForTrd2(const double *parTRAP, TString name);
+  // TRD2 - 2th design
+  void PbmoInTrd2(const AliEMCALGeometry * g, const Double_t parEMOD[5], Double_t parPBMO[5]);
+  void Division2X2InPbmo(const AliEMCALGeometry * g, const Double_t parPBMO[5]);
+
+  TList *fShishKebabModules; //! list of modules for twist geometries
   
  protected:
 
index 48e3ca0..be1a95b 100644 (file)
@@ -57,7 +57,7 @@ AliEMCALv1::AliEMCALv1(const char *name, const char *title):
     // Standard Creator.
 
     fHits= new TClonesArray("AliEMCALHit",1000);
-    gAlice->GetMCApp()->AddHitList(fHits);
+    //    gAlice->GetMCApp()->AddHitList(fHits); // 20-dec-04 - advice of Andreas
 
     fNhits = 0;
     fIshunt     =  2; // All hits are associated with particles entering the calorimeter
@@ -104,7 +104,7 @@ void AliEMCALv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t ipa
        new((*fHits)[fNhits]) AliEMCALHit(*newHit);
        fNhits++;
     } // end if
-
+    
     delete newHit;
 }
 //______________________________________________________________________
index 18d28f1..15bcab4 100644 (file)
@@ -47,13 +47,14 @@ public:
     return *this;}
  
     
-private:
+protected:
+  // Marco advice - 16-jan-05
   Int_t fCurPrimary;
   Int_t fCurParent;
   Int_t fCurTrack;
   Float_t fTimeCut;       // Cut to remove the background from the ALICE system
 
-  ClassDef(AliEMCALv1,8)//Implementation of EMCAL manager class to produce hits in a Central Calorimeter 
+  ClassDef(AliEMCALv1,9)//Implementation of EMCAL manager class to produce hits in a Central Calorimeter 
     
 };
 
diff --git a/EMCAL/AliEMCALv2.cxx b/EMCAL/AliEMCALv2.cxx
new file mode 100644 (file)
index 0000000..7234e0c
--- /dev/null
@@ -0,0 +1,569 @@
+/**************************************************************************
+ * Copyright(c) 1998-2004, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+//_________________________________________________________________________
+//*-- Implementation version v2 of EMCAL Manager class 
+//*-- An object of this class does not produce digits
+//*-- It is the one to use if you do want to produce outputs in TREEH 
+//*--                  
+//*-- Author: Sahal Yacoob (LBL /UCT)
+//*--       : Jennifer Klay (LBL)
+// This Class not stores information on all particles prior to EMCAL entry - in order to facilitate analysis.
+// This is done by setting fIShunt =2, and flagging all parents of particles entering the EMCAL.
+
+// 15/02/2002 .... Yves Schutz
+//  1. fSamplingFraction and fLayerToPreshowerRatio have been removed
+//  2. Timing signal is collected and added to hit
+
+// --- ROOT system ---
+#include "TParticle.h"
+#include "TVirtualMC.h"
+#include "TBrowser.h"
+#include "TH1.h"
+#include "TH2.h"
+
+// --- Standard library ---
+
+// --- AliRoot header files ---
+#include "AliEMCALv2.h"
+#include "AliEMCALHit.h"
+#include "AliEMCALGeometry.h"
+#include "AliRun.h"
+#include "AliHeader.h"
+#include "AliMC.h"
+#include "AliPoints.h"
+#include "AliEMCALGetter.h"
+// for TRD1,2 case
+
+ClassImp(AliEMCALv2)
+
+  //extern "C" void gdaxis_(float &x0, float &y0, float &z0, float &axsiz);
+
+
+//______________________________________________________________________
+AliEMCALv2::AliEMCALv2():AliEMCALv1(), fGeometry(0){
+  // ctor
+
+}
+
+//______________________________________________________________________
+AliEMCALv2::AliEMCALv2(const char *name, const char *title): AliEMCALv1(name,title) {
+    // Standard Creator.
+
+  //    fGeant3 = (TGeant3*)gMC;
+    fHits= new TClonesArray("AliEMCALHit",1000);
+    gAlice->GetMCApp()->AddHitList(fHits);
+
+    fNhits    = 0;
+    fIshunt   = 2; // All hits are associated with particles entering the calorimeter
+    fTimeCut  = 30e-09;
+
+    fGeometry = GetGeometry(); 
+    fHDe = 0;
+    //    if (gDebug>0){
+    if (1){
+      TH1::AddDirectory(0);
+      fHDe = new TH1F("fHDe","De in EMCAL", 1000, 0., 1.);
+      fHistograms->Add(fHDe);
+      TH1::AddDirectory(1);
+    }
+}
+
+//______________________________________________________________________
+AliEMCALv2::~AliEMCALv2(){
+    // dtor
+
+    if ( fHits) {
+       fHits->Delete();
+       delete fHits;
+       fHits = 0;
+    }
+}
+
+//______________________________________________________________________
+void AliEMCALv2::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t iparent, Float_t ienergy, 
+                       Int_t id, Float_t * hits,Float_t * p){
+    // Add a hit to the hit list.
+    // An EMCAL hit is the sum of all hits in a tower section
+    //   originating from the same entering particle 
+    static Int_t hitCounter;
+    static AliEMCALHit *newHit, *curHit;
+    static Bool_t deja;
+
+    deja = kFALSE;
+
+    newHit = new AliEMCALHit(shunt, primary, tracknumber, iparent, ienergy, id, hits, p);
+    for ( hitCounter = fNhits-1; hitCounter >= 0 && !deja; hitCounter-- ) {
+       curHit = (AliEMCALHit*) (*fHits)[hitCounter];
+       // We add hits with the same tracknumber, while GEANT treats
+       // primaries succesively
+       if(curHit->GetPrimary() != primary) 
+         break;
+       if( *curHit == *newHit ) {
+           *curHit = *curHit + *newHit;
+           deja = kTRUE;
+           //            break; // 30-aug-04 by PAI 
+       } // end if
+    } // end for hitCounter
+    
+    if ( !deja ) {
+       new((*fHits)[fNhits]) AliEMCALHit(*newHit);
+       fNhits++;
+    }
+    //    printf(" fNhits %i \n", fNhits); 
+    delete newHit;
+}
+//______________________________________________________________________
+void AliEMCALv2::StepManager(void){
+  // Accumulates hits as long as the track stays in a tower
+
+  // position wrt MRS and energy deposited
+  static Float_t        xyzte[5]={0.,0.,0.,0.,0.};// position wrt MRS, time and energy deposited
+  static Float_t        pmom[4]={0.,0.,0.,0.};
+  static TLorentzVector pos; // Lorentz vector of the track current position.
+  static TLorentzVector mom; // Lorentz vector of the track current momentum.
+  static Float_t ienergy = 0;
+  static TString curVolName;
+  static int supModuleNumber, moduleNumber, yNumber, xNumber, absid;
+  static int keyGeom=0;
+  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; 
+
+  if(keyGeom == 0) {
+    keyGeom = 2;
+    if(gMC->VolId("PBMO")==0 || gMC->VolId("WSUC")==1) {
+      vn      = "SCMX";   // old TRD2(TRD1) or WSUC
+      keyGeom = 1;
+    }    
+    printf("AliEMCALv2::StepManager():  keyGeom %i : Sensetive volume %s \n", 
+    keyGeom, vn); 
+    if(gMC->VolId("WSUC")==1) printf(" WSUC - cosmic ray stand geometry \n");
+  }
+  Int_t tracknumber =  gAlice->GetMCApp()->GetCurrentTrackNumber();
+
+  curVolName = gMC->CurrentVolName();
+  if(curVolName.Contains(vn)) { // We are in a scintillator layer
+    //    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
+      //       Info("StepManager "," entry %i DE %f",++ientry, depositedEnergy); // for testing
+       if (fCurPrimary==-1) 
+       fCurPrimary=gAlice->GetMCApp()->GetPrimary(tracknumber);
+
+      if (fCurParent==-1 || tracknumber != fCurTrack) {
+       // Check parentage
+       Int_t parent=tracknumber;
+       if (fCurParent != -1) {
+         while (parent != fCurParent && parent != -1) {
+           TParticle *part=gAlice->GetMCApp()->Particle(parent);
+           parent=part->GetFirstMother();
+         }
+       }
+       if (fCurParent==-1 || parent==-1) {
+         Int_t parent=tracknumber;
+         TParticle *part=gAlice->GetMCApp()->Particle(parent);
+         while (parent != -1 && fGeometry->IsInEMCAL(part->Vx(),part->Vy(),part->Vz())) {
+           parent=part->GetFirstMother();
+           if (parent!=-1) 
+             part=gAlice->GetMCApp()->Particle(parent);
+         } 
+         fCurParent=parent;
+         if (fCurParent==-1)
+           Error("StepManager","Cannot find parent");
+         else {
+           TParticle *part=gAlice->GetMCApp()->Particle(fCurParent);
+           ienergy = part->Energy(); 
+         }
+         while (parent != -1) {
+           part=gAlice->GetMCApp()->Particle(parent);
+           part->SetBit(kKeepBit);
+           parent=part->GetFirstMother();
+         }
+       }
+       fCurTrack=tracknumber;
+      }    
+      gMC->TrackPosition(pos);
+      xyzte[0] = pos[0];
+      xyzte[1] = pos[1];
+      xyzte[2] = pos[2];
+      xyzte[3] = gMC->TrackTime() ;       
+      
+      gMC->TrackMomentum(mom);
+      pmom[0] = mom[0];
+      pmom[1] = mom[1];
+      pmom[2] = mom[2];
+      pmom[3] = mom[3];
+      
+      //      if(ientry%200 > 0) return; // testing
+      supModuleNumber = moduleNumber = yNumber = xNumber = absid = 0;
+      if(keyGeom >= 1) { // old style
+        gMC->CurrentVolOffID(4, supModuleNumber);
+        gMC->CurrentVolOffID(3, moduleNumber);
+        gMC->CurrentVolOffID(1, yNumber);
+        gMC->CurrentVolOffID(0, xNumber); // really x number now
+      } else {
+        gMC->CurrentVolOffID(5, supModuleNumber);
+        gMC->CurrentVolOffID(4, moduleNumber);
+        gMC->CurrentVolOffID(1, yNumber);
+        gMC->CurrentVolOffID(0, xNumber); // really x number now
+        if    (strcmp(gMC->CurrentVolOffName(5),"SMOP")==0) supModuleNumber = nSMOP[supModuleNumber-1];
+        else if(strcmp(gMC->CurrentVolOffName(5),"SMON")==0) supModuleNumber = nSMON[supModuleNumber-1];
+        else   assert(0); // something wrong
+      }
+      absid = fGeometry->GetAbsCellId(supModuleNumber, moduleNumber, yNumber, xNumber);
+    
+      if (absid == -1) Fatal("StepManager()", "Wrong id ") ;
+
+      Float_t lightYield =  depositedEnergy ;
+      // Apply Birk's law (copied from G3BIRK)
+
+      if (gMC->TrackCharge()!=0) { // Check
+         Float_t BirkC1_mod = 0;
+       if (fBirkC0==1){ // Apply correction for higher charge states
+         if (TMath::Abs(gMC->TrackCharge())>=2) BirkC1_mod=fBirkC1*7.2/12.6;
+         else                                    BirkC1_mod=fBirkC1;
+       }
+
+       Float_t dedxcm;
+       if (gMC->TrackStep()>0)  dedxcm=1000.*gMC->Edep()/gMC->TrackStep();
+       else                     dedxcm=0;
+       lightYield=lightYield/(1.+BirkC1_mod*dedxcm+fBirkC2*dedxcm*dedxcm);
+      } 
+
+      // use sampling fraction to get original energy --HG
+      //      xyzte[4] = lightYield * fGeometry->GetSampling();
+      xyzte[4] = lightYield; // 15-dec-04
+        
+      if (gDebug == -2) 
+      printf("#sm %2i #m %3i #x %1i #z %1i -> absid %i : xyzte[4] = %f\n",
+      supModuleNumber,moduleNumber,yNumber,xNumber,absid, xyzte[4]);
+
+      AddHit(fIshunt, fCurPrimary,tracknumber, fCurParent, ienergy, absid,  xyzte, pmom);
+    } // there is deposited energy
+  }
+}
+
+void AliEMCALv2::FinishEvent()
+{ // 26-may-05
+  static double de=0.;
+  de = GetDepositEnergy(0);
+  if(fHDe) fHDe->Fill(de);
+}
+
+Double_t AliEMCALv2::GetDepositEnergy(int print)
+{ // 23-mar-05 - for testing
+  if(fHits == 0) return 0.;
+  AliEMCALHit  *hit=0;
+  Double_t de=0.;
+  for(int ih=0; ih<fHits->GetEntries(); ih++) {
+    hit = (AliEMCALHit*)fHits->UncheckedAt(ih);
+    de += hit->GetEnergy();
+  }
+  if(print>0) {
+    printf(" #hits %i de %f \n", fHits->GetEntries(), de);
+    if(print>1) {
+      printf(" #primary particles %i\n", gAlice->GetHeader()->GetNprimary()); 
+    }
+  }
+  return de;
+}
+
+void AliEMCALv2::Browse(TBrowser* b)
+{
+  TObject::Browse(b);
+}
+
+void AliEMCALv2::DrawCalorimeterCut(const char *name, int axis, double dcut)
+{ // Size of tower is 5.6x5.6x24.8 (25.0); cut on Z axiz
+  gMC->Gsatt("*", "seen", 0);
+
+  int fill = 1;
+
+  if(axis!=1) SetVolumeAttributes("STPL", 1, 1, fill);
+
+  int colo=4;
+  TString sn(name);
+  if(sn.Contains("SCM")) colo=5;
+  SetVolumeAttributes(name, 1, colo, fill);
+
+  TString st("Shish-Kebab, Compact, zcut , ");
+  st += name;
+
+  char *optShad = "on", *optHide="on";
+  double cxy=0.02;
+  if     (axis==1) {
+    dcut = 0.;
+    //    optHide = optShad = "off";
+  } else if(axis==2){
+    if(dcut==0.) SetVolumeAttributes("SCM0", 1, 5, fill); // yellow    
+  }
+
+  gMC->Gdopt("hide", optHide);
+  gMC->Gdopt("shad", optShad);
+
+  gROOT->ProcessLine("TGeant3 *g3 = (TGeant3*)gMC");
+  char cmd[200];
+  sprintf(cmd,"g3->Gdrawc(\"XEN1\", %i, %5.1f, 12., 8., %3.2f, %3.2f)\n", axis, dcut, cxy,cxy);
+  printf("\n%s\n",cmd); gROOT->ProcessLine(cmd);
+
+  sprintf(cmd,"gMC->Gdhead(1111, \"%s\")\n", st.Data());
+  printf("%s\n",cmd); gROOT->ProcessLine(cmd);
+}
+
+void AliEMCALv2::DrawSuperModuleCut(const char *name, int axis, double dcut, int fill)
+{ // Size of tower is 5.6x5.6x24.8 (25.0); cut on Z axiz
+  TString sn(GetGeometry()->GetName());
+  sn.ToUpper();
+  char *tit[3]={"xcut", "ycut", "zcut"};
+  if(axis<1) axis=1; if(axis>3) axis=3;
+
+  gMC->Gsatt("*", "seen", 0);
+  //  int fill = 6;
+
+  // SetVolumeAttributes("SMOD", 1, 1, fill);  // black
+  SetVolumeAttributes(name, 1, 5, fill);    // yellow 
+
+  double cxy=0.055, x0=10., y0=10.;
+  char *optShad = "on", *optHide="on";
+  if(sn.Contains("TRD1")) {
+    SetVolumeAttributes("STPL", 1, 3, fill);  // green 
+    if     (axis==1) {
+      gMC->Gsatt("STPL", "seen", 0);
+      dcut = 0.;
+      optHide = "off";
+      optShad = "off";
+    } else if(axis==3) cxy = 0.1;
+  } else if(sn.Contains("TRD2")) {
+    y0 = -10.;
+    if (axis==2) cxy=0.06;
+  }
+  gMC->Gdopt("hide", optHide);
+  gMC->Gdopt("shad", optShad);
+
+  TString st("Shish-Kebab, Compact, SMOD, ");
+  if(sn.Contains("TWIST")) st = "Shish-Kebab, Twist, SMOD, ";
+  st += tit[axis-1];;
+
+  gROOT->ProcessLine("TGeant3 *g3 = (TGeant3*)gMC");
+  char cmd[200];
+  sprintf(cmd,"g3->Gdrawc(\"SMOD\", %i, %6.3f, %5.1f, %5.1f, %5.3f, %5.3f)\n",axis,dcut,x0,y0,cxy,cxy);
+  printf("\n%s\n",cmd); gROOT->ProcessLine(cmd);
+
+  sprintf(cmd,"gMC->Gdhead(1111, \"%s\")\n", st.Data());
+  printf("%s\n",cmd); gROOT->ProcessLine(cmd);
+  // hint for testing
+  if(sn.Contains("TRD1")){
+    printf("Begin of super module\n");
+    printf("g3->Gdrawc(\"SMOD\", 2,  0.300, 89., 10., 0.5, 0.5)\n");
+    printf("Center of super modules\n");
+    printf("g3->Gdrawc(\"SMOD\", 2,  0.300, 0., 10., 0.5, 0.5)\n");
+    printf("end of super modules\n");
+    printf("g3->Gdrawc(\"SMOD\", 2,  0.300, -70., 10., 0.5, 0.5)\n");
+  } else if(sn.Contains("TRD2")){
+    printf("Begin of super module  ** TRD2 ** \n");
+    printf("g3->Gdrawc(\"SMOD\", 2,  0.00,  40., -80, 0.2, 0.2)\n");
+    printf("end of super modules\n");
+    printf("g3->Gdrawc(\"SMOD\", 2,  0.00, -20., -80, 0.2, 0.2)\n");
+
+    printf(" ***  Z cut (Y|X plane)\n scale 0.4\n");
+    printf("g3->Gdrawc(\"SMOD\", 3,  -170.,-165.,10.,0.4,0.4)\n");
+    printf(" scale 0.2\n");
+    printf("g3->Gdrawc(\"SMOD\", 3,  -170.,-80.,10.,0.2,0.2)\n");
+    printf(" scale 0.12\n");
+    printf("g3->Gdrawc(\"SMOD\", 3,   -170.,-45.,10.,0.12,0.12)\n");
+  }
+}
+
+void AliEMCALv2::DrawTowerCut(const char *name, int axis, double dcut, int fill, char *optShad)
+{ // Size of tower is 5.6x5.6x24.8 (25.0); cut on Z axiz
+  if(axis<1) axis=1; if(axis>3) axis=3;
+  TString mn(name); mn.ToUpper();
+  TString sn(GetGeometry()->GetName());
+
+  gMC->Gsatt("*", "seen", 0);
+  gMC->Gsatt("*", "fill", fill);
+
+  // int mainColo=5; // yellow
+  if(mn == "EMOD") {
+    SetVolumeAttributes(mn.Data(), 1, 1, fill);
+    SetVolumeAttributes("SCM0", 1, 5, fill); // yellow
+    SetVolumeAttributes("SCPA", 1, 3, fill); // green - 13-may-05
+  } else if(mn == "SCM0") { // first division 
+    SetVolumeAttributes(mn.Data(), 1, 1, fill);
+    SetVolumeAttributes("SCMY", 1, 5, fill); // yellow
+  } else if(mn == "SCMY") { // first division 
+    SetVolumeAttributes(mn.Data(), 1, 1, fill); 
+    if(sn.Contains("TEST") && sn.Contains("3X3")) {
+      SetVolumeAttributes("SCX1", 1, 5, fill); // yellow
+      SetVolumeAttributes("SCX2", 1, 2, fill); // red
+      SetVolumeAttributes("SCX3", 1, 3, fill); // green
+    } else {
+      SetVolumeAttributes("SCMX", 1, 5, fill); // yellow
+    }
+  } else if(mn == "SCMX" || mn.Contains("SCX")) {
+    SetVolumeAttributes(mn.Data(), 1, 5, fill);// yellow
+    SetVolumeAttributes("PBTI", 1, 4, fill);
+  } else {
+    printf("<W> for volume |%s| did not defined volume attributes\n", mn.Data());
+  }
+
+  //  TString st("Shish-Kebab, 2x2 mm sampling, 62 layers, ");
+  TString st("Shashlyk, 2x2 mm sampling, 62 layers, ");
+  if    (sn.Contains("25"))   st = "Shish-Kebab, 5x5 mm sampling, 25 layers, ";
+  else if(sn.Contains("MAY05"))   st = "Shish-Kebab, 5x5 mm sampling, 77 layers, ";
+  if(sn.Contains("TRD1")) st += " TRD1, ";
+  if(sn.Contains("3X3"))  st += " 3x3, ";
+  st += name;
+
+  gROOT->ProcessLine("TGeant3 *g3 = (TGeant3*)gMC");
+  double cx=0.78, cy=2.;
+  if     (axis==1 && sn.Contains("TRD")==0) {
+   cx = cy = 1.;
+   gMC->Gsatt("PBTI", "seen", 0);
+  } else if (sn.Contains("TEST") && sn.Contains("3X3")) {
+    cx = cy = 0.7;
+    if (axis==3) cx = cy = 1.;
+  } else if (sn.Contains("TRD2")) {
+    if (axis==3) cx = cy = 2.;
+  } else if (mn.Contains("EMOD")) {
+    cx = cy = 0.5;
+  }
+  char cmd[200];
+
+  gMC->Gdopt("hide", optShad);
+  gMC->Gdopt("shad", optShad);
+
+  sprintf(cmd,"g3->Gdrawc(\"%s\", %i, %6.2f, 10., 10., %5.3f, %5.3f)\n",name, axis,dcut,cx,cy);
+  printf("\n%s\n",cmd); gROOT->ProcessLine(cmd);
+
+  sprintf(cmd,"gMC->Gdhead(1111, \"%s\")\n", st.Data());
+  printf("%s\n",cmd); gROOT->ProcessLine(cmd);
+}
+  
+void AliEMCALv2::DrawAlicWithHits(int mode)
+{ // 20-sep-04; does not work now
+  static TH2F *h2;
+  if(h2==0) h2 = new TH2F("h2","test fo hits", 60,0.5,60.5, 28,0.5,28.5);
+  else      h2->Reset();
+  if(mode==0) {
+    gROOT->ProcessLine("TGeant3 *g3 = (TGeant3*)gMC");
+    gMC->Gsatt("*","seen",0);
+    gMC->Gsatt("scm0","seen",5);
+
+    gROOT->ProcessLine("g3->Gdrawc(\"ALIC\", 1,   2.0, 12., -130, 0.3, 0.3)");
+    // g3->Gdrawc("ALIC", 1,   2.0, 10., -14, 0.05, 0.05)
+  }
+
+  TClonesArray *hits = Hits();
+  Int_t nhits = hits->GetEntries(), absId, nSupMod, nTower, nIphi, nIeta, iphi, ieta;
+  AliEMCALHit *hit = 0;
+  Double_t de, des=0.;
+  for(Int_t i=0; i<nhits; i++) {
+    hit   = (AliEMCALHit*)hits->UncheckedAt(i);
+    absId = hit->GetId();
+    de    = hit->GetEnergy();
+    des += de;
+    if(fGeometry->GetCellIndex(absId, nSupMod, nTower, nIphi, nIeta)){
+      //      printf(" de %f abs id %i smod %i tower %i | cell iphi %i : ieta %i\n",
+      // de, absId, nSupMod, nTower, nIphi, nIeta); 
+      if(nSupMod==3) {
+        fGeometry->GetCellPhiEtaIndexInSModule(nTower,nIphi,nIeta, iphi,ieta);
+        // printf(" iphi %i : ieta %i\n", iphi,ieta);
+        h2->Fill(double(ieta),double(iphi), de);
+      }
+    }
+    //    hit->Dump();
+  }
+  printf(" #hits %i : sum de %f -> %f \n", nhits, des, h2->Integral());
+  if(mode>0 && h2->Integral()>0.) h2->Draw();
+}
+
+void AliEMCALv2::SetVolumeAttributes(const char *name,const int seen, const int color, const int fill)
+{
+ /* seen=-2:volume is visible but none of its descendants;
+         -1:volume is not visible together with all its descendants;
+          0:volume is not visible;
+         1:volume is     visible. */
+  gMC->Gsatt(name, "seen", seen);
+  gMC->Gsatt(name, "colo", color);
+  gMC->Gsatt(name, "fill", fill); 
+  printf(" %s : seen %i color %i fill %i \n", name, seen, color, fill);
+} 
+
+void AliEMCALv2::TestIndexTransition(int pri, int idmax)
+{ // for EMCAL_SHISH geometry
+  TString sn(fGeometry->GetName());
+  if(!sn.Contains("SHISH")) {
+    printf("Wrong geometry |%s| ! Bye \n", sn.Data());
+    return; 
+  }
+
+  Int_t nSupMod, nTower, nIphi, nIeta, idNew, nGood=0;
+  if(idmax==0) idmax = fGeometry->GetNCells();
+  for(Int_t id=1; id<=idmax; id++) {
+    if(!fGeometry->GetCellIndex(id, nSupMod, nTower, nIphi, nIeta)){
+      printf(" Wrong abs ID %i : #cells %i\n", id, fGeometry->GetNCells());
+      break;  
+    }
+    idNew = fGeometry->GetAbsCellId(nSupMod, nTower, nIphi, nIeta);
+    if(id != idNew || pri>0) {
+      printf(" ID %i : %i <- new id\n", id, idNew);
+      printf(" nSupMod   %i ",  nSupMod);
+      printf(" nTower    %i ", nTower);
+      printf(" nIphi     %i ", nIphi);
+      printf(" nIeta     %i \n", nIeta);
+     
+    } else nGood++;
+  }
+  printf(" Good decoding %i : %i <- #cells \n", nGood, fGeometry->GetNCells());
+}
+
+/* try to draw hits
+ gMC->Gsatt("*","seen",0);
+ gMC->Gsatt("scm0","seen",5);
+ g3->Gdrawc("alic", 1,   2.0, 12., -125, 0.3, 0.3);
+*/
+//void AliEMCALv2::Gdaxis(float x0, float y0, float z0, float <axsiz) {gdaxis_(x0,y0,z0,axsiz);}
+
+/*
+
+void AliEMCALv2::RemapTrackHitIDs(Int_t *map) {
+  // remap track index numbers for primary and parent indices
+  // (Called by AliStack::PurifyKine)
+  if (Hits()==0)
+    return;
+  TIter hit_it(Hits());
+  Int_t i_hit=0;
+  while (AliEMCALHit *hit=dynamic_cast<AliEMCALHit*>(hit_it()) ) {
+    if (map[hit->GetIparent()]==-99)
+      cout << "Remapping, found -99 for parent id " << hit->GetIparent() << ", " << map[hit->GetIparent()] << ", i_hit " << i_hit << endl;
+    hit->SetIparent(map[hit->GetIparent()]);
+    if (map[hit->GetPrimary()]==-99)
+      cout << "Remapping, found -99 for primary id " << hit->GetPrimary() << ", " << map[hit->GetPrimary()] << ", i_hit " << i_hit << endl;
+    hit->SetPrimary(map[hit->GetPrimary()]);
+    i_hit++;
+  }
+}
+
+void AliEMCALv2::FinishPrimary() {
+  fCurPrimary=-1;
+  fCurParent=-1;
+  fCurTrack=-1;
+}
+*/
diff --git a/EMCAL/AliEMCALv2.h b/EMCAL/AliEMCALv2.h
new file mode 100644 (file)
index 0000000..f3e2992
--- /dev/null
@@ -0,0 +1,75 @@
+#ifndef ALIEMCALV2_H
+#define ALIEMCALV2_H
+/* Copyright(c) 1998-2004, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                      
+         */
+/* $Id$ */
+
+//_________________________________________________________________________
+// Implementation version v2 of EMCAL Manager class for Shish-Kebab case 
+//*--                  
+//*-- Author:  Aleksei Pavlinov
+
+// --- ROOT system ---
+class TClonesArray;
+class TLorentzVector;
+class TFile;
+class TH1F;
+
+class TBrowser;
+class AliEMCALGeometry;
+
+// --- AliRoot header files ---
+#include "AliEMCALv1.h"
+
+// for TRD2 case
+//#include "TGeant3.h"
+
+class AliEMCALv2 : public AliEMCALv1 {
+  
+public:
+
+  AliEMCALv2(void) ; 
+  AliEMCALv2(const char *name, const char *title="") ;
+  // cpy ctor: no implementation yet
+  // requested by the Coding Convention
+  AliEMCALv2(const AliEMCALv1 & emcal):AliEMCALv1(emcal) {
+    Fatal("cpy ctor", "not implemented") ;  }
+  virtual ~AliEMCALv2(void) ;
+  virtual void  AddHit( Int_t shunt, Int_t primary, Int_t track, Int_t iparent, Float_t ienergy,
+                       Int_t id, Float_t *hits, Float_t *p);
+
+  virtual void StepManager(void) ;
+  virtual void FinishEvent();
+
+  // Gives the version number 
+  virtual Int_t  IsVersion(void) const {return 2;}
+  virtual const TString Version(void)const {return TString("v2");}
+  //  virtual void RemapTrackHitIDs(Int_t *map);
+  //virtual void FinishPrimary();
+  // virtual void SetTimeCut(Float_t tc){ fTimeCut = tc;}
+  // virtual Float_t GetTimeCut(){return fTimeCut;}
+  // assignement operator requested by coding convention but not needed  
+  AliEMCALv2 & operator = (const AliEMCALv1 & /*rvalue*/){
+    Fatal("operator =", "not implemented") ;  
+    return *this;}
+  // 23-mar-05
+  Double_t GetDepositEnergy(int print=1); // *MENU*
+  // 30-aug-04
+  virtual void Browse(TBrowser* b);
+  // drawing
+  void DrawCalorimeterCut(const char *name="SMOD", int axis=3, double dcut=1.); // *MENU*
+  void DrawSuperModuleCut(const char *name="EMOD", int axis=2, double dcut=0.03, int fill = 6);//  *MENU*
+  void DrawTowerCut(const char *name="SCMY", int axis=2, double dcut=0., int fill=1, char *optShad="on");   //  *MENU*
+  void DrawAlicWithHits(int mode=1);                            // *MENU*
+  void SetVolumeAttributes(const char *name="SCM0",const int seen=1, const int color=1, const int fill=1); // *MENU*
+  void TestIndexTransition(int pri=0, int idmax=0); // *MENU*
+
+  AliEMCALGeometry* fGeometry; //!
+  TH1F*             fHDe;      //!
+
+  ClassDef(AliEMCALv2,1)    //Implementation of EMCAL manager class to produce hits in a Shish-Kebab
+    
+};
+
+#endif // AliEMCALV2_H
index 14aa33f..f6a0b55 100644 (file)
@@ -16,6 +16,7 @@
 #pragma link C++ class AliEMCALGeometry+; 
 #pragma link C++ class AliEMCALv0+;
 #pragma link C++ class AliEMCALv1+;
+#pragma link C++ class AliEMCALv2+;
 #pragma link C++ class AliEMCALHit+;
 #pragma link C++ class AliEMCALDigit+;
 #pragma link C++ class AliEMCALTick+;
 #pragma link C++ class AliEMCALFastRecParticle+;               
 #pragma link C++ class AliEMCALPID+;           
 #pragma link C++ class AliEMCALPIDv1+;
-#pragma link C++ class AliEMCALTracker+;
 #pragma link C++ class AliEMCALLoader+;        
 #pragma link C++ class AliEMCALReconstructor+;
 #pragma link C++ class AliEMCALRawStream+;
+// 23-aug-04
+#pragma link C++ class AliEMCALGeneratorFactory+;
+// 9-sep-04
+#pragma link C++ class AliEMCALShishKebabModule+;
+#pragma link C++ class AliEMCALShishKebabTrd1Module+;
+#pragma link C++ class AliEMCALGeometryOfflineTrd1+;
+// 17-may-05
+#pragma link C++ class AliEMCALWsuCosmicRaySetUp+;
+//#pragma link C++ enum  PprRunFact_t;
+//#pragma link C++ enum  PprRadFact_t;
 
 #endif
index a2022e4..c10ee90 100644 (file)
@@ -5,6 +5,7 @@ AliEMCAL.cxx \
 AliEMCALGeometry.cxx \
 AliEMCALv0.cxx \
 AliEMCALv1.cxx \
+AliEMCALv2.cxx \
 AliEMCALHit.cxx \
 AliEMCALDigit.cxx \
 AliEMCALSDigitizer.cxx \
@@ -38,14 +39,18 @@ AliEMCALRecParticle.cxx \
 AliEMCALFastRecParticle.cxx \
 AliEMCALPID.cxx \
 AliEMCALPIDv1.cxx \
-AliEMCALTracker.cxx \
 AliEMCALLoader.cxx \
 AliEMCALReconstructor.cxx \
-AliEMCALRawStream.cxx
+AliEMCALRawStream.cxx \
+AliEMCALGeneratorFactory.cxx \
+AliEMCALShishKebabModule.cxx\
+AliEMCALShishKebabTrd1Module.cxx\
+AliEMCALGeometryOfflineTrd1.cxx\
+AliEMCALWsuCosmicRaySetUp.cxx
 
 HDRS= $(SRCS:.cxx=.h) 
 
 DHDR:=EMCALLinkDef.h
 
-EINCLUDE:=PYTHIA6 RAW
+EINCLUDE:=PYTHIA6 RAW EVGEN THijing