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]
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
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);
-
}
//____________________________________________________________________________
{
// 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 ");
}
}
-
-
-
-
-
virtual const TString Version() const {return TString(" ") ; }
AliEMCAL & operator = (const AliEMCAL & /*rvalue*/) {
Fatal("operator =", "not implemented") ; return *this ; }
-
+
protected:
friend class AliEMCALGetter;
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
#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;
}
//____________________________________________________________________________
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;
}
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;
printf("Exec took %f seconds for Clusterizing %f seconds per event",
gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nEvents ) ;
}
+
+ SaveHists();
+
}
//____________________________________________________________________________
fADCchannelECA = dig->GetECAchannel() ;
fADCpedestalECA = dig->GetECApedestal();
+ cout<<"ChannelECA, peds "<<fADCchannelECA<<" "<<fADCpedestalECA<<endl;
}
//____________________________________________________________________________
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;
}
//____________________________________________________________________________
{
// 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 ;
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] ) ;
}
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]) ;
branchECA->Fill() ;
- gime->WriteRecPoints("OVERWRITE");
- gime->WriteClusterizer("OVERWRITE");
+//Sub gime->WriteRecPoints("OVERWRITE");
+//Sub gime->WriteClusterizer("OVERWRITE");
}
//____________________________________________________________________________
// 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())
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 ;
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;
}
//____________________________________________________________________________
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];
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;
+ }
+}
+
// --- AliRoot header files ---
#include "AliEMCALClusterizer.h"
+#include "TH1F.h"
class AliEMCALRecPoint ;
class AliEMCALDigit ;
class AliEMCALDigitizer ;
// 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 ;
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)
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
//____________________________________________________________________________
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();
//____________________________________________________________________________
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();
// 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"
#include "AliEMCALSDigitizer.h"
#include "AliEMCALGeometry.h"
#include "AliEMCALTick.h"
+#include "AliEMCALJetMicroDst.h"
ClassImp(AliEMCALDigitizer)
// 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) ;
gime->Event(ReadEvent,"S");
sdigArray->AddAt(gime->SDigits(), i) ;
}
-
+
//Find the first tower with signal
Int_t nextSig = nEMC + 1 ;
TClonesArray * sdigits ;
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
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
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 ;
}
// 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 ;
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);
}
//____________________________________________________________________________
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
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")){
//____________________________________________________________________________
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;
}
//__________________________________________________________________
fEventNames[fInput] = eventFolderName ;
fInput++ ;
}
-
+
+void AliEMCALDigitizer::Print1(Option_t * option)
+{ // 19-nov-04 - just for convinience
+ Print();
+ PrintDigits(option);
+}
+
//__________________________________________________________________
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
}
}
}
+ printf("\n");
}
//__________________________________________________________________
// 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();
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);
+}
//
//*-- Author: Sahal Yacoob (LBL)
// based on : AliPHOSDigit
-// July 2003 Yves Schutz : NewIO
+// July 2003 Yves Schutz : NewIO
+// November 2003 Aleksei Pavlinov : Shish-Kebab geometry
//_________________________________________________________________________
#include "TObjString.h"
class TArrayI ;
class TClonesArray ;
+class TList;
+class TBrowser;
// --- Standard library ---
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
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();
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
-
};
--- /dev/null
+/**************************************************************************
+ * 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;
+}
--- /dev/null
+#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
//*-- Author: Sahal Yacoob (LBL / UCT)
// and : Yves Schutz (SUBATECH)
// and : Jennifer Klay (LBL)
+// SHASHLYK : Aleksei Pavlinov (WSU)
// --- AliRoot header files ---
#include <TMath.h>
AliEMCALGeometry *AliEMCALGeometry::fgGeom = 0;
Bool_t AliEMCALGeometry::fgInit = kFALSE;
+TString name; // contains name of geometry
//______________________________________________________________________
AliEMCALGeometry::~AliEMCALGeometry(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
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
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() );
}
}
}
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
+}
#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$ */
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));
}
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
--- /dev/null
+/**************************************************************************
+ * 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;
+}
--- /dev/null
+#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
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;
+}
class TBrowser;
class AliEMCALJetMicroDst: public TNamed {
-
-
public:
AliEMCALJetMicroDst(const char *name="jetMicroDst",
const char *tit="jet Micro Dst for preparation of proposal");
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 ; }
};
#endif // AliEMCALJETMICRODST_H
+
+typedef AliEMCALJetMicroDst sv; // for convinience
+
/*
What to do
1. Common info about event
ClassImp(AliEMCALRecPoint)
-
//____________________________________________________________________________
AliEMCALRecPoint::AliEMCALRecPoint()
: AliRecPoint()
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] ) ;
// 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");
}
//____________________________________________________________________________
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() ;
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) ;
d = 0. ;
fDispersion = TMath::Sqrt(d) ;
-
+ printf("Dispersion: = %f", fDispersion) ;
}
//____________________________________________________________________________
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 );
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 ;
}
AliEMCALDigit * digit ;
AliEMCALGeometry * geom = (AliEMCALGetter::Instance())->EMCALGeometry();
+ TString gn(geom->GetName());
Int_t iDigit;
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 ;
}
fLambda[1]= 0. ;
}
+ // printf("Evalaxis: lambdas = %f,%f", fLambda[0],fLambda[1]) ;
}
#include "AliEMCALClusterizerv1.h"
#include "AliEMCALPIDv1.h"
#include "AliEMCALGetter.h"
-#include "AliEMCALTracker.h"
#include "AliRawReaderFile.h"
ClassImp(AliEMCALReconstructor)
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();
delete et;
}
}
-
-AliTracker* AliEMCALReconstructor::CreateTracker(AliRunLoader* runLoader) const
-{
-// creates the EMCAL tracker
- if (!runLoader) return NULL;
- return new AliEMCALTracker(runLoader);
-}
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 ;
// --- ROOT system ---
#include "TBenchmark.h"
- //#include "TObjectTable.h"
+#include "TH1.h"
+#include "TBrowser.h"
+//#include "TObjectTable.h"
// --- Standard library ---
#include "stdlib.h"
#include "AliEMCALHit.h"
#include "AliEMCALSDigitizer.h"
#include "AliEMCALGeometry.h"
- //#include "AliMemoryWatcher.h"
+#include "AliEMCALJetMicroDst.h"
+//#include "AliMemoryWatcher.h"
ClassImp(AliEMCALSDigitizer)
AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("","")
{
// ctor
- fFirstEvent = fLastEvent = 0 ;
+ Info("AliEMCALSDigitizer()", " CTOR # 1 #");
+ fFirstEvent = fLastEvent = fControlHists = 0;
fDefaultInit = kTRUE ;
+ fHists = 0;
}
//____________________________________________________________________________
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 ;
//____________________________________________________________________________
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 ;
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() ;
//=========== 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 ;
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()) ;
}
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();
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 ) ;
}
}
+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") ;
-
}
//__________________________________________________________________
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++) {
}
}
delete tempo ;
- }
+ printf("\n** Sum %i : %10.3f GeV/c **\n ", isum, double(isum)*1.e-6);
+ } else printf("\n");
}
//____________________________________________________________________________
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);
+}
// --- ROOT system ---
#include "TTask.h"
class TFile ;
+class TList;
+class TBrowser;
+//class TBrowser;
// --- Standard library ---
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() ;
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
--- /dev/null
+/**************************************************************************
+ * 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();
+ }
+ }
+}
+
--- /dev/null
+#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
+ */
--- /dev/null
+/**************************************************************************
+ * 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());
+ }
+ }
+}
--- /dev/null
+#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
+ */
--- /dev/null
+/**************************************************************************
+ * 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");
+}
+
+
--- /dev/null
+#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
// --- 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;
}
//______________________________________________________________________
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;
}
//______________________________________________________________________
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 ;
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;
+ }
+}
// --- ROOT system ---
class TFile;
+class TList;
+class TNode;
// --- AliRoot header files ---
#include "AliEMCAL.h"
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 {
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:
// 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
new((*fHits)[fNhits]) AliEMCALHit(*newHit);
fNhits++;
} // end if
-
+
delete newHit;
}
//______________________________________________________________________
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
};
--- /dev/null
+/**************************************************************************
+ * 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;
+}
+*/
--- /dev/null
+#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
#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
AliEMCALGeometry.cxx \
AliEMCALv0.cxx \
AliEMCALv1.cxx \
+AliEMCALv2.cxx \
AliEMCALHit.cxx \
AliEMCALDigit.cxx \
AliEMCALSDigitizer.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