From 1963b290d5415854f3712a5255c9d2dd7ec461ef Mon Sep 17 00:00:00 2001 From: pavlinov Date: Mon, 1 Aug 2005 16:11:19 +0000 Subject: [PATCH] "Version --- EMCAL/AliEMCAL.cxx | 48 +- EMCAL/AliEMCAL.h | 3 +- EMCAL/AliEMCALClusterizerv1.cxx | 219 ++++- EMCAL/AliEMCALClusterizerv1.h | 23 +- EMCAL/AliEMCALDigit.cxx | 4 +- EMCAL/AliEMCALDigitizer.cxx | 156 +++- EMCAL/AliEMCALDigitizer.h | 21 +- EMCAL/AliEMCALGeneratorFactory.cxx | 669 +++++++++++++++ EMCAL/AliEMCALGeneratorFactory.h | 79 ++ EMCAL/AliEMCALGeometry.cxx | 229 +++++- EMCAL/AliEMCALGeometry.h | 67 +- EMCAL/AliEMCALGeometryOfflineTrd1.cxx | 197 +++++ EMCAL/AliEMCALGeometryOfflineTrd1.h | 78 ++ EMCAL/AliEMCALJetMicroDst.cxx | 58 ++ EMCAL/AliEMCALJetMicroDst.h | 12 +- EMCAL/AliEMCALRecPoint.cxx | 88 +- EMCAL/AliEMCALReconstruct.C | 1 - EMCAL/AliEMCALReconstructor.cxx | 21 +- EMCAL/AliEMCALReconstructor.h | 3 +- EMCAL/AliEMCALSDigitizer.cxx | 153 ++-- EMCAL/AliEMCALSDigitizer.h | 23 +- EMCAL/AliEMCALShishKebabModule.cxx | 183 +++++ EMCAL/AliEMCALShishKebabModule.h | 63 ++ EMCAL/AliEMCALShishKebabTrd1Module.cxx | 155 ++++ EMCAL/AliEMCALShishKebabTrd1Module.h | 75 ++ EMCAL/AliEMCALWsuCosmicRaySetUp.cxx | 132 +++ EMCAL/AliEMCALWsuCosmicRaySetUp.h | 29 + EMCAL/AliEMCALv0.cxx | 1030 +++++++++++++++++++++++- EMCAL/AliEMCALv0.h | 20 + EMCAL/AliEMCALv1.cxx | 4 +- EMCAL/AliEMCALv1.h | 5 +- EMCAL/AliEMCALv2.cxx | 569 +++++++++++++ EMCAL/AliEMCALv2.h | 75 ++ EMCAL/EMCALLinkDef.h | 12 +- EMCAL/libEMCAL.pkg | 11 +- 35 files changed, 4278 insertions(+), 237 deletions(-) create mode 100644 EMCAL/AliEMCALGeneratorFactory.cxx create mode 100644 EMCAL/AliEMCALGeneratorFactory.h create mode 100644 EMCAL/AliEMCALGeometryOfflineTrd1.cxx create mode 100644 EMCAL/AliEMCALGeometryOfflineTrd1.h create mode 100644 EMCAL/AliEMCALShishKebabModule.cxx create mode 100644 EMCAL/AliEMCALShishKebabModule.h create mode 100644 EMCAL/AliEMCALShishKebabTrd1Module.cxx create mode 100644 EMCAL/AliEMCALShishKebabTrd1Module.h create mode 100644 EMCAL/AliEMCALWsuCosmicRaySetUp.cxx create mode 100644 EMCAL/AliEMCALWsuCosmicRaySetUp.h create mode 100644 EMCAL/AliEMCALv2.cxx create mode 100644 EMCAL/AliEMCALv2.h diff --git a/EMCAL/AliEMCAL.cxx b/EMCAL/AliEMCAL.cxx index 2822ed5ff77..044c6646fc9 100644 --- a/EMCAL/AliEMCAL.cxx +++ b/EMCAL/AliEMCAL.cxx @@ -121,6 +121,12 @@ void AliEMCAL::CreateMaterials() AliMaterial(3, "Al$", 26.98, 13., 2.7, 8.9, 999., 0, 0) ; // --- Absorption length is ignored ^ + // 25-aug-04 by PAI - see PMD/AliPMDv0.cxx for STEEL definition + Float_t asteel[4] = { 55.847,51.9961,58.6934,28.0855 }; + Float_t zsteel[4] = { 26.,24.,28.,14. }; + Float_t wsteel[4] = { .715,.18,.1,.005 }; + AliMixture(4, "STAINLESS STEEL$", asteel, zsteel, 7.88, 4, wsteel); + // DEFINITION OF THE TRACKING MEDIA // for EMCAL: idtmed[1599->1698] equivalent to fIdtmed[0->100] @@ -145,6 +151,9 @@ void AliEMCAL::CreateMaterials() AliMedium(3, "Al parts $", 3, 0, isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.001, 0.001, 0, 0) ; + // 25-aug-04 by PAI : see PMD/AliPMDv0.cxx for STEEL definition -> idtmed[1603] + AliMedium(4, "S steel$", 4, 0, + isxfld, sxmgmx, 10.0, 0.1, 0.1, 0.001, 0.001, 0, 0) ; // --- Set decent energy thresholds for gamma and electron tracking @@ -170,12 +179,21 @@ void AliEMCAL::CreateMaterials() gMC->Gstpar(idtmed[1601],"CUTELE",0.001) ; gMC->Gstpar(idtmed[1601],"BCUTE",0.0001) ; + // the same parameters as for Pb - 10-sep-04 + gMC->Gstpar(idtmed[1603],"CUTGAM",0.00008) ; + gMC->Gstpar(idtmed[1603],"CUTELE",0.001) ; + gMC->Gstpar(idtmed[1603],"BCUTE", 0.0001); + // --- Generate explicitly delta rays in Lead --- + gMC->Gstpar(idtmed[1603], "LOSS",3.); + gMC->Gstpar(idtmed[1603], "DRAY",1.); + gMC->Gstpar(idtmed[1603], "DCUTE",0.00001); + gMC->Gstpar(idtmed[1603], "DCUTM",0.00001); + //set constants for Birk's Law implentation fBirkC0 = 1; fBirkC1 = 0.013/dP; fBirkC2 = 9.6e-6/(dP * dP); - } //____________________________________________________________________________ @@ -365,27 +383,21 @@ void AliEMCAL::SetTreeAddress() { // Linking Hits in Tree to Hits array TBranch *branch; - char branchname[20]; - sprintf(branchname,"%s",GetName()); - + // char branchname[20]; + // sprintf(branchname,"%s",GetName()); // Branch address for hit tree TTree *treeH = TreeH(); if (treeH) { - branch = treeH->GetBranch(branchname); - if (branch) - { + // treeH->Print(); + branch = treeH->GetBranch(GetName()); + if (branch) { if (fHits == 0x0) - fHits= new TClonesArray("AliEMCALHit",1000); + fHits= new TClonesArray("AliEMCALHit",1000); branch->SetAddress(&fHits); - } - else - { - Warning("SetTreeAddress","<%s> Failed",GetName()); - } + } else { + Warning("SetTreeAddress","<%s> Failed",GetName()); + } + } else { + // Warning("SetTreeAddress"," no treeH "); } } - - - - - diff --git a/EMCAL/AliEMCAL.h b/EMCAL/AliEMCAL.h index 13a103e858e..36093e43392 100644 --- a/EMCAL/AliEMCAL.h +++ b/EMCAL/AliEMCAL.h @@ -65,7 +65,7 @@ class AliEMCAL : public AliDetector { virtual const TString Version() const {return TString(" ") ; } AliEMCAL & operator = (const AliEMCAL & /*rvalue*/) { Fatal("operator =", "not implemented") ; return *this ; } - + protected: friend class AliEMCALGetter; @@ -76,6 +76,7 @@ protected: Int_t fBirkC0; // constants for Birk's Law implementation Double_t fBirkC1; // constants for Birk's Law implementation Double_t fBirkC2; // constants for Birk's Law implementation + static Double_t fgCapa ; // capacitor of the preamplifier for the raw RO signal Double_t fHighCharge ; // high charge (to convert energy to charge) for the raw RO signal Double_t fHighGain ; // high gain for the raw RO signal diff --git a/EMCAL/AliEMCALClusterizerv1.cxx b/EMCAL/AliEMCALClusterizerv1.cxx index 951c97aadbd..bfb8d83dd8c 100644 --- a/EMCAL/AliEMCALClusterizerv1.cxx +++ b/EMCAL/AliEMCALClusterizerv1.cxx @@ -69,14 +69,19 @@ #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"<MaxEvent() - 1; - + //fLastEvent = gime->MaxEvent() - 1; + fLastEvent = 1000 - 1; Int_t nEvents = fLastEvent - fFirstEvent + 1; Int_t ievent ; for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) { gime->Event(ievent,"D") ; - GetCalibrationParameters() ; fNumberOfECAClusters = 0; @@ -162,6 +176,9 @@ void AliEMCALClusterizerv1::Exec(Option_t * option) printf("Exec took %f seconds for Clusterizing %f seconds per event", gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nEvents ) ; } + + SaveHists(); + } //____________________________________________________________________________ @@ -267,6 +284,7 @@ void AliEMCALClusterizerv1::GetCalibrationParameters() fADCchannelECA = dig->GetECAchannel() ; fADCpedestalECA = dig->GetECApedestal(); + cout<<"ChannelECA, peds "<EMCALGeometry() ; + cout<<"gime,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 "<AbsToRelNumbering(d1->GetId(), relid1) ; - - Int_t relid2[2] ; - geom->AbsToRelNumbering(d2->GetId(), relid2) ; + //Sub geom->AbsToRelNumbering(d1->GetId(), relid1) ; + Int_t nSupMod=0; + Int_t nTower=0; + Int_t nIphi=0; + Int_t nIeta=0; + Int_t iphi=0; + Int_t ieta=0; + geom->GetCellIndex(d1->GetId(), nSupMod,nTower,nIphi,nIeta); + geom->GetCellPhiEtaIndexInSModule(nTower,nIphi,nIeta, iphi,ieta); + relid1[0]=ieta; + relid1[1]=iphi; + + + Int_t nSupMod1=0; + Int_t nTower1=0; + Int_t nIphi1=0; + Int_t nIeta1=0; + Int_t iphi1=0; + Int_t ieta1=0; + geom->GetCellIndex(d2->GetId(), nSupMod1,nTower1,nIphi1,nIeta1); + geom->GetCellPhiEtaIndexInSModule(nTower1,nIphi1,nIeta1, iphi1,ieta1); + Int_t relid2[2] ; + relid2[0]=ieta1; + relid2[1]=iphi1; + + //Sub geom->AbsToRelNumbering(d2->GetId(), relid2) ; Int_t rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ; @@ -336,7 +381,7 @@ Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d } if (gDebug == 2 ) - printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d ", +if(rv==1)printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n", rv, d1->GetId(), relid1[0], relid1[1], d2->GetId(), relid2[0], relid2[1]) ; @@ -388,8 +433,8 @@ void AliEMCALClusterizerv1::WriteRecPoints() branchECA->Fill() ; - gime->WriteRecPoints("OVERWRITE"); - gime->WriteClusterizer("OVERWRITE"); +//Sub gime->WriteRecPoints("OVERWRITE"); +//Sub gime->WriteClusterizer("OVERWRITE"); } //____________________________________________________________________________ @@ -412,18 +457,38 @@ void AliEMCALClusterizerv1::MakeClusters() // Clusterization starts TIter nextdigit(digitsC) ; AliEMCALDigit * digit; - + +//just for hist +/* + while ( (digit = dynamic_cast(nextdigit())) ) { // scan over the list of digitsC + digitAmp->Fill(digit->GetAmp()); + } + */ +/////////// + while ( (digit = dynamic_cast(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<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 "<GetAmp()<= aECARecPoints->GetSize()) @@ -432,12 +497,13 @@ void AliEMCALClusterizerv1::MakeClusters() aECARecPoints->AddAt(rp, fNumberOfECAClusters) ; clu = dynamic_cast(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 "<Remove(digit) ; if (gDebug == 2 ) - printf("MakeClusters: OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp()), fECAClusteringThreshold) ; + printf("MakeClusters: OK id = %d, ene = %f , thre = %f \n", digit->GetId(),Calibrate(digit->GetAmp()), fECAClusteringThreshold) ; nextdigit.Reset() ; AliEMCALDigit * digitN ; @@ -448,32 +514,46 @@ void AliEMCALClusterizerv1::MakeClusters() digit = (AliEMCALDigit*)digits->At(clusterECAdigitslist[index]) ; index++ ; while ( (digitN = (AliEMCALDigit *)nextdigit())) { // scan over the reduced list of digits +// cout<<"we have new digit "<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 "<AddDigit(*digitN, Calibrate( digitN->GetAmp()) ) ; +if(addOn[nSupMod-1][ieta-1][iphi-1]>0)cout<<"22 digit, add "<GetAmp()<AddDigit(*digitN, Calibrate( digitN->GetAmp()+addOn[nSupMod-1][ieta-1][iphi-1]) ) ; clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ; iDigitInECACluster++ ; +// cout<<"when neighbour: cluno, digNo "<GetId()<<" "<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"<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 "<GetAmp()<Remove(digit); } + //cout<<"after endofloop: cluno, digNo "<GetEntriesFast()<<" digits"<GetEntries() ; index++) { AliEMCALRecPoint * rp = dynamic_cast(aECARecPoints->At(index)) ; TVector3 globalpos; - rp->GetGlobalPosition(globalpos); + //rp->GetGlobalPosition(globalpos); TVector3 localpos; rp->GetLocalPosition(localpos); Float_t lambda[2]; @@ -588,10 +672,73 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ; + ///////////// + if(rp->GetEnergy()>maxE){ + maxE=rp->GetEnergy(); + maxL1=lambda[0]; + maxL2=lambda[1]; + maxDis=rp->GetDispersion(); + } + pointE->Fill(rp->GetEnergy()); + pointL1->Fill(lambda[0]); + pointL2->Fill(lambda[1]); + pointDis->Fill(rp->GetDispersion()); + pointMult->Fill(rp->GetMultiplicity()); + ///////////// for (Int_t iprimary=0; iprimaryFill(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<EMCALGeometry(); @@ -163,7 +163,7 @@ Float_t AliEMCALDigit::GetEta() const //____________________________________________________________________________ Float_t AliEMCALDigit::GetPhi() const -{ +{ // should be change in EMCALGeometry - 19-nov-04 Float_t eta=-10., phi=-10.; Int_t id = GetId(); const AliEMCALGeometry *g = AliEMCALGetter::Instance()->EMCALGeometry(); diff --git a/EMCAL/AliEMCALDigitizer.cxx b/EMCAL/AliEMCALDigitizer.cxx index ac37421e64a..54167f0ee8f 100644 --- a/EMCAL/AliEMCALDigitizer.cxx +++ b/EMCAL/AliEMCALDigitizer.cxx @@ -55,14 +55,18 @@ // Modif: // August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction // of new IO (à la PHOS) -///_________________________________________________________________________________ +// November 2003 Aleksei Pavlinov : adopted for Shish-Kebab geometry +//_________________________________________________________________________________ // --- ROOT system --- +#include "TROOT.h" #include "TTree.h" #include "TSystem.h" #include "TBenchmark.h" - -// --- Standard library --- +#include "TList.h" +#include "TH1.h" +#include "TBrowser.h" +#include "TObjectTable.h" // --- AliRoot header files --- #include "AliRunDigitizer.h" @@ -73,6 +77,7 @@ #include "AliEMCALSDigitizer.h" #include "AliEMCALGeometry.h" #include "AliEMCALTick.h" +#include "AliEMCALJetMicroDst.h" ClassImp(AliEMCALDigitizer) @@ -160,20 +165,35 @@ void AliEMCALDigitizer::Digitize(Int_t event) // after that adds contributions from SDigits. This design // helps to avoid scanning over the list of digits to add // contribution of any new SDigit. + static int isTrd1Geom = -1; // -1 - mean undefined + static int nEMC=0; //max number of digits possible + //<<<<<<< AliEMCALDigitizer.cxx + int deb=0; + AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName) ; + /* ======= AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle()) ; + >>>>>>> 1.59 */ Int_t ReadEvent = event ; + // fManager is data member from AliDigitizer if (fManager) ReadEvent = dynamic_cast(fManager->GetInputStream(0))->GetCurrentEventNumber() ; - Info("Digitize", "Adding event %d from input stream 0 %s %s", ReadEvent, GetTitle(), fEventFolderName.Data()) ; + if(deb>0) Info("Digitize", "Adding event %d from input stream 0 %s %s", + ReadEvent, GetTitle(), fEventFolderName.Data()) ; gime->Event(ReadEvent, "S") ; TClonesArray * digits = gime->Digits() ; digits->Clear() ; - const AliEMCALGeometry *geom = gime->EMCALGeometry() ; + const AliEMCALGeometry *geom = gime->EMCALGeometry() ; + if(isTrd1Geom < 0) { + TString ng(geom->GetName()); + isTrd1Geom = 0; + if(ng.Contains("SHISH") && ng.Contains("TRD1")) isTrd1Geom = 1; - Int_t nEMC = geom->GetNPhi()*geom->GetNZ(); //max number of digits possible - + if(isTrd1Geom == 0) nEMC = geom->GetNPhi()*geom->GetNZ(); + else nEMC = geom->GetNCells(); + printf(" nEMC %i (number cells in EMCAL) | %s | isTrd1Geom %i\n", nEMC, geom->GetName(), isTrd1Geom); + } Int_t absID ; digits->Expand(nEMC) ; @@ -200,7 +220,7 @@ void AliEMCALDigitizer::Digitize(Int_t event) gime->Event(ReadEvent,"S"); sdigArray->AddAt(gime->SDigits(), i) ; } - + //Find the first tower with signal Int_t nextSig = nEMC + 1 ; TClonesArray * sdigits ; @@ -211,7 +231,9 @@ void AliEMCALDigitizer::Digitize(Int_t event) Int_t curNext = dynamic_cast(sdigits->At(0))->GetId() ; if(curNext < nextSig) nextSig = curNext ; + if(deb>0) printf("input %i : #sdigits %i \n", i, sdigits->GetEntriesFast()); } + if(deb>0) printf("FIRST tower with signal %i \n", nextSig); TArrayI index(fInput) ; index.Reset() ; //Set all indexes to zero @@ -221,7 +243,7 @@ void AliEMCALDigitizer::Digitize(Int_t event) TClonesArray * ticks = new TClonesArray("AliEMCALTick",1000) ; - //Put Noise contribution + //Put Noise contributiony for(absID = 1; absID <= nEMC; absID++){ Float_t amp = 0 ; // amplitude set to zero, noise will be added later @@ -265,7 +287,7 @@ void AliEMCALDigitizer::Digitize(Int_t event) index[i]++ ; if( dynamic_cast(sdigArray->At(i))->GetEntriesFast() > index[i] ) - curSDigit = dynamic_cast(dynamic_cast(sdigArray->At(i))->At(index[i])) ; + curSDigit = dynamic_cast(dynamic_cast(sdigArray->At(i))->At(index[i])) ; else curSDigit = 0 ; } @@ -292,7 +314,8 @@ void AliEMCALDigitizer::Digitize(Int_t event) // add the noise now amp += TMath::Abs(gRandom->Gaus(0., fPinNoise)) ; digit->SetAmp(sDigitizer->Digitize(amp)) ; - } + if(deb>=10) printf(" absID %5i amp %f nextSig %5i\n", absID, amp, nextSig); + } // for(absID = 1; absID <= nEMC; absID++) ticks->Delete() ; delete ticks ; @@ -314,13 +337,20 @@ void AliEMCALDigitizer::Digitize(Int_t event) Int_t ndigits = digits->GetEntriesFast() ; digits->Expand(ndigits) ; - //Set indexes in list of digits + //Set indexes in list of digits and fill hists. + sv::FillH1(fHists, 0, Double_t(ndigits)); + Float_t energy=0., esum=0.; for (i = 0 ; i < ndigits ; i++) { digit = dynamic_cast( digits->At(i) ) ; digit->SetIndexInList(i) ; - Float_t energy = sDigitizer->Calibrate(digit->GetAmp()) ; - digit->SetAmp(DigitizeEnergy(energy) ) ; + energy = sDigitizer->Calibrate(digit->GetAmp()) ; + esum += energy; + digit->SetAmp(DigitizeEnergy(energy) ) ; // for what ?? + sv::FillH1(fHists, 2, double(digit->GetAmp())); + sv::FillH1(fHists, 3, double(energy)); + sv::FillH1(fHists, 4, double(digit->GetId())); } + sv::FillH1(fHists, 1, esum); } //____________________________________________________________________________ @@ -364,14 +394,17 @@ void AliEMCALDigitizer::Exec(Option_t *option) if (fLastEvent == -1) fLastEvent = gime->MaxEvent() - 1 ; else if (fManager) - fLastEvent = fFirstEvent ; + fLastEvent = fFirstEvent ; // what is this ?? Int_t nEvents = fLastEvent - fFirstEvent + 1; - - Int_t ievent ; - - for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) { - + Int_t ievent, nfr=50; + + fLastEvent = TMath::Min(fLastEvent, gime->MaxEvent()); + printf("AliEMCALDigitizer::Exec : option: %s | %i -> %i events : Max events %i \n", + option, fFirstEvent, fLastEvent, gime->MaxEvent()); + for (ievent = fFirstEvent; ievent < fLastEvent; ievent++) { + if(ievent%nfr==0 || ievent==fLastEvent-1); + printf(" processed event %i\n", ievent); gime->Event(ievent,"S") ; Digitize(ievent) ; //Add prepared SDigits to digits and add the noise @@ -380,11 +413,13 @@ void AliEMCALDigitizer::Exec(Option_t *option) if(strstr(option,"deb")) PrintDigits(option); + if(strstr(option,"table")) gObjectTable->Print(); //increment the total number of Digits per run fDigitsInRun += gime->Digits()->GetEntriesFast() ; } - + // gime->WriteDigitizer("OVERWRITE"); + gime->EmcalLoader()->CleanDigitizer() ; if(strstr(option,"tim")){ @@ -459,21 +494,24 @@ Bool_t AliEMCALDigitizer::Init() //____________________________________________________________________________ void AliEMCALDigitizer::InitParameters() -{ - fMeanPhotonElectron = 18200 ; // electrons per GeV - fPinNoise = 0.003 ; // noise equivalent GeV (random choice) +{ // Tune parameters - 24-nov-04 + + fMeanPhotonElectron = 3300 ; // electrons per GeV + fPinNoise = 0.004; if (fPinNoise == 0. ) Warning("InitParameters", "No noise added\n") ; - fDigitThreshold = fPinNoise * 3; //2 sigma - fTimeResolution = 1.0e-9 ; + fDigitThreshold = fPinNoise * 3; // 3 * sigma + fTimeResolution = 0.3e-9 ; // 300 psc fTimeSignalLength = 1.0e-9 ; - fADCchannelEC = 0.00305; // width of one ADC channel in GeV - HG fix so that we see 200 GeV gammas - fADCpedestalEC = 0.005 ; // GeV - fNADCEC = (Int_t) TMath::Power(2,16) ; // number of channels in Tower ADC + fADCchannelEC = 0.00305; // 200./65536 - width of one ADC channel in GeV + fADCpedestalEC = 0.009 ; // GeV + fNADCEC = (Int_t) TMath::Power(2,16) ; // number of channels in Tower ADC - 65536 - fTimeThreshold = 0.001*10000000 ; //Means 1 MeV in terms of SDigits amplitude - + fTimeThreshold = 0.001*10000000 ; // Means 1 MeV in terms of SDigits amplitude ?? + // hists. for control; no hists on default + fControlHists = 0; + fHists = 0; } //__________________________________________________________________ @@ -527,7 +565,13 @@ void AliEMCALDigitizer::MixWith(TString alirunFileName, TString eventFolderName) fEventNames[fInput] = eventFolderName ; fInput++ ; } - + +void AliEMCALDigitizer::Print1(Option_t * option) +{ // 19-nov-04 - just for convinience + Print(); + PrintDigits(option); +} + //__________________________________________________________________ void AliEMCALDigitizer::Print()const { @@ -569,11 +613,11 @@ void AliEMCALDigitizer::Print()const void AliEMCALDigitizer::PrintDigits(Option_t * option){ AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName) ; - TClonesArray * digits = gime->Digits() ; + TClonesArray * digits = gime->Digits() ; + TClonesArray * sdigits = gime->SDigits() ; - printf("PrintDigits: %d", digits->GetEntriesFast()) ; - printf("\nevent %d", gAlice->GetEvNumber()) ; - printf("\n Number of entries in Digits list %d", digits->GetEntriesFast() ) ; + printf("\n #Digits: %d : sdigits %d ", digits->GetEntriesFast(), sdigits->GetEntriesFast()) ; + printf("\n event %d", gAlice->GetEvNumber()) ; if(strstr(option,"all")){ //loop over digits @@ -591,6 +635,7 @@ void AliEMCALDigitizer::PrintDigits(Option_t * option){ } } } + printf("\n"); } //__________________________________________________________________ @@ -629,7 +674,12 @@ void AliEMCALDigitizer::WriteDigits() // and branch "AliEMCALDigitizer", with the same title to keep all the parameters // and names of files, from which digits are made. + //<<<<<<< AliEMCALDigitizer.cxx + static Int_t writeSdigitizer=0; + AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName) ; + /* ======= AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle()) ; + >>>>>>> 1.59*/ const TClonesArray * digits = gime->Digits() ; TTree * treeD = gime->TreeD(); @@ -640,9 +690,41 @@ void AliEMCALDigitizer::WriteDigits() digitsBranch->Fill() ; gime->WriteDigits("OVERWRITE"); - gime->WriteDigitizer("OVERWRITE"); + if(writeSdigitizer==0) { + gime->WriteDigitizer("OVERWRITE"); + writeSdigitizer = 1; + } Unload() ; } +void AliEMCALDigitizer::Browse(TBrowser* b) +{ + if(fHists) b->Add(fHists); + TTask::Browse(b); +} + +TList *AliEMCALDigitizer::BookControlHists(const int var) +{ // 22-nov-04 + Info("BookControlHists"," started "); + gROOT->cd(); + AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName); + const AliEMCALGeometry *geom = gime->EMCALGeometry() ; + if(var>=1){ + new TH1F("hDigiN", "#EMCAL digits with fAmp > fDigitThreshold", + fNADCEC+1, -0.5, Double_t(fNADCEC)); + new TH1F("HDigiSumEnergy","Sum.EMCAL energy from digi", 1000, 0.0, 200.); + new TH1F("hDigiAmp", "EMCAL digital amplitude", fNADCEC+1, -0.5, Double_t(fNADCEC)); + new TH1F("hDigiEnergy","EMCAL cell energy", 2000, 0.0, 200.); + new TH1F("hDigiAbsId","EMCAL absId cells with fAmp > fDigitThreshold ", + geom->GetNCells(), 0.5, Double_t(geom->GetNCells())+0.5); + } + fHists = sv::MoveHistsToList("EmcalDigiControlHists", kFALSE); + return fHists; +} + +void AliEMCALDigitizer::SaveHists(const char* name, Bool_t kSingleKey, const char* opt) +{ + sv::SaveListOfHists(fHists, name, kSingleKey, opt); +} diff --git a/EMCAL/AliEMCALDigitizer.h b/EMCAL/AliEMCALDigitizer.h index 970b82c4bf7..61e2874a0f0 100644 --- a/EMCAL/AliEMCALDigitizer.h +++ b/EMCAL/AliEMCALDigitizer.h @@ -10,7 +10,8 @@ // //*-- Author: Sahal Yacoob (LBL) // based on : AliPHOSDigit -// July 2003 Yves Schutz : NewIO +// July 2003 Yves Schutz : NewIO +// November 2003 Aleksei Pavlinov : Shish-Kebab geometry //_________________________________________________________________________ @@ -18,6 +19,8 @@ #include "TObjString.h" class TArrayI ; class TClonesArray ; +class TList; +class TBrowser; // --- Standard library --- @@ -54,7 +57,8 @@ public: Int_t GetDigitsInRun() const { return fDigitsInRun; } void MixWith(TString alirunFileName, TString eventFolderName = AliConfig::GetDefaultEventFolderName()) ; // Add another one file to mix - void Print()const ; + void Print() const ; + void Print1(Option_t * option); // *MENU* AliEMCALDigitizer & operator = (const AliEMCALDigitizer & /*rvalue*/) { // assignement operator requested by coding convention but not needed @@ -62,6 +66,15 @@ public: return *this ; } + virtual void Browse(TBrowser* b); + // hists + void SetControlHists(const Int_t var=0) {fControlHists=var;} + Int_t GetControlHist() const {return fControlHists;} + TList *GetListOfHists() {return fHists;} + TList* BookControlHists(const int var=0); + void SaveHists(const char* name="RF/TRD1/Digitizations/DigiVar?", + Bool_t kSingleKey=kTRUE, const char* opt="RECREATE"); // *MENU* + private: Bool_t Init(); @@ -100,9 +113,11 @@ private: TString fEventFolderName; // skowron: name of EFN to read data from in stand alone mode Int_t fFirstEvent; // first event to process Int_t fLastEvent; // last event to process + // Control hists + Int_t fControlHists; //! + TList *fHists; //! ClassDef(AliEMCALDigitizer,5) // description - }; diff --git a/EMCAL/AliEMCALGeneratorFactory.cxx b/EMCAL/AliEMCALGeneratorFactory.cxx new file mode 100644 index 00000000000..012fbbbecd7 --- /dev/null +++ b/EMCAL/AliEMCALGeneratorFactory.cxx @@ -0,0 +1,669 @@ +/************************************************************************** + * Copyright(c) 1998-2002, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* +$Log$ +*/ +#include "assert.h" +#include "TString.h" +#include "AliEMCALGeneratorFactory.h" +#include "AliGenerator.h" +#include "AliGenFixed.h" +#include "AliGenBox.h" +#include "AliGenHIJINGpara.h" +#include "AliGenHIJINGparaBa.h" +#include "AliGenHijing.h" +#include "AliGenCocktail.h" +#include "AliGenPythia.h" + +//*-- Authors: Aleksei Pavlinov (WSU) +//* Initial variant is in Config.C for EMCAL production. +ClassImp(AliEMCALGeneratorFactory) + +AliEMCALGeneratorFactory::AliEMCALGeneratorFactory(PprRunFact_t run, PprRadFact_t rad) +{ + fGenerator = fBgGenerator = fSignalGenerator = 0; + fRunType = run; + fRadiation = rad; + fMomentum = 0; + + Int_t isw = 3; + if (rad == kNoGluonRadiation) isw = 0; + Float_t thmin=0, thmax=0; + + AliGenFixed *genGG=0; + AliGenBox *genB=0; + AliGenHIJINGparaBa *genHijParaB=0, *bg=0; + AliGenHIJINGpara *genHijPara=0; + AliGenHijing *genHij=0; + AliGenCocktail *genCoct=0; + AliGenPythia *genPy=0, *jets=0; + // AliPythia *aliPy = 0; + + fComment = new TString; + + switch (run) { + case kGammaGun: + genGG = new AliGenFixed(1); + genGG->SetMomentum(100.); + genGG->SetPhi(240.); + genGG->SetTheta(91.); + genGG->SetPart(kGamma); + fGenerator = (AliGenerator*)genGG; + break; + + case kGammaBox: + genB = new AliGenBox(100); + genB->SetMomentumRange(100., 100.1); + genB->SetPhiRange(0,360); + genB->SetThetaRange(45., 135.); + genB->SetPart(kGamma); + fGenerator = (AliGenerator*)genB; + break; + + case kTest50: + genHijParaB = new AliGenHIJINGparaBa(50); + genHijParaB->SetMomentumRange(0, 999999.); + genHijParaB->SetPhiRange(-180., 180.); + // Set pseudorapidity range from -8 to 8. + thmin = EtaToTheta(8); // theta min. <---> eta max + thmax = EtaToTheta(-8); // theta max. <---> eta min + genHijParaB->SetThetaRange(thmin,thmax); + fGenerator = (AliGenerator*)genHijParaB; + break; + + case kParam_8000: + //coment= fComment.Append(":HIJINGparam N=8000"); + genHijPara = new AliGenHIJINGpara(86030); + genHijPara->SetMomentumRange(0, 999999.); + genHijPara->SetPhiRange(-180., 180.); + // Set pseudorapidity range from -8 to 8. + thmin = EtaToTheta(8); // theta min. <---> eta max + thmax = EtaToTheta(-8); // theta max. <---> eta min + genHijPara->SetThetaRange(thmin,thmax); + fGenerator = (AliGenerator*)genHijPara; + break; + case kParam_4000: + genHijPara = new AliGenHIJINGpara(43015); + genHijPara->SetMomentumRange(0, 999999.); + genHijPara->SetPhiRange(-180., 180.); + // Set pseudorapidity range from -8 to 8. + thmin = EtaToTheta(8); // theta min. <---> eta max + thmax = EtaToTheta(-8); // theta max. <---> eta min + genHijPara->SetThetaRange(thmin,thmax); + fGenerator = (AliGenerator*)genHijPara; + break; + case kParam_2000: + (*fComment) = "HIJINGparam N=2000"; + genHijPara = new AliGenHIJINGpara(21507); + genHijPara->SetMomentumRange(0, 999999.); + genHijPara->SetPhiRange(-180., 180.); + // Set pseudorapidity range from -8 to 8. + thmin = EtaToTheta(8); // theta min. <---> eta max + thmax = EtaToTheta(-8); // theta max. <---> eta min + genHijPara->SetThetaRange(thmin,thmax); + fGenerator = (AliGenerator*)genHijPara; + break; + + case kParam_8000_Ecal: + genHijParaB = new AliGenHIJINGparaBa(82534); + genHijParaB->SetMomentumRange(0, 999999.); + genHijParaB->SetPhiRange(-180., 180.); + // Set pseudorapidity range from -8 to 8. + thmin = EtaToTheta( 5); // theta min. <---> eta max + thmax = EtaToTheta(-5); // theta max. <---> eta min + genHijParaB->SetThetaRange(thmin,thmax); + fGenerator = (AliGenerator*)genHijParaB; + break; + + case kParam_4000_Ecal: + genHijParaB = new AliGenHIJINGparaBa(82534/2); + genHijParaB->SetMomentumRange(0, 999999.); + genHijParaB->SetPhiRange(-180., 180.); + // Set pseudorapidity range from -8 to 8. + thmin = EtaToTheta( 5); // theta min. <---> eta max + thmax = EtaToTheta(-5); // theta max. <---> eta min + genHijParaB->SetThetaRange(thmin,thmax); + fGenerator = (AliGenerator*)genHijParaB; + break; +// +// Hijing Central +// + case kHijing_cent1: + genHij = HijingStandard(); +// impact parameter range + genHij->SetImpactParameterRange(0., 5.); + fGenerator = (AliGenerator*)genHij; + break; + case kHijing_cent2: + genHij = HijingStandard(); +// impact parameter range + genHij->SetImpactParameterRange(0., 2.); + fGenerator = (AliGenerator*)genHij; + break; +// +// Hijing Peripheral +// + case kHijing_per1: + genHij = HijingStandard(); +// impact parameter range + genHij->SetImpactParameterRange(5., 8.6); + fGenerator = (AliGenerator*)genHij; + break; + case kHijing_per2: + //coment= comment.Append("HIJING per2"); + genHij = HijingStandard(); +// impact parameter range + genHij->SetImpactParameterRange(8.6, 11.2); + fGenerator = (AliGenerator*)genHij; + break; + case kHijing_per3: + genHij = HijingStandard(); +// impact parameter range + genHij->SetImpactParameterRange(11.2, 13.2); + fGenerator = (AliGenerator*)genHij; + break; + case kHijing_per4: + genHij = HijingStandard(); +// impact parameter range + genHij->SetImpactParameterRange(13.2, 15.); + fGenerator = (AliGenerator*)genHij; + break; + case kHijing_per5: + //coment= comment.Append("HIJING per5"); + genHij = HijingStandard(); +// impact parameter range + genHij->SetImpactParameterRange(15., 100.); + fGenerator = (AliGenerator*)genHij; + break; +// +// Jet-Jet +// + case kHijing_jj25: + //coment= comment.Append("HIJING Jet 25 GeV"); + genHij = HijingStandard(); +// impact parameter range + genHij->SetImpactParameterRange(0., 5.); + // trigger + genHij->SetTrigger(1); + genHij->SetPtJet(25.); + genHij->SetSimpleJets(1); + genHij->SetRadiation(isw); + genHij->SetJetEtaRange(-0.3,0.3); + genHij->SetJetPhiRange(15.,105.); + fGenerator = (AliGenerator*)genHij; + break; + + case kHijing_jj50: + //coment= comment.Append("HIJING Jet 50 GeV"); + genHij = HijingStandard(); +// impact parameter range + genHij->SetImpactParameterRange(0., 5.); + // trigger + genHij->SetTrigger(1); + genHij->SetPtJet(50.); + genHij->SetSimpleJets(1); + genHij->SetRadiation(isw); + genHij->SetJetEtaRange(-0.3,0.3); + genHij->SetJetPhiRange(15.,105.); + fGenerator = (AliGenerator*)genHij; + break; + + case kHijing_jj75: + //coment= comment.Append("HIJING Jet 75 GeV"); + genHij = HijingStandard(); +// impact parameter range + genHij->SetImpactParameterRange(0., 5.); + // trigger + genHij->SetTrigger(1); + genHij->SetPtJet(75.); + genHij->SetSimpleJets(1); + genHij->SetRadiation(isw); + genHij->SetJetEtaRange(-0.3,0.3); + genHij->SetJetPhiRange(15.,105.); + fGenerator = (AliGenerator*)genHij; + break; + + case kHijing_jj100: + //coment= comment.Append("HIJING Jet 100 GeV"); + genHij = HijingStandard(); +// impact parameter range + genHij->SetImpactParameterRange(0., 5.); + // trigger + genHij->SetTrigger(1); + genHij->SetPtJet(100.); + genHij->SetSimpleJets(1); + genHij->SetRadiation(isw); + genHij->SetJetEtaRange(-0.3,0.3); + genHij->SetJetPhiRange(15.,105.); + fGenerator = (AliGenerator*)genHij; + break; + + case kHijing_jj125: + //coment= comment.Append("HIJING Jet 125 GeV"); + genHij = HijingStandard(); +// impact parameter range + genHij->SetImpactParameterRange(0., 5.); + // trigger + genHij->SetTrigger(1); + genHij->SetPtJet(125.); + genHij->SetSimpleJets(1); + genHij->SetRadiation(isw); + genHij->SetJetEtaRange(-0.3,0.3); + genHij->SetJetPhiRange(15.,105.); + fGenerator = (AliGenerator*)genHij; + break; +// +// Gamma-Jet +// + case kHijing_gj25: + //coment= comment.Append("HIJING Gamma 25 GeV"); + genHij = HijingStandard(); +// impact parameter range + genHij->SetImpactParameterRange(0., 5.); + // trigger + genHij->SetTrigger(2); + genHij->SetPtJet(25.); + genHij->SetSimpleJets(1); + genHij->SetRadiation(isw); + genHij->SetJetEtaRange(-0.3,0.3); + genHij->SetJetPhiRange(15.,105.); + fGenerator = (AliGenerator*)genHij; + break; + + case kHijing_gj50: + //coment= comment.Append("HIJING Gamma 50 GeV"); + genHij = HijingStandard(); +// impact parameter range + genHij->SetImpactParameterRange(0., 5.); + // trigger + genHij->SetTrigger(2); + genHij->SetPtJet(50.); + genHij->SetSimpleJets(1); + genHij->SetRadiation(isw); + genHij->SetJetEtaRange(-0.3,0.3); + genHij->SetJetPhiRange(15.,105.); + fGenerator = (AliGenerator*)genHij; + break; + + case kHijing_gj75: + //coment= comment.Append("HIJING Gamma 75 GeV"); + genHij = HijingStandard(); +// impact parameter range + genHij->SetImpactParameterRange(0., 5.); + // trigger + genHij->SetTrigger(2); + genHij->SetPtJet(75.); + genHij->SetSimpleJets(1); + genHij->SetRadiation(isw); + genHij->SetJetEtaRange(-0.3,0.3); + genHij->SetJetPhiRange(15.,105.); + fGenerator = (AliGenerator*)genHij; + break; + + case kHijing_gj100: + //coment= comment.Append("HIJING Gamma 100 GeV"); + genHij = HijingStandard(); +// impact parameter range + genHij->SetImpactParameterRange(0., 5.); + // trigger + genHij->SetTrigger(2); + genHij->SetPtJet(100.); + genHij->SetSimpleJets(1); + genHij->SetRadiation(isw); + genHij->SetJetEtaRange(-0.3,0.3); + genHij->SetJetPhiRange(15.,105.); + fGenerator = (AliGenerator*)genHij; + break; + + case kHijing_gj125: + //coment= comment.Append("HIJING Gamma 125 GeV"); + genHij = HijingStandard(); +// impact parameter range + genHij->SetImpactParameterRange(0., 5.); + // trigger + genHij->SetTrigger(2); + genHij->SetPtJet(125.); + genHij->SetSimpleJets(1); + genHij->SetRadiation(isw); + genHij->SetJetEtaRange(-0.3,0.3); + genHij->SetJetPhiRange(15.,105.); + fGenerator = (AliGenerator*)genHij; + break; + case kJetPlusBg: + genCoct = new AliGenCocktail(); + genCoct->SetMomentumRange(0, 999999.); + genCoct->SetPhiRange(-180., 180.); + // Set pseudorapidity range from -8 to 8. + thmin = EtaToTheta( 5.); // theta min. <---> eta max + thmax = EtaToTheta(-5.); // theta max. <---> eta min + genCoct->SetThetaRange(thmin,thmax); + +// +// Underlying Event +// +// AliGenHIJINGparaBa *bg = new AliGenHIJINGparaBa(82534); + bg = new AliGenHIJINGparaBa(10); + fBgGenerator = (AliGenerator*)bg; +// +// Jets from Pythia +// + jets = new AliGenPythia(-1); + fSignalGenerator = (AliGenerator*)jets; +// Centre of mass energy + jets->SetEnergyCMS(5500.); +// Process type + jets->SetProcess(kPyJets); +// final state kinematic cuts + jets->SetJetEtaRange(-0.3, 0.3); + jets->SetJetPhiRange(15., 105.); +// Structure function + jets->SetStrucFunc(kGRVLO98); +// +// Pt transfer of the hard scattering + jets->SetPtHard(100.,100.1); +// Decay type (semielectronic, semimuonic, nodecay) + jets->SetForceDecay(kAll); +// +// Add all to cockail ... +// + genCoct->AddGenerator(jets,"Jets",1); + genCoct->AddGenerator(bg,"Underlying Event", 1); + fGenerator = (AliGenerator*)genCoct; + + break; + case kGammaPlusBg: + genCoct = new AliGenCocktail(); + genCoct->SetMomentumRange(0, 999999.); + genCoct->SetPhiRange(-180., 180.); + // Set pseudorapidity range from -8 to 8. + thmin = EtaToTheta( 5.); // theta min. <---> eta max + thmax = EtaToTheta(-5.); // theta max. <---> eta min + genCoct->SetThetaRange(thmin,thmax); +// +// Underlying Event +// + bg = new AliGenHIJINGparaBa(82534); + fBgGenerator = (AliGenerator*)bg; +// +// Jets from Pythia +// + jets = new AliGenPythia(-1); + fSignalGenerator = (AliGenerator*)jets; +// Centre of mass energy + jets->SetEnergyCMS(5500.); +// Process type + jets->SetProcess(kPyDirectGamma); +// final state kinematic cuts + jets->SetJetEtaRange(-0.3, 0.3); + jets->SetJetPhiRange(15., 105.); + jets->SetGammaEtaRange(-0.12, 0.12); + jets->SetGammaPhiRange(220., 320.); +// Structure function + jets->SetStrucFunc(kGRVLO98); +// +// Pt transfer of the hard scattering + jets->SetPtHard(100.,100.1); +// Decay type (semielectronic, semimuonic, nodecay) + jets->SetForceDecay(kAll); +// +// Add all to cockail ... +// + genCoct->AddGenerator(jets,"Jets",1); + genCoct->AddGenerator(bg,"Underlying Event", 1); + fGenerator = (AliGenerator*)genCoct; + + break; + case kJets_50: +// 50 GeV Jets + genPy = PythiaJets(50.); + fGenerator = (AliGenerator*)genPy; + break; + case kJets_75: +// 75 GeV Jets + genPy = PythiaJets(75.); + fGenerator = (AliGenerator*)genPy; + break; + case kJets_100: +// 100 GeV Jets + genPy = PythiaJets(100.); + fGenerator = (AliGenerator*)genPy; + break; + case kJets_200: +// 200 GeV Jets + genPy = PythiaJets(200.); + fGenerator = (AliGenerator*)genPy; + break; + + case kJets_100RadOn: +// 100 GeV Jets with radiation on - 22-mar-2002 +// See AliPythia.cxx for default + genPy = PythiaJets(100.); + // genPy->SetKeyPartonJets(1); // for jet partons + // genPy->DefineParametersForPartonsJets(); + + fGenerator = (AliGenerator*)genPy; + break; + + case kGammaJets_50: +// 50 GeV Jets + Gamma + genPy = PythiaJets(-1); + genPy->SetEnergyCMS(5500.); + genPy->SetProcess(kPyDirectGamma); + genPy->SetJetEtaRange(-0.3,+0.3); + genPy->SetJetPhiRange(15.,105.); + genPy->SetGammaEtaRange(-0.12, 0.12); + genPy->SetGammaPhiRange(220., 320.); + genPy->SetStrucFunc(kGRVLO98); + genPy->SetPtHard(50.,50.001); + genPy->SetForceDecay(kAll); + fGenerator = (AliGenerator*)genPy; + break; + case kGammaJets_75: +// 75 GeV Jets + Gamma + genPy = PythiaJets(-1); + genPy->SetEnergyCMS(5500.); + genPy->SetProcess(kPyDirectGamma); + genPy->SetJetEtaRange(-0.3,+0.3); + genPy->SetJetPhiRange(15.,105.); + genPy->SetGammaEtaRange(-0.12, 0.12); + genPy->SetGammaPhiRange(220., 320.); + genPy->SetStrucFunc(kGRVLO98); + genPy->SetPtHard(75.,75.001); + genPy->SetForceDecay(kAll); + fGenerator = (AliGenerator*)genPy; + break; + case kGammaJets_100: +// 100 GeV Jets + Gamma + genPy = PythiaJets(-1); + genPy->SetEnergyCMS(5500.); + genPy->SetProcess(kPyDirectGamma); + genPy->SetJetEtaRange(-0.3,+0.3); + genPy->SetJetPhiRange(15.,105.); + genPy->SetGammaEtaRange(-0.12, 0.12); + genPy->SetGammaPhiRange(220., 320.); + genPy->SetStrucFunc(kGRVLO98); + genPy->SetPtHard(100.,100.001); + genPy->SetForceDecay(kAll); + fGenerator = (AliGenerator*)genPy; + break; + case kGammaJets_200: +// 200 GeV Jets + Gamma + genPy = PythiaJets(-1); + genPy->SetEnergyCMS(5500.); + genPy->SetProcess(kPyDirectGamma); + genPy->SetJetEtaRange(-0.3,+0.3); + genPy->SetJetPhiRange(15.,105.); + genPy->SetGammaEtaRange(-0.12, 0.12); + genPy->SetGammaPhiRange(220., 320.); + genPy->SetStrucFunc(kGRVLO98); + genPy->SetPtHard(200.,200.001); + genPy->SetForceDecay(kAll); + fGenerator = (AliGenPythia*)genPy; + break; + case kGammaJets_250: +// 250 GeV Jets + Gamma + genPy = PythiaJets(-1); + genPy->SetEnergyCMS(5500.); + genPy->SetProcess(kPyDirectGamma); + genPy->SetJetEtaRange(-0.3,+0.3); + genPy->SetJetPhiRange(15.,105.); + genPy->SetGammaEtaRange(-0.12, 0.12); + genPy->SetGammaPhiRange(220., 320.); + genPy->SetStrucFunc(kGRVLO98); + genPy->SetPtHard(250.,250.001); + genPy->SetForceDecay(kAll); + fGenerator = (AliGenerator*)genPy; + break; + case kGammaJets_300: +// 300 GeV Jets + Gamma + genPy = PythiaJets(-1); + genPy->SetEnergyCMS(5500.); + genPy->SetProcess(kPyDirectGamma); + genPy->SetJetEtaRange(-0.3,+0.3); + genPy->SetJetPhiRange(15.,105.); + genPy->SetGammaEtaRange(-0.12, 0.12); + genPy->SetGammaPhiRange(220., 320.); + genPy->SetStrucFunc(kGRVLO98); + genPy->SetPtHard(300.,300.001); + genPy->SetForceDecay(kAll); + fGenerator = (AliGenerator*)genPy; + break; + default: + printf(" 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(" 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(" AliEMCALGeneratorFactory::OneParticleWithFixedEnergy \n"); + printf(" type of particle -> %i \n", type); + printf(" %6.4f < eta < %6.4f\n", thmin, thmax); + printf(" %6.4f < phi < %6.4f\n", phimin, phimax); + printf(" %7.2f < Momentum < %7.2f\n", pmin, pmax); + printf(" TestBit(kPtRange) %i | TestBit(kMomentumRange) %i\n", + gen->TestBit(BIT(17)), gen->TestBit(BIT(19))); + return gen; +} diff --git a/EMCAL/AliEMCALGeneratorFactory.h b/EMCAL/AliEMCALGeneratorFactory.h new file mode 100644 index 00000000000..252882dcb33 --- /dev/null +++ b/EMCAL/AliEMCALGeneratorFactory.h @@ -0,0 +1,79 @@ +#ifndef ALIEMCALGENERATORFACTORY_H +#define ALIEMCALGENERATORFACTORY_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id$ */ + +//_________________________________________________________________________ +// Class for generator factory which used in production for EMCAL. +// Based on Config.C file +//*-- Author: Aleksei Pavlinov (WSU) +#include +#include +#include "AliPDG.h" +// 23-aug-04 for kGamma and kPi0 +#include + +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 diff --git a/EMCAL/AliEMCALGeometry.cxx b/EMCAL/AliEMCALGeometry.cxx index f6f08169e3b..a641e9437bd 100644 --- a/EMCAL/AliEMCALGeometry.cxx +++ b/EMCAL/AliEMCALGeometry.cxx @@ -26,6 +26,7 @@ //*-- Author: Sahal Yacoob (LBL / UCT) // and : Yves Schutz (SUBATECH) // and : Jennifer Klay (LBL) +// SHASHLYK : Aleksei Pavlinov (WSU) // --- AliRoot header files --- #include @@ -41,6 +42,7 @@ ClassImp(AliEMCALGeometry) AliEMCALGeometry *AliEMCALGeometry::fgGeom = 0; Bool_t AliEMCALGeometry::fgInit = kFALSE; +TString name; // contains name of geometry //______________________________________________________________________ AliEMCALGeometry::~AliEMCALGeometry(void){ @@ -66,9 +68,21 @@ void AliEMCALGeometry::Init(void){ // WX inform about the composition of the EM calorimeter section: // thickness in mm of Pb radiator (W) and of scintillator (X), and number of scintillator layers (N) // New geometry: EMCAL_55_25 - + // 24-aug-04 for shish-kebab + // SHISH_25 or SHISH_62 fgInit = kFALSE; // Assume failed until proven otherwise. - TString name(GetName()) ; + name = GetName(); + name.ToUpper(); + + fNZ = 114; // granularity along Z (eta) + fNPhi = 168; // granularity in phi (azimuth) + fArm1PhiMin = 60.0; // degrees, Starting EMCAL Phi position + fArm1PhiMax = 180.0; // degrees, Ending EMCAL Phi position + fArm1EtaMin = -0.7; // pseudorapidity, Starting EMCAL Eta position + fArm1EtaMax = +0.7; // pseudorapidity, Ending EMCAL Eta position + fIPDistance = 454.0; // cm, Radial distance to inner surface of EMCAL + + // geometry if (name == "EMCAL_55_25") { fECPbRadThickness = 0.5; // cm, Thickness of the Pb radiators fECScintThick = 0.5; // cm, Thickness of the scintillator @@ -82,21 +96,114 @@ void AliEMCALGeometry::Init(void){ else if( name == "G56_2_55_19" || name == "EMCAL_5655_21" || name == "G56_2_55_19_104_14"|| name == "G65_2_64_19" || name == "EMCAL_6564_21"){ Fatal("Init", "%s is an old geometry! Please update your Config file", name.Data()) ; } + else if(name.Contains("SHISH")){ + fNumberOfSuperModules = 12; // 12 = 6 * 2 (6 in phi, 2 in Z); + fSteelFrontThick = 2.54; // 9-sep-04 + fIPDistance = 460.0; + fFrontSteelStrip = fPassiveScintThick = 0.0; // 13-may-05 + fLateralSteelStrip = 0.025; // before MAY 2005 + fPhiModuleSize = fEtaModuleSize = 11.4; + fPhiTileSize = fEtaTileSize = 5.52; // (11.4-5.52*2)/2. = 0.18 cm (wall thickness) + fNPhi = 14; + fNZ = 30; + fAlFrontThick = fGap2Active = 0; + fNPHIdiv = fNETAdiv = 2; + + fNECLayers = 62; + fECScintThick = fECPbRadThickness = 0.2; + fSampling = 1.; // 30-aug-04 - should be calculated + if(name.Contains("TWIST")) { // all about EMCAL module + fNZ = 27; // 16-sep-04 + } else if(name.Contains("TRD")) { + fIPDistance = 428.0; // 11-may-05 + fSteelFrontThick = 0.0; // 3.17 -> 0.0; 28-mar-05 : no stell plate + fNPhi = 12; + fSampling = 12.327; + fPhiModuleSize = fEtaModuleSize = 12.26; + fNZ = 26; // 11-oct-04 + fTrd1Angle = 1.3; // in degree +// 18-nov-04; 1./0.08112=12.327 +// http://pdsfweb01.nersc.gov/~pavlinov/ALICE/SHISHKEBAB/RES/linearityAndResolutionForTRD1.html + if(name.Contains("TRD1")) { // 30-jan-05 + // for final design + if(name.Contains("MAY05") || name.Contains("WSUC")){ + fNumberOfSuperModules = 12; // 20-may-05 + if(name.Contains("WSUC")) fNumberOfSuperModules = 1; // 27-may-05 + fNECLayers = 77; // (13-may-05 from V.Petrov) + fPhiModuleSize = 12.5; // 20-may-05 - rectangular shape + fEtaModuleSize = 11.9; + fECScintThick = fECPbRadThickness = 0.16;// (13-may-05 from V.Petrov) + fFrontSteelStrip = 0.025;// 0.025cm = 0.25mm (13-may-05 from V.Petrov) + fLateralSteelStrip = 0.01; // 0.01cm = 0.1mm (13-may-05 from V.Petrov) - was 0.025 + fPassiveScintThick = 0.8; // 0.8cm = 8mm (13-may-05 from V.Petrov) + fNZ = 24; + fTrd1Angle = 1.5; // 1.3 or 1.5 + } + } else if(name.Contains("TRD2")) { // 30-jan-05 + fSteelFrontThick = 0.0; // 11-mar-05 + fIPDistance+= fSteelFrontThick; // 1-feb-05 - compensate absence of steel plate + fTrd1Angle = 1.64; // 1.3->1.64 + fTrd2AngleY = fTrd1Angle; // symmetric case now + fEmptySpace = 0.2; // 2 mm + fTubsR = fIPDistance; // 31-jan-05 - as for Fred case + + fPhiModuleSize = fTubsR*2.*TMath::Tan(fTrd2AngleY*TMath::DegToRad()/2.); + fPhiModuleSize -= fEmptySpace/2.; // 11-mar-05 + fEtaModuleSize = fPhiModuleSize; // 20-may-05 + fTubsTurnAngle = 3.; + } + fNPHIdiv = fNETAdiv = 2; // 13-oct-04 - division again + if(name.Contains("3X3")) { // 23-nov-04 + fNPHIdiv = fNETAdiv = 3; + } else if(name.Contains("4X4")) { + fNPHIdiv = fNETAdiv = 4; + } + } + fPhiTileSize = fPhiModuleSize/2. - fLateralSteelStrip; // 13-may-05 + fEtaTileSize = fEtaModuleSize/2. - fLateralSteelStrip; // 13-may-05 + + if(name.Contains("25")){ + fNECLayers = 25; + fECScintThick = fECPbRadThickness = 0.5; + } + if(name.Contains("WSUC")){ // 18-may-05 - about common structure + fShellThickness = 30.; // should be change + fNPhi = fNZ = 4; + } + // constant for transition absid <--> indexes + fNCellsInTower = fNPHIdiv*fNETAdiv; + fNCellsInSupMod = fNCellsInTower*fNPhi*fNZ; + fNCells = fNCellsInSupMod*fNumberOfSuperModules; + + fLongModuleSize = fNECLayers*(fECScintThick + fECPbRadThickness); + if(name.Contains("MAY05")) fLongModuleSize += (fFrontSteelStrip + fPassiveScintThick); + + // 30-sep-04 + if(name.Contains("TRD")) { + f2Trd1Dx2 = fEtaModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd1Angle*TMath::DegToRad()/2.); + if(name.Contains("TRD2")) { // 27-jan-05 + f2Trd2Dy2 = fPhiModuleSize + 2.*fLongModuleSize*TMath::Tan(fTrd2AngleY*TMath::DegToRad()/2.); + } + } + } else Fatal("Init", "%s is an undefined geometry!", name.Data()) ; - // geometry - fNZ = 114; // granularity along Z (eta) - fNPhi = 168; // granularity in phi (azimuth) - fArm1PhiMin = 60.0; // degrees, Starting EMCAL Phi position - fArm1PhiMax = 180.0; // degrees, Ending EMCAL Phi position - fArm1EtaMin = -0.7; // pseudorapidity, Starting EMCAL Eta position - fArm1EtaMax = +0.7; // pseudorapidity, Ending EMCAL Eta position - - fIPDistance = 454.0; // cm, Radial distance to inner surface of EMCAL + fNPhiSuperModule = fNumberOfSuperModules/2; + if(fNPhiSuperModule<1) fNPhiSuperModule = 1; //There is always one more scintillator than radiator layer because of the first block of aluminium fShellThickness = fAlFrontThick + fGap2Active + fNECLayers*GetECScintThick()+(fNECLayers-1)*GetECPbRadThick(); + if(name.Contains("SHISH")) { + fShellThickness = fSteelFrontThick + fLongModuleSize; + if(name.Contains("TWIST")) { // 13-sep-04 + fShellThickness = TMath::Sqrt(fLongModuleSize*fLongModuleSize + fPhiModuleSize*fEtaModuleSize); + fShellThickness += fSteelFrontThick; + } else if(name.Contains("TRD")) { // 1-oct-04 + fShellThickness = TMath::Sqrt(fLongModuleSize*fLongModuleSize + f2Trd1Dx2*f2Trd1Dx2); + fShellThickness += fSteelFrontThick; + } + } fZLength = 2.*ZFromEtaR(fIPDistance+fShellThickness,fArm1EtaMax); // Z coverage fEnvelop[0] = fIPDistance; // mother volume inner radius @@ -105,12 +212,41 @@ void AliEMCALGeometry::Init(void){ fgInit = kTRUE; - if (gDebug) { - printf("Init: geometry of EMCAL named %s is as follows:", name.Data()); + if (kTRUE) { + printf("Init: geometry of EMCAL named %s is as follows:\n", name.Data()); printf( " ECAL : %d x (%f mm Pb, %f mm Sc) \n", GetNECLayers(), GetECPbRadThick(), GetECScintThick() ) ; + if(name.Contains("SHISH")){ + printf(" fIPDistance %6.3f cm \n", fIPDistance); + if(fSteelFrontThick>0.) + printf(" fSteelFrontThick %6.3f cm \n", fSteelFrontThick); + printf(" fNPhi %i | fNZ %i \n", fNPhi, fNZ); + if(name.Contains("MAY05")){ + printf(" fFrontSteelStrip %6.4f cm (thickness of front steel strip)\n", + fFrontSteelStrip); + printf(" fLateralSteelStrip %6.4f cm (thickness of lateral steel strip)\n", + fLateralSteelStrip); + printf(" fPassiveScintThick %6.4f cm (thickness of front passive Sc tile)\n", + fPassiveScintThick); + } + printf(" X:Y module size %6.3f , %6.3f cm \n", fPhiModuleSize, fEtaModuleSize); + printf(" X:Y tile size %6.3f , %6.3f cm \n", fPhiTileSize, fEtaTileSize); + printf(" fLongModuleSize %6.3f cm \n", fLongModuleSize); + printf(" #supermodule in phi direction %i \n", fNPhiSuperModule ); + } + if(name.Contains("TRD")) { + printf(" fTrd1Angle %7.4f\n", fTrd1Angle); + printf(" f2Trd1Dx2 %7.4f\n", f2Trd1Dx2); + if(name.Contains("TRD2")) { + printf(" fTrd2AngleY %7.4f\n", fTrd2AngleY); + printf(" f2Trd2Dy2 %7.4f\n", f2Trd2Dy2); + printf(" fTubsR %7.2f\n", fTubsR); + printf(" fTubsTurnAngle %7.4f\n", fTubsTurnAngle); + printf(" fEmptySpace %7.4f\n", fEmptySpace); + } + } printf("Granularity: %d in eta and %d in phi\n", GetNZ(), GetNPhi()) ; - printf("Layout: phi = (%f, %f), eta = (%f, %f), y = %f\n", - GetArm1PhiMin(), GetArm1PhiMax(),GetArm1EtaMin(), GetArm1EtaMax(), GetIPDistance() ) ; + printf("Layout: phi = (%7.1f, %7.1f), eta = (%5.2f, %5.2f), IP = %7.2f\n", + GetArm1PhiMin(), GetArm1PhiMax(),GetArm1EtaMin(), GetArm1EtaMax(), GetIPDistance() ); } } @@ -431,3 +567,66 @@ Bool_t AliEMCALGeometry::IsInEMCAL(Double_t x, Double_t y, Double_t z) const { } return 0; } + +// +// == Shish-kebab cases == +// +Int_t AliEMCALGeometry::GetAbsCellId(const int nSupMod, const int nTower, const int nIphi, const int nIeta) +{ // 27-aug-04; corr. 21-sep-04 + static Int_t id; // have to change from 1 to fNCells + id = fNCellsInSupMod*(nSupMod-1); + id += fNCellsInTower *(nTower-1); + id += fNPHIdiv *(nIphi-1); + id += nIeta; + if(id<=0 || id > fNCells) { + printf(" wrong numerations !!\n"); + printf(" id %6i(will be force to -1)\n", id); + printf(" fNCells %6i\n", fNCells); + printf(" nSupMod %6i\n", nSupMod); + printf(" nTower %6i\n", nTower); + printf(" nIphi %6i\n", nIphi); + printf(" nIeta %6i\n", nIeta); + id = -1; + } + return id; +} + +Bool_t AliEMCALGeometry::CheckAbsCellId(Int_t ind) +{ // 17-niv-04 - analog of IsInECA + if(name.Contains("TRD")) { + if(ind<=0 || ind > fNCells) return kFALSE; + else return kTRUE; + } else return IsInECA(ind); +} + +Bool_t AliEMCALGeometry::GetCellIndex(const Int_t absId,Int_t &nSupMod,Int_t &nTower,Int_t &nIphi,Int_t &nIeta) +{ // 21-sep-04 + static Int_t tmp=0; + if(absId<=0 || absId>fNCells) { + Info("GetCellIndex"," wrong abs Id %i !! \n", absId); + return kFALSE; + } + nSupMod = (absId-1) / fNCellsInSupMod + 1; + tmp = (absId-1) % fNCellsInSupMod; + + nTower = tmp / fNCellsInTower + 1; + tmp = tmp % fNCellsInTower; + + nIphi = tmp / fNPHIdiv + 1; + nIeta = tmp % fNPHIdiv + 1; + + return kTRUE; +} + +void AliEMCALGeometry::GetCellPhiEtaIndexInSModule(const int nTower, const int nIphi, const int nIeta, +int &iphi, int &ieta) +{ // don't check validity of nTower, nIphi and nIeta index + // have to change - 1-nov-04 ?? + static Int_t iphit, ietat; + + ietat = (nTower-1)/fNPhi; + ieta = ietat*fNETAdiv + nIeta; // change from 1 to fNZ*fNETAdiv + + iphit = (nTower-1)%fNPhi; + iphi = iphit*fNPHIdiv + nIphi; // change from 1 to fNPhi*fNPHIdiv +} diff --git a/EMCAL/AliEMCALGeometry.h b/EMCAL/AliEMCALGeometry.h index 3dd0be00586..cb17c68d88e 100644 --- a/EMCAL/AliEMCALGeometry.h +++ b/EMCAL/AliEMCALGeometry.h @@ -1,6 +1,6 @@ #ifndef ALIEMCALGEOMETRY_H #define ALIEMCALGEOMETRY_H -/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * +/* Copyright(c) 1998-2004, ALICE Experiment at CERN, All rights reserved. * * See cxx source for full Copyright notice */ /* $Id$ */ @@ -78,7 +78,34 @@ public: Float_t GetECScintThick() const {return fECScintThick;} Float_t GetSampling() const {return fSampling ; } Bool_t IsInECA(Int_t index) const { if ( (index > 0 && (index <= GetNZ() * GetNPhi()))) return kTRUE; else return kFALSE ;} - + + Int_t GetNumberOfSuperModules() {return fNumberOfSuperModules;} + Float_t GetPhiModuleSize() const {return fPhiModuleSize;} + Float_t GetEtaModuleSize() const {return fEtaModuleSize;} + Float_t GetFrontSteelStrip() const {return fFrontSteelStrip;} + Float_t GetLateralSteelStrip() const {return fLateralSteelStrip;} + Float_t GetPassiveScintThick() const {return fPassiveScintThick;} + Float_t GetPhiTileSize() const {return fPhiTileSize;} + Float_t GetEtaTileSize() const {return fEtaTileSize;} + Int_t GetNPhiSuperModule() const {return fNPhiSuperModule;} + Int_t GetNPHIdiv() const {return fNPHIdiv ;} + Int_t GetNETAdiv() const {return fNETAdiv ;} + Int_t GetNCells() const {return fNCells;} + Float_t GetSteelFrontThickness() const { return fSteelFrontThick;} + Float_t GetLongModuleSize() const {return fLongModuleSize;} + + Float_t GetTrd1Angle() const {return fTrd1Angle;} + Float_t Get2Trd1Dx2() const {return f2Trd1Dx2;} + Float_t GetTrd2AngleY()const {return fTrd2AngleY;} + Float_t Get2Trd2Dy2() const {return f2Trd2Dy2;} + Float_t GetTubsR() const {return fTubsR;} + Float_t GetTubsTurnAngle() const {return fTubsTurnAngle;} + // Dabs id <-> indexes; Shish-kebab case + Int_t GetAbsCellId(const Int_t nSupMod, const Int_t nTower, const Int_t nIphi, const Int_t nIeta); + Bool_t GetCellIndex(const Int_t absId, Int_t &nSupMod, Int_t &nTower, Int_t &nIphi, Int_t &nIeta); + void GetCellPhiEtaIndexInSModule(const Int_t nTower, const Int_t nIphi, const Int_t nIeta, Int_t &iphi, Int_t &ieta); + Bool_t CheckAbsCellId(Int_t ind); // replace the IsInECA + // --- Float_t AngleFromEta(Float_t eta){ // returns theta in radians for a given pseudorapidity return 2.0*TMath::ATan(TMath::Exp(-eta)); } @@ -132,11 +159,39 @@ private: Float_t fZLength; // Total length in z direction Float_t fGap2Active; // Gap between the envelop and the active material Int_t fNZ; // Number of Towers in the Z direction - Int_t fNPhi; // Number of Towers in the Phi Direction + Int_t fNPhi; // Number of Towers in the PHI direction Float_t fSampling; // Sampling factor - - ClassDef(AliEMCALGeometry,8) // EMCAL geometry class - + + // Shish-kebab option - 23-aug-04 by PAI; COMPACT, TWIST, TRD1 and TRD2 + Int_t fNumberOfSuperModules; // default is 12 = 6 * 2 + Float_t fSteelFrontThick; // Thickness of the front stell face of the support box - 9-sep-04 + Float_t fFrontSteelStrip; // 13-may-05 + Float_t fLateralSteelStrip; // 13-may-05 + Float_t fPassiveScintThick; // 13-may-05 + Float_t fPhiModuleSize; // Phi -> X + Float_t fEtaModuleSize; // Eta -> Y + Float_t fPhiTileSize; // + Float_t fEtaTileSize; // + Float_t fLongModuleSize; // + Int_t fNPhiSuperModule; // 6 - number supermodule in phi direction + Int_t fNPHIdiv; // number phi dvizion + Int_t fNETAdiv; // number eta divizion + // + Int_t fNCells; // number of cells in calo + Int_t fNCellsInSupMod; // number cell in super module + Int_t fNCellsInTower; // number cell in tower + // TRD1 options - 30-sep-04 + Float_t fTrd1Angle; // angle in x-z plane (in degree) + Float_t f2Trd1Dx2; // 2*dx2 for TRD1 + // TRD2 options - 27-jan-07 + Float_t fTrd2AngleY; // angle in y-z plane (in degree) + Float_t f2Trd2Dy2; // 2*dy2 for TRD2 + Float_t fEmptySpace; // 2mm om fred drawing + // Sumper module as TUBS + Float_t fTubsR; // radius of tubs + Float_t fTubsTurnAngle; // turn angle of tubs in degree + + ClassDef(AliEMCALGeometry,9) // EMCAL geometry class }; #endif // AliEMCALGEOMETRY_H diff --git a/EMCAL/AliEMCALGeometryOfflineTrd1.cxx b/EMCAL/AliEMCALGeometryOfflineTrd1.cxx new file mode 100644 index 00000000000..82706c41322 --- /dev/null +++ b/EMCAL/AliEMCALGeometryOfflineTrd1.cxx @@ -0,0 +1,197 @@ +/************************************************************************** + * Copyright(c) 1998-2004, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ +// GeometryOfflineTrd1 class for EMCAL : singleton +//*-- Author: Aleksei Pavlinov (WSU) + +/* $Id$*/ + +#include +#include +#include +#include +#include + +#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; iGetCenterOfCell(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; iGetNPhi(); 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; iGetArm1PhiMin() + 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; iGetCellIndex(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; i0) b->Add(fTrd1Modules[i]); + + for(Int_t i=0; iAdd(&fSuperModuleRotationZ[i], fNameSuperModuleRotationZ[i].Data()); + for(Int_t j=0; j<2; j++) { + TString name("rotationM_"); name += (2*i+j); + b->Add(&fSuperModuleRotation[2*i+j], name.Data()); + } + } + + b->Add(&fSuperModuleRotationX, "rotationX_180"); + + b->Add(fXYZofCells); +} + +Bool_t AliEMCALGeometryOfflineTrd1::IsFolder() const +{ + return kTRUE; +} diff --git a/EMCAL/AliEMCALGeometryOfflineTrd1.h b/EMCAL/AliEMCALGeometryOfflineTrd1.h new file mode 100644 index 00000000000..5b83ba3fcfa --- /dev/null +++ b/EMCAL/AliEMCALGeometryOfflineTrd1.h @@ -0,0 +1,78 @@ +#ifndef ALIEMCALGEOMETRYOFFLINETRD1_H +#define ALIEMCALGEOMETRYOFFLINETRD1_H +/* Copyright(c) 1998-2004, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id$ */ + +//_________________________________________________________________________ +// GeometryOfflineTrd1 class for EMCAL : singleton +//*-- Author: Aleksei Pavlinov (WSU) + +#include +#include +#include + +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 diff --git a/EMCAL/AliEMCALJetMicroDst.cxx b/EMCAL/AliEMCALJetMicroDst.cxx index 79328b60282..496cc9f6b10 100644 --- a/EMCAL/AliEMCALJetMicroDst.cxx +++ b/EMCAL/AliEMCALJetMicroDst.cxx @@ -759,3 +759,61 @@ TList* AliEMCALJetMicroDst::MoveHistsToList(const char* name, Bool_t putToBrowse if(putToBrowser) gROOT->GetListOfBrowsables()->Add((TObject*)list); return list; } + +void AliEMCALJetMicroDst::FillH1(TList *l, Int_t ind, Double_t x, Double_t w) +{ + static TH1* hid=0; + if(l == 0) return; + if(ind < l->GetSize()){ + hid = (TH1*)l->At(ind); + hid->Fill(x,w); + } +} + +void AliEMCALJetMicroDst::FillH2(TList *l, Int_t ind, Double_t x, Double_t y, Double_t w) +{ + static TH2* hid=0; + if(l == 0) return; + if(ind < l->GetSize()){ + hid = (TH2*)l->At(ind); + hid->Fill(x,y,w); + } +} + +int AliEMCALJetMicroDst::SaveListOfHists(TList *list,const char* name,Bool_t kSingleKey,const char* opt) +{ + printf(" Name of out file |%s|\n", name); + int save = 0; + if(list && list->GetSize() && strlen(name)){ + TString nf(name); + if(nf.Contains(".root") == kFALSE) nf += ".root"; + TFile file(nf.Data(),opt); + TIter nextHist(list); + TObject* objHist=0; + int nh=0; + if(kSingleKey) { + file.cd(); + list->Write(list->GetName(),TObject::kSingleKey); + list->ls(); + save = 1; + } else { + while((objHist=nextHist())) { // loop over list + if(objHist->InheritsFrom("TH1")) { + TH1* hid = (TH1*)objHist; + file.cd(); + hid->Write(); + nh++; + printf("Save hist. %s \n",hid ->GetName()); + } + } + printf("%i hists. save to file -> %s\n", nh,file.GetName()); + if(nh>0) save = 1; + } + file.Close(); + } else { + printf("TAliasPAI::saveListOfHists : N O S A V I N G \n"); + if(list==0) printf("List of object 0 : %i \n", (int)list); + else printf("Size of list %i \n", list->GetSize()); + } + return save; +} diff --git a/EMCAL/AliEMCALJetMicroDst.h b/EMCAL/AliEMCALJetMicroDst.h index fea92a0be03..8ddc6bff341 100644 --- a/EMCAL/AliEMCALJetMicroDst.h +++ b/EMCAL/AliEMCALJetMicroDst.h @@ -23,8 +23,6 @@ class TVector3; class TBrowser; class AliEMCALJetMicroDst: public TNamed { - - public: AliEMCALJetMicroDst(const char *name="jetMicroDst", const char *tit="jet Micro Dst for preparation of proposal"); @@ -80,7 +78,12 @@ class AliEMCALJetMicroDst: public TNamed { virtual Bool_t IsFolder() const; virtual void Browse(TBrowser* b) const ; - static TList *MoveHistsToList(const char* name="List of Hist", Bool_t putToBrowser=kTRUE); + // service routine + static TList *MoveHistsToList(const char* name="ListOfHists", Bool_t putToBrowser=kTRUE); + static void FillH1(TList *l=0, Int_t ind=0, Double_t x=-99999., Double_t w=1.); + static void FillH2(TList *l=0, Int_t ind=0, Double_t x=-99999., Double_t w=-99999., Double_t w=1.); + static int SaveListOfHists(TList *list=0, const char* name="test", Bool_t kSingleKey=kFALSE, + const char* opt="RECREATE"); AliEMCALJetMicroDst & operator = (const AliEMCALJetMicroDst &) { Fatal("operator =", "not implemented") ; return *this ; } @@ -142,6 +145,9 @@ class AliEMCALJetMicroDst: public TNamed { }; #endif // AliEMCALJETMICRODST_H + +typedef AliEMCALJetMicroDst sv; // for convinience + /* What to do 1. Common info about event diff --git a/EMCAL/AliEMCALRecPoint.cxx b/EMCAL/AliEMCALRecPoint.cxx index 98791d5928e..57507367ed2 100644 --- a/EMCAL/AliEMCALRecPoint.cxx +++ b/EMCAL/AliEMCALRecPoint.cxx @@ -38,7 +38,6 @@ ClassImp(AliEMCALRecPoint) - //____________________________________________________________________________ AliEMCALRecPoint::AliEMCALRecPoint() : AliRecPoint() @@ -135,10 +134,26 @@ Bool_t AliEMCALRecPoint::AreNeighbours(AliEMCALDigit * digit1, AliEMCALDigit * d AliEMCALGeometry * geom = (AliEMCALGetter::Instance())->EMCALGeometry(); Int_t relid1[2] ; - geom->AbsToRelNumbering(digit1->GetId(), relid1) ; + //copied for shish-kebab geometry, ieta,iphi is cast as float as eta,phi conversion + // for this geometry does not exist + int nSupMod=0, nTower=0, nIphi=0, nIeta=0; + int iphi=0, ieta=0; + geom->GetCellIndex(digit1->GetId(), nSupMod,nTower,nIphi,nIeta); + geom->GetCellPhiEtaIndexInSModule(nTower,nIphi,nIeta, iphi,ieta); + relid1[0]=ieta; + relid1[1]=iphi; +// geom->AbsToRelNumbering(digit1->GetId(), relid1) ; Int_t relid2[2] ; - geom->AbsToRelNumbering(digit2->GetId(), relid2) ; + //copied for shish-kebab geometry, ieta,iphi is cast as float as eta,phi conversion + // for this geometry does not exist + int nSupMod1=0, nTower1=0, nIphi1=0, nIeta1=0; + int iphi1=0, ieta1=0; + geom->GetCellIndex(digit2->GetId(), nSupMod1,nTower1,nIphi1,nIeta1); + geom->GetCellPhiEtaIndexInSModule(nTower1,nIphi1,nIeta1, iphi1,ieta1); + relid2[0]=ieta1; + relid2[1]=iphi1; +// geom->AbsToRelNumbering(digit2->GetId(), relid2) ; Int_t rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ; Int_t coldiff = TMath::Abs( relid1[1] - relid2[1] ) ; @@ -293,13 +308,20 @@ void AliEMCALRecPoint::EvalAll(Float_t logWeight,TClonesArray * digits) // Evaluates all shower parameters EvalLocalPosition(logWeight, digits) ; + printf("eval position done\n"); EvalElipsAxis(logWeight, digits) ; + printf("eval axis done\n"); EvalDispersion(logWeight, digits) ; - EvalCoreEnergy(logWeight, digits); + printf("eval dispersion done\n"); + // EvalCoreEnergy(logWeight, digits); + // printf("eval energy done\n"); EvalTime(digits) ; + printf("eval time done\n"); EvalPrimaries(digits) ; + printf("eval pri done\n"); EvalParents(digits); + printf("eval parent done\n"); } //____________________________________________________________________________ @@ -320,7 +342,7 @@ void AliEMCALRecPoint::EvalDispersion(Float_t logWeight, TClonesArray * digits) if (!fLocPos.X() || !fLocPos.Y()) EvalLocalPosition(logWeight, digits) ; - const Float_t kDeg2Rad = TMath::DegToRad() ; + //Sub const Float_t kDeg2Rad = TMath::DegToRad() ; Float_t cluEta = fLocPos.X() ; Float_t cluPhi = fLocPos.Y() ; @@ -335,8 +357,24 @@ void AliEMCALRecPoint::EvalDispersion(Float_t logWeight, TClonesArray * digits) digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit]) ; Float_t etai = 0.; Float_t phii = 0.; - geom->EtaPhiFromIndex(digit->GetId(), etai, phii); + //copied for shish-kebab geometry, ieta,iphi is cast as float as eta,phi conversion + // for this geometry does not exist + int nSupMod=0, nTower=0, nIphi=0, nIeta=0; + int iphi=0, ieta=0; + geom->GetCellIndex(digit->GetId(), nSupMod,nTower,nIphi,nIeta); + geom->GetCellPhiEtaIndexInSModule(nTower,nIphi,nIeta, iphi,ieta); + etai=(Float_t)ieta; + phii=(Float_t)iphi; + printf("%f,%d,%d \n", fAmp, ieta, iphi) ; + +/* Sub + geom->EtaPhiFromIndex(digit->GetId(), etai, phii); phii = phii * kDeg2Rad; +*/ +/////////////////////////// + if(fAmp>0)printf("%f %d %f", fAmp,iDigit,fEnergyList[iDigit]) ; +///////////////////////// + if (gDebug == 2) printf("EvalDispersion: id = %d, etai,phii = %f,%f", digit->GetId(), etai, phii) ; @@ -351,7 +389,7 @@ void AliEMCALRecPoint::EvalDispersion(Float_t logWeight, TClonesArray * digits) d = 0. ; fDispersion = TMath::Sqrt(d) ; - + printf("Dispersion: = %f", fDispersion) ; } //____________________________________________________________________________ @@ -367,15 +405,23 @@ void AliEMCALRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digit Int_t iDigit; Float_t cluEta = 0; Float_t cluPhi = 0; - const Float_t kDeg2Rad = TMath::DegToRad(); + //Sub const Float_t kDeg2Rad = TMath::DegToRad(); for(iDigit=0; iDigit(digits->At(fDigitsList[iDigit])) ; Float_t etai ; Float_t phii ; - geom->EtaPhiFromIndex(digit->GetId(), etai, phii); - phii = phii * kDeg2Rad; + //copied for shish-kebab geometry, ieta,iphi is cast as float as eta,phi conversion + // for this geometry does not exist + int nSupMod=0, nTower=0, nIphi=0, nIeta=0; + int iphi=0, ieta=0; + geom->GetCellIndex(digit->GetId(), nSupMod,nTower,nIphi,nIeta); + geom->GetCellPhiEtaIndexInSModule(nTower,nIphi,nIeta, iphi,ieta); + etai=(Float_t)ieta; + phii=(Float_t)iphi; +//Sub geom->EtaPhiFromIndex(digit->GetId(), etai, phii); +//Sub phii = phii * kDeg2Rad; Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ; cluEta += (etai * w) ; cluPhi += (phii * w ); @@ -394,7 +440,7 @@ void AliEMCALRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digit fLocPos.SetY(cluPhi); fLocPos.SetZ(geom->GetIP2ECASection()); - if (gDebug==2) +// if (gDebug==2) printf("EvalLocalPosition: eta,phi,r = %f,%f,%f", fLocPos.X(), fLocPos.Y(), fLocPos.Z()) ; fLocPosM = 0 ; } @@ -445,6 +491,7 @@ void AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits) AliEMCALDigit * digit ; AliEMCALGeometry * geom = (AliEMCALGetter::Instance())->EMCALGeometry(); + TString gn(geom->GetName()); Int_t iDigit; @@ -452,15 +499,29 @@ void AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits) digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit]) ; Float_t etai = 0. ; Float_t phii = 0. ; - geom->EtaPhiFromIndex(digit->GetId(), etai, phii); - phii = phii * kDeg2Rad; + if(gn.Contains("SHISH")) { + //copied for shish-kebab geometry, ieta,iphi is cast as float as eta,phi conversion + // for this geometry does not exist + int nSupMod=0, nTower=0, nIphi=0, nIeta=0; + int iphi=0, ieta=0; + geom->GetCellIndex(digit->GetId(), nSupMod,nTower,nIphi,nIeta); + geom->GetCellPhiEtaIndexInSModule(nTower,nIphi,nIeta, iphi,ieta); + etai=(Float_t)ieta; + phii=(Float_t)iphi; + } else { + geom->EtaPhiFromIndex(digit->GetId(), etai, phii); + phii = phii * kDeg2Rad; + } + Double_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; dxx += w * etai * etai ; x += w * etai ; dzz += w * phii * phii ; z += w * phii ; + dxz += w * etai * phii ; + wtot += w ; } @@ -491,6 +552,7 @@ void AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits) fLambda[1]= 0. ; } + // printf("Evalaxis: lambdas = %f,%f", fLambda[0],fLambda[1]) ; } diff --git a/EMCAL/AliEMCALReconstruct.C b/EMCAL/AliEMCALReconstruct.C index b7421be83d1..d3823cce9bc 100644 --- a/EMCAL/AliEMCALReconstruct.C +++ b/EMCAL/AliEMCALReconstruct.C @@ -238,5 +238,4 @@ void EMCALHits2Digits (Bool_t split=kFALSE, TString fileName = "galice.root") { } } - diff --git a/EMCAL/AliEMCALReconstructor.cxx b/EMCAL/AliEMCALReconstructor.cxx index 204a0dcccd4..d5caaded9a3 100644 --- a/EMCAL/AliEMCALReconstructor.cxx +++ b/EMCAL/AliEMCALReconstructor.cxx @@ -31,7 +31,6 @@ #include "AliEMCALClusterizerv1.h" #include "AliEMCALPIDv1.h" #include "AliEMCALGetter.h" -#include "AliEMCALTracker.h" #include "AliRawReaderFile.h" ClassImp(AliEMCALReconstructor) @@ -100,6 +99,19 @@ void AliEMCALReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const Int_t eventNumber = runLoader->GetEventNumber() ; + TString headerFile(runLoader->GetFileName()) ; + TString branchName(runLoader->GetEventFolder()->GetName()) ; + + AliEMCALPIDv1 pid(headerFile, branchName); + + + // do current event; the loop over events is done by AliReconstruction::Run() + pid.SetEventRange(eventNumber, eventNumber) ; + if ( Debug() ) + pid.ExecuteTask("deb all") ; + else + pid.ExecuteTask("") ; + // Creates AliESDtrack from AliEMCALRecParticles AliEMCALGetter::Instance()->Event(eventNumber, "P") ; TClonesArray *recParticles = AliEMCALGetter::Instance()->RecParticles(); @@ -123,10 +135,3 @@ void AliEMCALReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const delete et; } } - -AliTracker* AliEMCALReconstructor::CreateTracker(AliRunLoader* runLoader) const -{ -// creates the EMCAL tracker - if (!runLoader) return NULL; - return new AliEMCALTracker(runLoader); -} diff --git a/EMCAL/AliEMCALReconstructor.h b/EMCAL/AliEMCALReconstructor.h index dbe71ebf744..fad05183984 100644 --- a/EMCAL/AliEMCALReconstructor.h +++ b/EMCAL/AliEMCALReconstructor.h @@ -40,8 +40,7 @@ public: virtual ~AliEMCALReconstructor() ; //dtor Bool_t Debug() const { return fDebug ; } - AliTracker *CreateTracker(AliRunLoader* runLoader) const; - virtual void FillESD(AliRunLoader* runLoader, AliESD* esd) const ; + virtual void FillESD(AliRunLoader* runLoader, AliESD* esd) const ; virtual void Reconstruct(AliRunLoader* runLoader) const ; virtual void Reconstruct(AliRunLoader* runLoader, AliRawReaderFile * rawreader) const ; diff --git a/EMCAL/AliEMCALSDigitizer.cxx b/EMCAL/AliEMCALSDigitizer.cxx index 89c786dcd66..cf18702aeb5 100644 --- a/EMCAL/AliEMCALSDigitizer.cxx +++ b/EMCAL/AliEMCALSDigitizer.cxx @@ -50,7 +50,9 @@ // --- ROOT system --- #include "TBenchmark.h" - //#include "TObjectTable.h" +#include "TH1.h" +#include "TBrowser.h" +//#include "TObjectTable.h" // --- Standard library --- #include "stdlib.h" @@ -62,7 +64,8 @@ #include "AliEMCALHit.h" #include "AliEMCALSDigitizer.h" #include "AliEMCALGeometry.h" - //#include "AliMemoryWatcher.h" +#include "AliEMCALJetMicroDst.h" +//#include "AliMemoryWatcher.h" ClassImp(AliEMCALSDigitizer) @@ -70,8 +73,10 @@ ClassImp(AliEMCALSDigitizer) AliEMCALSDigitizer::AliEMCALSDigitizer():TTask("","") { // ctor - fFirstEvent = fLastEvent = 0 ; + Info("AliEMCALSDigitizer()", " CTOR # 1 #"); + fFirstEvent = fLastEvent = fControlHists = 0; fDefaultInit = kTRUE ; + fHists = 0; } //____________________________________________________________________________ @@ -81,16 +86,19 @@ AliEMCALSDigitizer::AliEMCALSDigitizer(const char * alirunFileName, fEventFolderName(eventFolderName) { // ctor - fFirstEvent = fLastEvent = 0 ; // runs one event by defaut + Info("AliEMCALSDigitizer()", " CTOR # 2 #"); + fFirstEvent = fLastEvent = fControlHists = 0 ; // runs one event by defaut Init(); InitParameters() ; fDefaultInit = kFALSE ; + fHists = 0; } //____________________________________________________________________________ AliEMCALSDigitizer::AliEMCALSDigitizer(const AliEMCALSDigitizer & sd) : TTask(sd) { //cpy ctor + Info("AliEMCALSDigitizer()", " CPY CTOR # 3 #"); fFirstEvent = sd.fFirstEvent ; fLastEvent = sd.fLastEvent ; @@ -144,26 +152,25 @@ void AliEMCALSDigitizer::Init(){ //____________________________________________________________________________ void AliEMCALSDigitizer::InitParameters() { - // Initializes parameters - fA = 0; - fB = 10000000.; - AliEMCALGetter * gime = AliEMCALGetter::Instance() ; const AliEMCALGeometry * geom = gime->EMCALGeometry() ; if (geom->GetSampling() == 0.) { Fatal("InitParameters", "Sampling factor not set !") ; } -// else -// Info("InitParameters", "Sampling factor set to %f", geom->GetSampling()) ; - - // this threshold corresponds approximately to 100 MeV - fECPrimThreshold = 100E-3; + // Initializes parameters + fA = 0; + fB = 1.e+7; + fSampling = geom->GetSampling(); + + // threshold for deposit energy of hit + fECPrimThreshold = 0.; // 24-nov-04 - was 1.e-6; } //____________________________________________________________________________ void AliEMCALSDigitizer::Exec(Option_t *option) { // Collects all hit of the same tower into digits + TString o(option); o.ToUpper(); if (strstr(option, "print") ) { Print() ; return ; @@ -183,16 +190,19 @@ void AliEMCALSDigitizer::Exec(Option_t *option) if (fLastEvent == -1) fLastEvent = gime->MaxEvent() - 1 ; - else - fLastEvent = TMath::Min(fFirstEvent, gime->MaxEvent()); // only ine event at the time + else { + if(o != "EXACT") fLastEvent = TMath::Min(fFirstEvent, gime->MaxEvent()); + fLastEvent = TMath::Min(fLastEvent, gime->MaxEvent()); + printf("AliEMCALSDigitizer::Exec : option: %s | %i -> %i events : Max events %i \n \n", + option, fFirstEvent, fLastEvent, gime->MaxEvent()); + } Int_t nEvents = fLastEvent - fFirstEvent + 1; - Int_t ievent ; - - //AliMemoryWatcher memwatcher; - - for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) { - gime->Event(ievent,"H") ; + Int_t ievent; + Float_t energy=0.; // de * fSampling - 23-nov-04 + for (ievent = fFirstEvent; ievent < fLastEvent; ievent++) { + if(ievent%100==0 || ievent==fLastEvent-1) printf(" processed event %i \n", ievent); + gime->Event(ievent,"XH"); // read primaries and hits onl TTree * treeS = gime->TreeS(); TClonesArray * hits = gime->Hits() ; TClonesArray * sdigits = gime->SDigits() ; @@ -207,28 +217,30 @@ void AliEMCALSDigitizer::Exec(Option_t *option) //=========== Get the EMCAL branch from Hits Tree for the Primary iprim gime->Track(iprim) ; Int_t i; + AliEMCALGeometry *geom = gime->EMCALGeometry(); for ( i = 0 ; i < hits->GetEntries() ; i++ ) { AliEMCALHit * hit = dynamic_cast(hits->At(i)) ; AliEMCALDigit * curSDigit = 0 ; AliEMCALDigit * sdigit = 0 ; Bool_t newsdigit = kTRUE; + // hit->GetId() - Absolute Id number EMCAL segment + if(geom->CheckAbsCellId(hit->GetId())) { // was IsInECA(hit->GetId()) + energy = hit->GetEnergy() * fSampling; // 23-nov-04 + if(energy > fECPrimThreshold ) // Assign primary number only if deposited energy is significant - AliEMCALGeometry * geom = gime->EMCALGeometry() ; - if( geom->IsInECA(hit->GetId()) ){ - if( hit->GetEnergy() > fECPrimThreshold ) curSDigit = new AliEMCALDigit( hit->GetPrimary(), hit->GetIparent(), hit->GetId(), - Digitize(hit->GetEnergy()), hit->GetTime() ) ; + Digitize(energy), hit->GetTime() ) ; else curSDigit = new AliEMCALDigit( -1 , -1 , hit->GetId(), - Digitize(hit->GetEnergy()), hit->GetTime() ) ; - } - else{ - newsdigit = kFALSE; - } + Digitize(energy), hit->GetTime() ) ; + } else { + Warning("Exec"," abs id %i is bad \n", hit->GetId()); + newsdigit = kFALSE; + } Int_t check = 0 ; @@ -256,12 +268,21 @@ void AliEMCALSDigitizer::Exec(Option_t *option) Int_t nPrimarymax = -1 ; Int_t i ; + Double_t e=0.,esum=0.; + sv::FillH1(fHists, 0, double(sdigits->GetEntriesFast())); for (i = 0 ; i < sdigits->GetEntriesFast() ; i++) { AliEMCALDigit * sdigit = dynamic_cast(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(sdigits->At(i)))->GetNprimary()) > nPrimarymax) nPrimarymax = ((dynamic_cast(sdigits->At(i)))->GetNprimary()) ; } @@ -277,15 +298,12 @@ void AliEMCALSDigitizer::Exec(Option_t *option) gime->WriteSDigits("OVERWRITE"); //NEXT - SDigitizer - - gime->WriteSDigitizer("OVERWRITE"); + if(ievent == fFirstEvent) gime->WriteSDigitizer("OVERWRITE"); // why in event cycle ? if(strstr(option,"deb")) PrintSDigits(option) ; - - //gObjectTable->Print() ; - //memwatcher.Watch(ievent); - }// event loop + } + // gime->WriteSDigitizer("OVERWRITE"); // 12-jan-04 Unload(); @@ -293,7 +311,7 @@ void AliEMCALSDigitizer::Exec(Option_t *option) if(strstr(option,"tim")){ gBenchmark->Stop("EMCALSDigitizer"); - printf("Exec: took %f seconds for SDigitizing %f seconds per event", + printf("\n Exec: took %f seconds for SDigitizing %f seconds per event\n", gBenchmark->GetCpuTime("EMCALSDigitizer"), gBenchmark->GetCpuTime("EMCALSDigitizer")/nEvents ) ; } @@ -303,17 +321,26 @@ void AliEMCALSDigitizer::Exec(Option_t *option) } +void AliEMCALSDigitizer::Print1(Option_t * option) +{ + Print(); + PrintSDigits(option); +} + //__________________________________________________________________ -void AliEMCALSDigitizer::Print()const +void AliEMCALSDigitizer::Print() const { // Prints parameters of SDigitizer - printf("Print: \n------------------- %s -------------", GetName() ) ; + printf("Print: \n------------------- %s -------------\n", GetName() ) ; + printf(" fInit %i\n", int(fInit)); + printf(" fFirstEvent %i\n", fFirstEvent); + printf(" fLastEvent %i\n", fLastEvent); printf(" Writing SDigits to branch with title %s\n", fEventFolderName.Data()) ; - printf(" with digitization parameters A = %f\n", fA) ; - printf(" B = %f\n", fB) ; - printf(" Threshold for EC Primary assignment= %f\n", fECPrimThreshold) ; + printf(" with digitization parameters A = %f\n", fA) ; + printf(" B = %f\n", fB) ; + printf(" Threshold for EC Primary assignment = %f\n", fECPrimThreshold) ; + printf(" Sampling = %f\n", fSampling); printf("---------------------------------------------------\n") ; - } //__________________________________________________________________ @@ -339,19 +366,20 @@ void AliEMCALSDigitizer::PrintSDigits(Option_t * option) printf("\n") ; printf("event %i", gAlice->GetEvNumber()) ; - printf("\n Number of entries in SDigits list %i", sdigits->GetEntriesFast()); + printf(" Number of entries in SDigits list %i", sdigits->GetEntriesFast()); if(strstr(option,"all")||strstr(option,"EMC")){ //loop over digits AliEMCALDigit * digit; printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ; - Int_t index ; + Int_t index, isum=0; char * tempo = new char[8192]; for (index = 0 ; index < sdigits->GetEntries() ; index++) { digit = dynamic_cast( 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; iprimaryGetNprimary(); iprimary++) { @@ -360,7 +388,8 @@ void AliEMCALSDigitizer::PrintSDigits(Option_t * option) } } delete tempo ; - } + printf("\n** Sum %i : %10.3f GeV/c **\n ", isum, double(isum)*1.e-6); + } else printf("\n"); } //____________________________________________________________________________ @@ -372,3 +401,31 @@ void AliEMCALSDigitizer::Unload() const loader->UnloadHits() ; loader->UnloadSDigits() ; } + +void AliEMCALSDigitizer::Browse(TBrowser* b) +{ + if(fHists) b->Add(fHists); + TTask::Browse(b); +} + +TList *AliEMCALSDigitizer::BookControlHists(const int var) +{ // 22-nov-04 + gROOT->cd(); + AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName); + const AliEMCALGeometry *geom = gime->EMCALGeometry() ; + if(var>=1){ + new TH1F("HSDigiN", "#EMCAL sdigits ", 1001, -0.5, 1000.5); + new TH1F("HSDigiSumEnergy","Sum.EMCAL energy", 1000, 0.0, 100.); + new TH1F("HSDigiAmp", "EMCAL sdigits amplitude", 1000, 0., 2.e+9); + new TH1F("HSDigiEnergy","EMCAL cell energy", 1000, 0.0, 100.); + new TH1F("HSDigiAbsId","EMCAL absID for sdigits", + geom->GetNCells(), 0.5, Double_t(geom->GetNCells())+0.5); + } + fHists = sv::MoveHistsToList("EmcalSDigiControlHists", kFALSE); + return fHists; +} + +void AliEMCALSDigitizer::SaveHists(const char* name, Bool_t kSingleKey, const char* opt) +{ + sv::SaveListOfHists(fHists, name, kSingleKey, opt); +} diff --git a/EMCAL/AliEMCALSDigitizer.h b/EMCAL/AliEMCALSDigitizer.h index 4429ce5c5b1..ec5065bb12b 100644 --- a/EMCAL/AliEMCALSDigitizer.h +++ b/EMCAL/AliEMCALSDigitizer.h @@ -17,6 +17,9 @@ // --- ROOT system --- #include "TTask.h" class TFile ; +class TList; +class TBrowser; +//class TBrowser; // --- Standard library --- @@ -32,16 +35,25 @@ public: virtual ~AliEMCALSDigitizer(); // dtor Float_t Calibrate(Int_t amp)const {return (amp - fA)/fB ; } - Int_t Digitize(Float_t Energy)const { return (Int_t ) ( fA + Energy*fB); } + Int_t Digitize(Float_t energy)const { return (Int_t ) (fA + energy*fB); } virtual void Exec(Option_t *option); Int_t GetSDigitsInRun() const {return fSDigitsInRun ;} - virtual void Print() const ; + virtual void Print() const; + void Print1(Option_t *option="all"); // *MENU* void SetEventFolderName(TString name) { fEventFolderName = name ; } void SetEventRange(Int_t first=0, Int_t last=-1) {fFirstEvent=first; fLastEvent=last; } Bool_t operator == (const AliEMCALSDigitizer & sd) const ; const AliEMCALSDigitizer & operator = (const AliEMCALSDigitizer & /*sd*/) {return *this ;} + virtual void Browse(TBrowser* b); + // hists + void SetControlHists(const Int_t var=0) {fControlHists=var;} + Int_t GetControlHist() const {return fControlHists;} + TList *GetListOfHists() {return fHists;} + TList* BookControlHists(const int var=0); + void SaveHists(const char* name="RF/TRD1/Digitizations/SDigiVar?", + Bool_t kSingleKey=kTRUE, const char* opt="RECREATE"); // *MENU* private: void Init() ; @@ -55,13 +67,16 @@ private: Float_t fECPrimThreshold ; // To store primary if EC Shower Elos > threshold Bool_t fDefaultInit; //! Says if the task was created by defaut ctor (only parameters are initialized) TString fEventFolderName; // event folder name - Bool_t fInit ; //! tells if initialisation wennt OK, will revent exec if not + Bool_t fInit ; //! tells if initialisation went OK, will revent exec if not Int_t fSDigitsInRun ; //! Total number of sdigits in one run Int_t fFirstEvent; // first event to process Int_t fLastEvent; // last event to process + Float_t fSampling; // See AliEMCALGeometry + // Control hists + Int_t fControlHists; //! + TList *fHists; //! ClassDef(AliEMCALSDigitizer,5) // description - }; #endif // AliEMCALSDigitizer_H diff --git a/EMCAL/AliEMCALShishKebabModule.cxx b/EMCAL/AliEMCALShishKebabModule.cxx new file mode 100644 index 00000000000..94e89cdc8df --- /dev/null +++ b/EMCAL/AliEMCALShishKebabModule.cxx @@ -0,0 +1,183 @@ +/************************************************************************** + * Copyright(c) 1998-2004, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ + +//*-- Author: Aleksei Pavlinov(WSU) +#include "AliEMCALShishKebabModule.h" +#include "AliEMCALGeometry.h" +#include +#include + +#include + +ClassImp(AliEMCALShishKebabModule) + + AliEMCALGeometry *AliEMCALShishKebabModule::fgGeometry=0; + Double_t AliEMCALShishKebabModule::fga=0.; + Double_t AliEMCALShishKebabModule::fgb=0.; + Double_t AliEMCALShishKebabModule::fgr=0.; + +AliEMCALShishKebabModule::AliEMCALShishKebabModule(const double theta) : TNamed() +{ // theta in radians ; first object shold be with theta=pi/2. + fTheta = theta; + if(fgGeometry==0) { + if(GetParameters()) { + DefineFirstModule(); + } + } else Warning("AliEMCALShishKebabModule(theta)","You should call this constractor just once !!"); + DefineName(fTheta); +} + +AliEMCALShishKebabModule::AliEMCALShishKebabModule(AliEMCALShishKebabModule &leftNeighbor) : TNamed() +{ // 22-sep-04 + TObject::SetUniqueID(leftNeighbor.GetUniqueID()+1); + Init(leftNeighbor.GetA(),leftNeighbor.GetB()); +} + +void AliEMCALShishKebabModule::Init(const double A,const double B) +{ + Double_t thetaMin, thetaMax, par[4]; + Int_t npar=0; + if(A<0){ + // DefineSecondModuleFirstAssumption(); + thetaMax = TMath::ATan2(fgr, 1.4*fga); + thetaMin = TMath::ATan2(fgr, 1.6*fga); + fTheta = Solve(AliEMCALShishKebabModule::Y2, thetaMin,thetaMax,npar,par); + } else{ + npar = 4; + par[0] = fga; + par[1] = fgr; + par[2] = A; + par[3] = B; + Double_t x = fgr/A; + thetaMax = TMath::ATan2(fgr,x + 0.5*fga); + thetaMin = TMath::ATan2(fgr,x + 1.5*fga); + fTheta = Solve(AliEMCALShishKebabModule::YALL, thetaMin,thetaMax,npar,par); + } + + Double_t rOK = fgr / TMath::Sin(fTheta) + (fga/2.)/TMath::Tan(fTheta) + fgb/2.; + fOK.SetMagPhi(rOK, fTheta); + // have to define A and B + fA = TMath::Tan(fTheta); + fB = -fga/(2.*TMath::Cos(fTheta)); + DefineName(fTheta); +} + +void AliEMCALShishKebabModule::DefineFirstModule() +{ + fOK.Set(fga/2., fgr + fgb/2.); // position the center of module vs o + + fB = fga/2.; // z=fB + fA = -999.; // fA=infinite in this case + TObject::SetUniqueID(1); // +} + +void AliEMCALShishKebabModule::DefineSecondModuleFirstAssumption() +{ // Keep for testing and checking + // cos(theta) << 1, theta ~ pi/2.; a/r = 11.4/462.54 = 0.0246465 << 1; + // theta=1.53382 from this case; theta=1.533869 from TGraph::Zero + Double_t x = (3.*fga)/(2.*fgr); + fTheta = TMath::ACos(x); + /* + Double_t rOK = fgr / TMath::Sin(fTheta) + (fga/2.)/TMath::Tan(fTheta) + fgb/2.; + fOK.SetMagPhi(rOK, fTheta); + // have to define A and B + fA = TMath::Tan(fTheta); + fB = -fga/(2.*TMath::Cos(fTheta)); + DefineName(fTheta); + */ +} + +Double_t AliEMCALShishKebabModule::Solve(Double_t (*fcn)(Double_t*,Double_t*), +Double_t xmin, Double_t xmax, Int_t npar, Double_t *par, Double_t eps, Int_t maxIter) +{ + if(npar); // unused now + TGraph gr; + double X,Y; + Int_t k = 0; + gr.Zero(k, xmin,xmax, eps, X,Y, maxIter); // remember initial interval + while(k!=2) { + Y = fcn(&X, par); + gr.Zero(k, xmin,xmax, eps, X,Y, maxIter); + } + return X; +} + +Double_t AliEMCALShishKebabModule::Y2(double *x, double *par) +{ // For position calulation of second module + if(par); + Double_t theta = x[0]; + Double_t cos = TMath::Cos(theta); + Double_t sin = TMath::Sin(theta); + Double_t y1 = fgr*cos/sin + fga/(2.*sin) - fga*sin; + Double_t y2 = fga, y = y1-y2;; + // printf(" theta %f Y %12.5e \n", theta, y); + return y; +} + +Double_t AliEMCALShishKebabModule::YALL(double *x, double *par) +{ // For position calulation of 3th, 4th to 30th modules + Double_t a=par[0], r=par[1], A=par[2], B=par[3]; + Double_t theta = x[0]; + Double_t cos = TMath::Cos(theta); + Double_t sin = TMath::Sin(theta); + + Double_t y1 = r + a*cos; + Double_t y2 = A*(r*cos/sin + a/(2.*sin) - a*sin) + B; + Double_t y = y1-y2; + // printf(" theta %f Y %12.5e \n", theta, y); + return y; +} + +void AliEMCALShishKebabModule::DefineName(const double theta) +{ + char name[100]; + // sprintf(name,"theta_%5.2f",theta*180./TMath::Pi()); + sprintf(name,"%2i(%5.2f)", TObject::GetUniqueID(), theta*180./TMath::Pi()); + SetName(name); +} + +Bool_t AliEMCALShishKebabModule::GetParameters() +{ + fgGeometry = AliEMCALGeometry::GetInstance(); + // if(!fgGeometry) assert(0); + if(!fgGeometry) { + Warning("GetParameters()"," No geometry "); + return kFALSE; + } + + fga = (Double_t)fgGeometry->GetPhiModuleSize(); + fgb = (Double_t)fgGeometry->GetLongModuleSize(); + fgr = (Double_t)(fgGeometry->GetIPDistance() + fgGeometry->GetSteelFrontThickness()); + Print(0); + return kTRUE; +} + +// service methods +void AliEMCALShishKebabModule::Print(const int pri) const +{ + if(pri>=0) { + Info("Print()", " a %7.2f | b %7.2f | r %7.2f ", fga, fgb, fgr); + printf(" fTheta %f : %5.2f : cos(theta) %f\n", fTheta, GetThetaInDegree(),TMath::Cos(fTheta)); + if(pri>0) { + printf("%i %s | theta %f -> %f\n", GetUniqueID(), GetName(), fTheta, fOK.Phi()); + printf(" A %f B %f \n", fA, fB); + + fOK.Dump(); + } + } +} + diff --git a/EMCAL/AliEMCALShishKebabModule.h b/EMCAL/AliEMCALShishKebabModule.h new file mode 100644 index 00000000000..59a9e3629c9 --- /dev/null +++ b/EMCAL/AliEMCALShishKebabModule.h @@ -0,0 +1,63 @@ +#ifndef ALIEMCALSHISHKEBABMODULE_H +#define ALIEMCALSHISHKEBABMODULE_H + +/* Copyright(c) 1998-2004, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id$ */ +#include "TNamed.h" +#include "TMath.h" +#include "TVector2.h" + +class AliEMCALGeometry; + +class AliEMCALShishKebabModule : public TNamed { + public: + AliEMCALShishKebabModule(const double theta=TMath::Pi()/2.); + AliEMCALShishKebabModule(AliEMCALShishKebabModule &leftNeighbor); + void Init(const double A,const double B); + + virtual ~AliEMCALShishKebabModule(void) {} + Bool_t GetParameters(); + void DefineName(const double theta); + void DefineFirstModule(); + void DefineSecondModuleFirstAssumption(); // need for testing + + Double_t Solve(Double_t (*fcn)(Double_t*, Double_t*), Double_t xmin=0., Double_t xmax=1., + Int_t npar=0, Double_t *par=0, Double_t eps=1.0e-8, Int_t maxIter=1000); + + static Double_t Y2(double *x, double *par); + static Double_t YALL(double *x, double *par); + + Double_t GetTheta() const {return fTheta;} + Double_t GetThetaInDegree() const {return fTheta*180./TMath::Pi();} + + Double_t GetPosX() {return fOK.Y();} + Double_t GetPosZ() {return fOK.X();} + Double_t GetPosXfromR() {return fOK.Y() - fgr;} + Double_t GetA() {return fA;} + Double_t GetB() {return fB;} + + // geometry info + static AliEMCALGeometry *fgGeometry; //! + static Double_t fga; // default 11.2cm + static Double_t fgb; // + // radius to IP + static Double_t fgr; + + TVector2 fOK; // position the module center x->y; z->x; + Double_t fA; // parameters of line = y = A*z + B + Double_t fB; // + // service methods + void Print(const int pri=1) const; // *MENU* + protected: + // size of SK module + Double_t fTheta; // theta for SK module + + ClassDef(AliEMCALShishKebabModule,0) // Turned Shish-Kebab module +}; + +#endif +/* To do + 1. Insert position the center of towers - 2 additional TVector2 + */ diff --git a/EMCAL/AliEMCALShishKebabTrd1Module.cxx b/EMCAL/AliEMCALShishKebabTrd1Module.cxx new file mode 100644 index 00000000000..f0d07d39cd5 --- /dev/null +++ b/EMCAL/AliEMCALShishKebabTrd1Module.cxx @@ -0,0 +1,155 @@ +/************************************************************************** + * Copyright(c) 1998-2004, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ + +//*-- Author: Aleksei Pavlinov(WSU) + +#include "AliEMCALShishKebabTrd1Module.h" +#include +#include +#include "AliEMCALGeometry.h" + +ClassImp(AliEMCALShishKebabTrd1Module) + + AliEMCALGeometry *AliEMCALShishKebabTrd1Module::fgGeometry=0; + Double_t AliEMCALShishKebabTrd1Module::fga=0.; + Double_t AliEMCALShishKebabTrd1Module::fga2=0.; + Double_t AliEMCALShishKebabTrd1Module::fgb=0.; + Double_t AliEMCALShishKebabTrd1Module::fgr=0.; + Double_t AliEMCALShishKebabTrd1Module::fgangle=0.; // one degrre + Double_t AliEMCALShishKebabTrd1Module::fgtanBetta=0; // + +AliEMCALShishKebabTrd1Module::AliEMCALShishKebabTrd1Module(const double theta) : TNamed() +{ // theta in radians ; first object shold be with theta=pi/2. + fTheta = theta; + if(fgGeometry==0) { + if(GetParameters()) { + DefineFirstModule(); + } + } else Warning("AliEMCALShishKebabTrd1Module(theta)","You should call this constractor just once !!"); + DefineName(fTheta); +} + +AliEMCALShishKebabTrd1Module::AliEMCALShishKebabTrd1Module(AliEMCALShishKebabTrd1Module &leftNeighbor) : TNamed() +{ // 22-sep-04 + // printf("** Left Neighbor : %s **\n", leftNeighbor.GetName()); + TObject::SetUniqueID(leftNeighbor.GetUniqueID()+1); + fTheta = leftNeighbor.GetTheta() - fgangle; + Init(leftNeighbor.GetA(),leftNeighbor.GetB()); +} + +void AliEMCALShishKebabTrd1Module::Init(const double A,const double B) +{ // Define parameter module from parameters A,B from previos. + Double_t yl = (fgb/2)*TMath::Sin(fTheta) + (fga/2)*TMath::Cos(fTheta) + fgr, y = yl; + Double_t xl = (yl - B) / A; // y=A*x+B + + // Double_t xp1 = (fga/2. + fgb/2.*fgtanBetta)/(TMath::Sin(fTheta) + fgtanBetta*TMath::Cos(fTheta)); + // printf(" xp1 %9.3f \n ", xp1); + // xp1 == xp => both methods give the same results - 3-feb-05 + Double_t alpha = TMath::Pi()/2. + fgangle/2; + Double_t xt = (fga+fga2)*TMath::Tan(fTheta)*TMath::Tan(alpha)/(4.*(1.-TMath::Tan(fTheta)*TMath::Tan(alpha))); + Double_t yt = xt / TMath::Tan(fTheta), xp = TMath::Sqrt(xt*xt + yt*yt); + Double_t x = xl + xp; + fOK.Set(x, y); + // printf(" yl %9.3f | xl %9.3f | xp %9.3f \n", yl, xl, xp); + + // have to define A and B; + Double_t yCprev = fgr + fga*TMath::Cos(fTheta); + Double_t xCprev = (yCprev - B) / A; + Double_t xA = xCprev + fga*TMath::Sin(fTheta), yA = fgr; + + fThetaA = fTheta - fgangle/2.; + fA = TMath::Tan(fThetaA); // !! + fB = yA - fA*xA; + + DefineName(fTheta); + // Centers of modules + Double_t kk1 = (fga+fga2)/(2.*4.); // kk1=kk2 + + Double_t xk1 = fOK.X() - kk1*TMath::Sin(fTheta); + Double_t yk1 = fOK.Y() + kk1*TMath::Cos(fTheta); + fOK1.Set(xk1,yk1); + + Double_t xk2 = fOK.X() + kk1*TMath::Sin(fTheta); + Double_t yk2 = fOK.Y() - kk1*TMath::Cos(fTheta); + fOK2.Set(xk2,yk2); +} + +void AliEMCALShishKebabTrd1Module::DefineFirstModule() +{ + fOK.Set(fga2/2., fgr + fgb/2.); // position the center of module vs o + + // parameters of right line : y = A*z + B in system where zero point is IP. + fThetaA = fTheta - fgangle/2.; + fA = TMath::Tan(fThetaA); + Double_t xA = fga/2. + fga2/2., yA = fgr; + fB = yA - fA*xA; + + Double_t kk1 = (fga+fga2)/(2.*4.); // kk1=kk2 + fOK1.Set(fOK.X() - kk1, fOK.Y()); + fOK2.Set(fOK.X() + kk1, fOK.Y()); + + TObject::SetUniqueID(1); // +} + +void AliEMCALShishKebabTrd1Module::DefineName(const double theta) +{ + char name[100]; + // sprintf(name,"theta_%5.2f",theta*180./TMath::Pi()); + sprintf(name,"%2i(%5.2f)", TObject::GetUniqueID(), theta*180./TMath::Pi()); + SetName(name); +} + +Bool_t AliEMCALShishKebabTrd1Module::GetParameters() +{ + fgGeometry = AliEMCALGeometry::GetInstance(); + TString sn(fgGeometry->GetName()); // 2-Feb-05 + sn.ToUpper(); + if(!fgGeometry) { + Warning("GetParameters()"," No geometry "); + return kFALSE; + } + + fga = (Double_t)fgGeometry->GetEtaModuleSize(); + fgb = (Double_t)fgGeometry->GetLongModuleSize(); + fgangle = Double_t(fgGeometry->GetTrd1Angle())*TMath::DegToRad(); + fgtanBetta = TMath::Tan(fgangle/2.); + fgr = (Double_t)fgGeometry->GetIPDistance(); + if(!sn.Contains("TRD2")) fgr += fgGeometry->GetSteelFrontThickness(); + fga2 = Double_t(fgGeometry->Get2Trd1Dx2()); + Print(0); + return kTRUE; +} + +// service methods +void AliEMCALShishKebabTrd1Module::Print(const int pri) const +{ + if(pri>=0) { + Info("Print()", "\n a %7.3f:%7.3f | b %7.2f | r %7.2f \n TRD1 angle %7.6f(%5.2f) | tanBetta %7.6f", + fga, fga2, fgb, fgr, fgangle, fgangle*TMath::RadToDeg(), fgtanBetta); + printf(" fTheta %f : %5.2f : cos(theta) %f\n", + fTheta, GetThetaInDegree(),TMath::Cos(fTheta)); + if(pri>0) { + printf("\n%i |%s| theta %f : fOK.Phi = %f(%5.2f)\n", + GetUniqueID(), GetName(), fTheta, fOK.Phi(),fOK.Phi()*TMath::RadToDeg()); + printf(" A %f B %f | fThetaA %7.6f(%5.2f)\n", fA,fB, fThetaA,fThetaA*TMath::RadToDeg()); + fOK.Dump(); + printf(" fOK : X %8.3f: Y %8.3f \n", fOK.X(), fOK.Y()); + printf(" fOK1 : X %8.3f: Y %8.3f \n", fOK1.X(), fOK1.Y()); + printf(" fOK2 : X %8.3f: Y %8.3f \n", fOK2.X(), fOK2.Y()); + } + } +} diff --git a/EMCAL/AliEMCALShishKebabTrd1Module.h b/EMCAL/AliEMCALShishKebabTrd1Module.h new file mode 100644 index 00000000000..e23240aeb6a --- /dev/null +++ b/EMCAL/AliEMCALShishKebabTrd1Module.h @@ -0,0 +1,75 @@ +#ifndef ALIEMCALSHISHKEBABTRD1MODULE_H +#define ALIEMCALSHISHKEBABTRD1MODULE_H + +/* Copyright(c) 1998-2004, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id$ */ +//*-- Author: Aleksei Pavlinov (WSU) +// TO DO : create class for Super module Geometry - 4-nov-04 + +#include "TNamed.h" +#include "TMath.h" +#include "TVector2.h" + +class AliEMCALGeometry; + +class AliEMCALShishKebabTrd1Module : public TNamed { + public: + AliEMCALShishKebabTrd1Module(const double theta=TMath::Pi()/2.); + AliEMCALShishKebabTrd1Module(AliEMCALShishKebabTrd1Module &leftNeighbor); + void Init(const double A,const double B); + + virtual ~AliEMCALShishKebabTrd1Module(void) {} + Bool_t GetParameters(); + void DefineName(const double theta); + void DefineFirstModule(); + + Double_t GetTheta() const{return fTheta;} + Double_t GetThetaInDegree() const {return fTheta*180./TMath::Pi();} + TVector2& GetCenterOfModule() {return fOK;} + Double_t GetEtaOfCenterOfModule(){return -TMath::Log(TMath::Tan(fOK.Phi()/2.));} + + Double_t GetPosX() {return fOK.Y();} + Double_t GetPosZ() {return fOK.X();} + Double_t GetPosXfromR() {return fOK.Y() - fgr;} + Double_t GetA() {return fA;} + Double_t GetB() {return fB;} + // Additional offline staff + TVector2& GetCenterOfCell(const Int_t ieta) + { if(ieta<=1) return fOK1; + else return fOK2;} + // + Double_t GetTanBetta() {return fgtanBetta;} + Double_t Getb() {return fgb;} + // service methods + void Print(const int pri=1) const; // *MENU* + + // geometry info + static AliEMCALGeometry *fgGeometry; //! + static Double_t fga; // 2*dx1=2*dy1 + static Double_t fga2; // 2*dx2 + static Double_t fgb; // 2*dz1 + static Double_t fgangle; // ~1 degree + static Double_t fgtanBetta; // tan(fgangle/2.) + // radius to IP + static Double_t fgr; + + protected: + TVector2 fOK; // position the module center x->y; z->x; + Double_t fA; // parameters of right line : y = A*z + B + Double_t fB; // system where zero point is IP. + Double_t fThetaA; // angle coresponding fA - for convinience + Double_t fTheta; // theta angle of perependicular to SK module + // position of towers with differents ieta (1 or 2) - 4-nov-04 + TVector2 fOK1; + TVector2 fOK2; + + public: + ClassDef(AliEMCALShishKebabTrd1Module,0) // Turned Shish-Kebab module +}; + +#endif +/* To do + 1. Insert position the center of towers - 2 additional TVector2 + */ diff --git a/EMCAL/AliEMCALWsuCosmicRaySetUp.cxx b/EMCAL/AliEMCALWsuCosmicRaySetUp.cxx new file mode 100644 index 00000000000..47df67e3435 --- /dev/null +++ b/EMCAL/AliEMCALWsuCosmicRaySetUp.cxx @@ -0,0 +1,132 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ + +/////////////////////////////////////////////////////////////////////////////// +// // +// Wsu Cosmic Ray SetUp // +// This class contains the description of the Wsu Cosmic Ray SetUp // +// external volume // +// // +//Begin_Html +/* + + +
+ +

The responsible person for this module is +Aleksei Pavlino, WSU. + +

+*/
+//End_Html
+//                                                                           //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include 
+
+#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
+  /*
+    
+  */
+  //End_Html
+
+  Float_t dASUC[3];
+  Int_t *idtmed = fIdtmed->GetArray()+1;
+  int idSC = idtmed[0];
+  //
+  dASUC[0]=50;
+  dASUC[1]=50;
+  dASUC[2]=50;
+  //  TString tmp(GetTitle());
+  gMC->Gsvolu(GetName(),"BOX",idSC, dASUC,3); // WSUC - Wsu Cosmic Ray SetUp
+}
+ 
+//_____________________________________________________________________________
+void AliEMCALWsuCosmicRaySetUp::CreateMaterials()
+{
+// Create materials and media
+  Int_t   isxfld = 0;
+  Float_t sxmgmx = 0.;
+  
+  // AIR
+  Float_t aAir[4]={12.0107,14.0067,15.9994,39.948};
+  Float_t zAir[4]={6.,7.,8.,18.};
+  Float_t wAir[4]={0.000124,0.755267,0.231781,0.012827};
+  Float_t dAir = 1.20479E-3;
+  //  Float_t dAir1 = 1.20479E-10;
+  //
+  AliMixture(1,"Air     $",aAir,zAir,dAir,4,wAir);
+  //
+  AliMedium(1,"Air     $",1,0,isxfld,sxmgmx,10,-1,-0.1,0.1 ,-10);
+}
+ 
+//_____________________________________________________________________________
+void AliEMCALWsuCosmicRaySetUp::DrawWSUC(float cxy) const
+{
+  //
+  // Draw a view of the Wsu Cosmic Ray SetUp 
+  //
+  // Set everything unseen
+  gMC->Gsatt("*", "seen", -1);
+  // 
+  // Set WSUC mother visible
+  gMC->Gsatt("WSUC","SEEN",1);
+  //
+  // Set the volumes visible
+  //
+  gMC->Gdopt("hide","off");
+
+  gMC->Gdraw("WSUC", 40, 30, 0, 10, 9, cxy, cxy);
+  gMC->Gdhead(1111, "WSU Cosmic Ray Setup ");
+
+  gMC->Gdman(18, 4, "MAN");
+}
+ 
+
diff --git a/EMCAL/AliEMCALWsuCosmicRaySetUp.h b/EMCAL/AliEMCALWsuCosmicRaySetUp.h
new file mode 100644
index 00000000000..b3fb8ce2b7b
--- /dev/null
+++ b/EMCAL/AliEMCALWsuCosmicRaySetUp.h
@@ -0,0 +1,29 @@
+#ifndef ALIBODY_H
+#define ALIBODY_H
+/* Copyright(c) 1998-2005, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////////////////
+//  Manager class for detector:  AliEMCALWsuCosmicRaySetUp        //
+//   This is the envelop for Alice                                //
+///////////////////////////////////////////////////////////////////
+ 
+#include "AliModule.h"
+
+class AliEMCALWsuCosmicRaySetUp : public AliModule {
+ 
+public:
+  AliEMCALWsuCosmicRaySetUp();
+  AliEMCALWsuCosmicRaySetUp(const char *name, const char *title);
+  virtual     ~AliEMCALWsuCosmicRaySetUp() {}
+  virtual void  CreateGeometry();
+  virtual void  CreateMaterials();
+  virtual Int_t IsVersion() const {return 0;}
+  void  DrawWSUC(float cxy=0.025) const; // *MENU*
+
+  ClassDef(AliEMCALWsuCosmicRaySetUp,1)  // Class manager for the Wsu Cosmic Ray SetUp
+};
+
+#endif
diff --git a/EMCAL/AliEMCALv0.cxx b/EMCAL/AliEMCALv0.cxx
index 5e17ae6cb88..9d4e2a9dd4e 100644
--- a/EMCAL/AliEMCALv0.cxx
+++ b/EMCAL/AliEMCALv0.cxx
@@ -30,12 +30,24 @@
 
 // --- 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  // can not include - I don't know why
 
 // --- Standard library ---
 
@@ -50,13 +62,20 @@
 
 ClassImp(AliEMCALv0)
 
+TArrayF ENVELOP1; // now global for BuildGeometry() - 30-aug-2004
+int idAIR=1599, idPB = 1600, idSC = 1601, idSTEEL = 1603; // global
+Int_t *idtmed=0, idrotm=0;
+double sampleWidth=0.;
+double parEMOD[5], smodPar0=0., smodPar1=0., smodPar2=0.;
+
+
 //______________________________________________________________________
 AliEMCALv0::AliEMCALv0(const char *name, const char *title):
   AliEMCAL(name,title)
 {
   // ctor : title is used to identify the layout
   GetGeometry() ; 
-
+  fShishKebabModules = 0;
 }
 
 //______________________________________________________________________
@@ -67,23 +86,153 @@ void AliEMCALv0::BuildGeometry()
     const Int_t kColorArm1   = kBlue ;
 
     AliEMCALGeometry * geom = GetGeometry() ; 
+    TString gn(geom->GetName());
+    gn.ToUpper(); 
 
     // Define the shape of the Calorimeter 
-    TNode * top = gAlice->GetGeometry()->GetNode("alice") ;
-    new TTUBS("Envelop1", "Tubs that contains arm 1", "void", 
-	      geom->GetEnvelop(0),     // rmin 
-	      geom->GetEnvelop(1) +30 ,// rmax
+    TNode * top = gAlice->GetGeometry()->GetNode("alice") ; // See AliceGeom/Nodes
+    TNode * envelopNode = 0;
+    char *envn = "Envelop1";
+    if(!gn.Contains("SHISH") || gn.Contains("TRD2")){
+      new TTUBS(envn, "Tubs that contains arm 1", "void", 
+	      geom->GetEnvelop(0) -10, // rmin 
+	      geom->GetEnvelop(1) +40 ,// rmax
 	      geom->GetEnvelop(2)/2.0, // half length in Z
 	      geom->GetArm1PhiMin(),   // minimum phi angle
 	      geom->GetArm1PhiMax()    // maximum phi angle
 	);
+      top->cd();
+      envelopNode = new TNode(envn, "Arm1 Envelop", "Envelop1", 0., 0., 0., "") ;
+    } else {
+      if(gn.Contains("WSUC")) {
+        envelopNode = BuildGeometryOfWSUC();
+      } else { // Shish-kebab now for compact, twist and TRD1 cases (ALIC)
+        envn="Envelop2";
+        TPGON *pgon = new TPGON(envn, "PGON that contains arm 1", "void", 
+        geom->GetArm1PhiMin(),geom->GetArm1PhiMax()-geom->GetArm1PhiMin(),geom->GetNPhiSuperModule(), 2);
+      // define section
+        pgon->DefineSection(0, ENVELOP1[4],  ENVELOP1[5], ENVELOP1[6]);
+        pgon->DefineSection(1, ENVELOP1[7],  ENVELOP1[5], ENVELOP1[6]);
+        top->cd();
+        envelopNode = new TNode(envn, "Arm1 Envelop2", envn, 0., 0., 0., "") ;
+      }
+    }
+				
+    envelopNode->SetLineColor(kColorArm1) ;
+    fNodes->Add(envelopNode);
+}
+
+TNode *AliEMCALv0::BuildGeometryOfWSUC()
+{ // June 8, 2005; see directory geant3/TGeant3/G3toRoot.cxx
+  // enum EColor { kWhite, kBlack, kRed, kGreen, kBlue, kYellow, kMagenta, kCyan } - see $ROOTSYS/include/Gtypes.h
+   AliEMCALGeometry * g = GetGeometry(); 
+   TNode * top = gAlice->GetGeometry()->GetNode("alice") ; // See AliceGeom/Nodes
+   top->cd();
+
+   TNode *envelopNode = 0;
+   char *name = "";
+   /*
+    name = "WSUC";
+   new TBRIK(name, "WSUC(XEN1 in Geant)","void",ENVELOP1[0],ENVELOP1[1],ENVELOP1[2]);
+    envelopNode = new TNode(name, "envelope for WSUC", name, 0., 0., 0., "");
+   envelopNode->SetVisibility(0);
+   */
+
+   TNode *emod=0, *scmx=0;
+   name = "SMOD"; // super module
+   new TBRIK(name, "SMOD(SMOD in Geant)","void", smodPar0,smodPar1,smodPar2);
+   if(envelopNode) envelopNode->cd();
+   TNode *smod = new TNode(name, "SMOD", name, 0., 0., 0., "");   
+   smod->SetLineColor(kBlue) ;
+   if(envelopNode==0) envelopNode = smod;
+
+   name = "EMOD"; // see CreateEMOD 
+   TTRD1 *EMOD = new TTRD1(name, "EMOD(EMOD in Geant)","void", float(parEMOD[0]),
+   float(parEMOD[1]),float(parEMOD[2]),float(parEMOD[3]));
+
+   // SCMX = EMOD/4 for simplicity of drawing
+   name = "SCMX";
+   Float_t dz=0.,theta=0.,phi=0.,h1=0.,bl1=0.,tl1=0.,alpha1=0.,h2=0.,bl2=0.,tl2=0.,alpha2=0.;
+   h1     = EMOD->GetDy()/2.;
+   bl1    = EMOD->GetDx()/2.;
+   tl1    = bl1;
+   alpha1 = 0.;
+   h2     = EMOD->GetDy()/2.;
+   bl2    = EMOD->GetDx2()/2.;
+   tl2    = bl2;
+   alpha2 = 0.;
+
+   dz       = EMOD->GetDz();
+   double dr = TMath::Sqrt((h2-h1)*(h2-h1)+(bl2-bl1)*(bl2-bl1));
+   theta    = TMath::ATan2(dr,2.*dz) * TMath::RadToDeg(); 
+   phi      = 180.;
+
+   TTRAP *SCMX = new TTRAP(name, "SCMX(SCMX as in Geant)","void",
+   	     dz,theta,phi, h1,bl1,tl1,alpha1, h2,bl2,tl2,alpha2);
+   //   SCMX->Dump();
+   Float_t xShiftSCMX = (EMOD->GetDx() + EMOD->GetDx2())/4.;
+   Float_t yShiftSCMX = EMOD->GetDy()/2.;
+   printf(" xShiftSCMX %7.4f yShiftSCMX %7.4f \n",xShiftSCMX,yShiftSCMX);
+
+   name = "EMOD"; // see CreateEMOD 
+   smod->cd();
+   
+   AliEMCALShishKebabTrd1Module *mod=0;
+   Double_t  angle=90., xpos=0.,ypos=0.,zpos=0., xposSCMX=0.,yposSCMX=0.,zposSCMX=0.;
+   char rtmn[100], rtmt[100];
+   TRotMatrix *rtm=0, *rtmSCMX=0;
+   int numEmod=0;
+   for(int iz=0; izGetNZ(); 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; ixGetNPhi(); ix++) { // flat in phi
+       xpos = g->GetPhiModuleSize()*(2*ix+1 - g->GetNPhi())/2.;
+       sprintf(rtmt,"EMOD, iz %i, ix %i, angle %6.3f",iz,ix, angle);
+       TString namNode=name;
+       namNode += numEmod++;
+       smod->cd();
+       emod = new TNode(namNode.Data(), rtmt, (TShape*)EMOD, xpos,ypos,zpos,rtm);   
+       //       emod->SetLineColor(kGreen) ;
+       emod->SetVisibility(0); // SCMX will bi visible 
+       if(SCMX) { // 4(2x2) sensetive volume inside EMOD
+         emod->cd();
+         zposSCMX = 0.;
+         for(int jy=0; jy<2; jy++){ // division on y
+           yposSCMX = yShiftSCMX *(2*jy - 1); 
+           for(int jx=0; jx<2; jx++){ // division on x
+             Double_t theta1=90.,phi1=0., theta2=90.,phi2=90., theta3=0.,phi3=0 ; 
+             xposSCMX = xShiftSCMX *(2*jx - 1);              
+             namNode = "SCMX";
+             namNode += jy;
+             namNode += jx;
+             sprintf(rtmn,"rm%s",namNode.Data());
+             sprintf(rtmt,"rotation matrix for %s inside EMOD",namNode.Data());
+             rtmSCMX = tg->GetRotMatrix(rtmn);
+             if(jx == 1) {
+               phi1 = 180.; // x' = -x
+               phi2 = 270.; // y' = -y
+             }
+             if(rtmSCMX == 0) rtmSCMX = new TRotMatrix(rtmn,rtmt, theta1,phi1, theta2,phi2, theta3,phi3);
+             sprintf(rtmt,"%s inside %s", namNode.Data(), emod->GetName());
+             scmx = new TNode(namNode.Data(), rtmt, (TShape*)SCMX, xposSCMX,yposSCMX,zposSCMX,rtmSCMX);   
+             scmx->SetLineColor(kMagenta);
+           }
+         }
+       }
+     }
+   }
+   //   emod->Draw(); // for testing
+
+   return envelopNode;
 }
 
 //______________________________________________________________________
@@ -118,44 +267,82 @@ void AliEMCALv0::CreateGeometry()
     Float_t *dum=0;
 
     AliEMCALGeometry * geom = GetGeometry() ; 
+    TString gn(geom->GetName());
+    gn.ToUpper(); 
 
     if(!(geom->IsInitialized())){
 	Error("CreateGeometry","EMCAL Geometry class has not been set up.");
     } // end if
 
     // Get pointer to the array containing media indices
-    Int_t *idtmed = fIdtmed->GetArray() - 1599 ;
+    idtmed = fIdtmed->GetArray() - 1599 ;
 
-    Int_t idrotm = 1;
-    AliMatrix(idrotm, 90.0, 0., 90.0, 90.0, 0.0, 0.0) ;
+    idrotm = 1;
+    //  gMC->Matrix(nmat, theta1, phi1, theta2, phi2, theta3, phi3) - see AliModule
+    AliMatrix(idrotm, 90.0, 0., 90.0, 90.0, 0.0, 0.0) ; 
 
     // Create the EMCAL Mother Volume (a polygone) within which to place the Detector and named XEN1 
 
     Float_t envelopA[10];
-    envelopA[0] = geom->GetArm1PhiMin();                         // minimum phi angle
-    envelopA[1] = geom->GetArm1PhiMax() - geom->GetArm1PhiMin(); // angular range in phi
-    envelopA[2] = geom->GetNPhi();                               // number of sections in phi
-    envelopA[3] = 2;                                             // 2 z coordinates
-    envelopA[4] = geom->ZFromEtaR(geom->GetEnvelop(1),
+    if(gn.Contains("TRD2")) { // TUBS
+       envelopA[0] = geom->GetEnvelop(0) - 10.; // rmin 
+       envelopA[1] = geom->GetEnvelop(1) + 12.; // rmax
+       //       envelopA[2] = geom->ZFromEtaR(geom->GetEnvelop(1), geom->GetArm1EtaMin());
+       envelopA[2] = 390.; // 6-feb-05
+       envelopA[3] = geom->GetArm1PhiMin();
+       envelopA[4] = geom->GetArm1PhiMax();
+       gMC->Gsvolu("XEN1", "TUBS", idtmed[1599], envelopA, 5) ;   // Tubs filled with air 
+       ENVELOP1.Set(5, envelopA);
+    // Position the EMCAL Mother Volume (XEN1) in Alice (ALIC)  
+       gMC->Gspos("XEN1", 1, "ALIC", 0.0, 0.0, 0.0, idrotm, "ONLY") ;
+    } else if(gn.Contains("TRD1") && gn.Contains("WSUC") ) { // TRD1 for WSUC facility
+      // 17-may-05 - just BOX
+      envelopA[0] = 26;
+      envelopA[1] = 15;
+      envelopA[2] = 30;
+      gMC->Gsvolu("XEN1", "BOX", idtmed[idSC], envelopA, 3) ;
+      ENVELOP1.Set(3);
+      for(int i=0; i<3; i++) ENVELOP1[i] = envelopA[i]; // 23-may-05  
+      // Position the EMCAL Mother Volume (XEN1) in WSUC  
+      gMC->Gspos("XEN1", 1, "WSUC", 0.0, 0.0, 0.0, idrotm, "ONLY") ;
+    } else { 
+      envelopA[0] = geom->GetArm1PhiMin();                         // minimum phi angle
+      envelopA[1] = geom->GetArm1PhiMax() - geom->GetArm1PhiMin(); // angular range in phi
+      envelopA[2] = geom->GetNPhi();                               // number of sections in phi
+      envelopA[3] = 2;                                             // 2 z coordinates
+      envelopA[4] = geom->ZFromEtaR(geom->GetEnvelop(1),
 				   geom->GetArm1EtaMin());       // z coordinate 1
     //add some padding for mother volume
-    envelopA[5] = geom->GetEnvelop(0) ;                          // rmin at z1
-    envelopA[6] = geom->GetEnvelop(1) ;                          // rmax at z1
-    envelopA[7] = geom->ZFromEtaR(geom->GetEnvelop(1),
+      envelopA[5] = geom->GetEnvelop(0) ;                          // rmin at z1
+      envelopA[6] = geom->GetEnvelop(1) ;                          // rmax at z1
+      envelopA[7] = geom->ZFromEtaR(geom->GetEnvelop(1),
 				  geom->GetArm1EtaMax());        // z coordinate 2
-    envelopA[8] = envelopA[5] ;                                  // radii are the same.
-    envelopA[9] = envelopA[6] ;                                  // radii are the same.
+      envelopA[8] = envelopA[5] ;                                  // radii are the same.
+      envelopA[9] = envelopA[6] ;                                  // radii are the same.
 
-    gMC->Gsvolu("XEN1", "PGON ", idtmed[1599], envelopA, 10) ;   // Polygone filled with air 
+      if(gn.Contains("SHISH")) envelopA[2] = geom->GetNPhiSuperModule();
 
+      gMC->Gsvolu("XEN1", "PGON ", idtmed[1599], envelopA, 10) ;   // Polygone filled with air 
+      ENVELOP1.Set(10, envelopA);
+      if (gDebug==2) {
+        printf("CreateGeometry: XEN1 = %f, %f\n", envelopA[5], envelopA[6]); 
+        printf("CreateGeometry: XU0 = %f, %f\n", envelopA[5], envelopA[6]); 
+      }
     // Position the EMCAL Mother Volume (XEN1) in Alice (ALIC)  
+      gMC->Gspos("XEN1", 1, "ALIC", 0.0, 0.0, 0.0, idrotm, "ONLY") ;
+    }
+
+    if(gn.Contains("SHISH")){
+      // Compact and twist
+      CreateShishKebabGeometry();
+      return;
+    }
 
-    gMC->Gspos("XEN1", 1, "ALIC", 0.0, 0.0, 0.0, idrotm, "ONLY") ;
-    
     if (AliLog::GetGlobalDebugLevel()>=2) {
       printf("CreateGeometry: XEN1 = %f, %f\n", envelopA[5], envelopA[6]); 
       printf("CreateGeometry: XU0 = %f, %f\n", envelopA[5], envelopA[6]); 
     }
+
     // Create mini-envelopes which will contain the Tower scintillator-radiator
     
     TString label ;
@@ -343,3 +530,788 @@ void AliEMCALv0::Init(void)
     printf(message.Data() ) ; 
   }
 }
+
+// 24-aug-04 by PAI
+void AliEMCALv0::CreateShishKebabGeometry()
+{  // TWIST, TRD1 and TRD2 
+  AliEMCALGeometry * g = GetGeometry(); 
+  TString gn(g->GetName()); gn.ToUpper(); 
+  // see AliModule::fIdtmed
+  //  idtmed = fIdtmed->GetArray() - 1599 ; // see AliEMCAL::::CreateMaterials()
+  //  int idAIR=1599, idPB = 1600, idSC = 1601, idSTEEL = 1603;
+  //  idAL = 1602;
+  Double_t par[10], xpos=0., ypos=0., zpos=0.;
+
+  CreateSmod("XEN1"); // 18-may-05 
+
+  CreateEmod("SMOD","EMOD"); // 18-may-95
+
+  // Sensitive SC  (2x2 tiles)
+  double parSCM0[5], *dummy = 0, parTRAP[11];
+  Double_t trd1Angle = g->GetTrd1Angle()*TMath::DegToRad(), tanTmp = TMath::Tan(trd1Angle/2.);
+   if(!gn.Contains("TRD")) { // standard module
+    par[0] = (g->GetECPbRadThick()+g->GetECScintThick())*g->GetNECLayers()/2.;
+    par[1] = par[2] = g->GetPhiTileSize();   // Symetric case
+    gMC->Gsvolu("SCM0", "BOX", idtmed[idSC], par, 3); // 2x2 tiles
+    gMC->Gspos("SCM0", 1, "EMOD", 0., 0., 0., 0, "ONLY") ;
+    // Division to tile size
+    gMC->Gsdvn("SCM1","SCM0", g->GetNPHIdiv(), 2); // y-axis
+    gMC->Gsdvn("SCM2","SCM1", g->GetNETAdiv(), 3); // z-axis
+  // put LED to the SCM2 
+    par[0] = g->GetECPbRadThick()/2.;
+    par[1] = par[2] = g->GetPhiTileSize()/2.; // Symetric case
+    gMC->Gsvolu("PBTI", "BOX", idtmed[idPB], par, 3);
+
+    printf(" Pb tiles \n");
+    int nr=0;
+    ypos = zpos = 0.0;
+    xpos = -sampleWidth*g->GetNECLayers()/2. + g->GetECPbRadThick()/2.;
+    for(int ix=0; ixGetNECLayers(); 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(" 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; izGetNECLayers(); 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; iGetArm1PhiMin() + 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=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; izGetNZ(); 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; iyGetNPhi(); 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; izGetNZ(); 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; izGetNZ(); 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; iyGetNPhi(); 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; ixGetNPhi(); 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; iyGetNPhi(); 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; iyGetNPhi(); iy++)
+      }
+    } 
+  } else {
+    xpos = g->GetSteelFrontThickness()/2.;
+    for(int iz=0; izGetNZ(); iz++) {
+      zpos = -smodPar2 + g->GetEtaModuleSize()*(2*iz+1)/2.;
+      for(int iy=0; iyGetNPhi(); 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()"," 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; izGetNECLayers(); 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; izGetNECLayers(); 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; izGetNECLayers(); 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; izGetNECLayers(); 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;
+  }
+}
diff --git a/EMCAL/AliEMCALv0.h b/EMCAL/AliEMCALv0.h
index a7227092aa1..7be48204361 100644
--- a/EMCAL/AliEMCALv0.h
+++ b/EMCAL/AliEMCALv0.h
@@ -15,6 +15,8 @@
 // --- ROOT system ---
 
 class TFile;
+class TList;
+class TNode;
 
 // --- AliRoot header files ---
 #include "AliEMCAL.h"
@@ -35,6 +37,7 @@ class AliEMCALv0 : public AliEMCAL {
   virtual ~AliEMCALv0(){} 
 
   virtual void BuildGeometry();// creates the geometry for the ROOT display
+  TNode *BuildGeometryOfWSUC();  // WSUC - test environment
   virtual void CreateGeometry() ;// creates the geometry for GEANT
   virtual void   Init(void) ;                                       // does nothing
   virtual Int_t  IsVersion(void) const { 
@@ -51,6 +54,23 @@ class AliEMCALv0 : public AliEMCAL {
     Fatal("operator =", "not implemented") ;  
     return *this ; 
   }
+  // ShishKebab 
+  void CreateShishKebabGeometry();
+  void CreateSmod(const char* mother="XEN1");
+  void CreateEmod(const char* mother="SMOD", const char* child="EMOD");
+  // TRD1
+  void Trd1Tower3X3(const double parSCM0[5]);
+  void Trd1Tower4X4();
+  void PbInTrap(const double parTRAP[11], TString n);
+  // TRD2 - 1th design
+  void Scm0InTrd2(const AliEMCALGeometry * g, const Double_t parEMOD[5], Double_t parSCM0[5]);
+  void Division2X2InScm0(const AliEMCALGeometry * g, const Double_t parSCM0[5]);
+  void PbInTrapForTrd2(const double *parTRAP, TString name);
+  // TRD2 - 2th design
+  void PbmoInTrd2(const AliEMCALGeometry * g, const Double_t parEMOD[5], Double_t parPBMO[5]);
+  void Division2X2InPbmo(const AliEMCALGeometry * g, const Double_t parPBMO[5]);
+
+  TList *fShishKebabModules; //! list of modules for twist geometries
   
  protected:
 
diff --git a/EMCAL/AliEMCALv1.cxx b/EMCAL/AliEMCALv1.cxx
index 48e3ca0baef..be1a95b19dc 100644
--- a/EMCAL/AliEMCALv1.cxx
+++ b/EMCAL/AliEMCALv1.cxx
@@ -57,7 +57,7 @@ AliEMCALv1::AliEMCALv1(const char *name, const char *title):
     // Standard Creator.
 
     fHits= new TClonesArray("AliEMCALHit",1000);
-    gAlice->GetMCApp()->AddHitList(fHits);
+    //    gAlice->GetMCApp()->AddHitList(fHits); // 20-dec-04 - advice of Andreas
 
     fNhits = 0;
     fIshunt     =  2; // All hits are associated with particles entering the calorimeter
@@ -104,7 +104,7 @@ void AliEMCALv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t ipa
 	new((*fHits)[fNhits]) AliEMCALHit(*newHit);
 	fNhits++;
     } // end if
-
+    
     delete newHit;
 }
 //______________________________________________________________________
diff --git a/EMCAL/AliEMCALv1.h b/EMCAL/AliEMCALv1.h
index 18d28f19374..15bcab4418c 100644
--- a/EMCAL/AliEMCALv1.h
+++ b/EMCAL/AliEMCALv1.h
@@ -47,13 +47,14 @@ public:
     return *this;}
  
     
-private:
+protected:
+  // Marco advice - 16-jan-05
   Int_t fCurPrimary;
   Int_t fCurParent;
   Int_t fCurTrack;
   Float_t fTimeCut;       // Cut to remove the background from the ALICE system
 
-  ClassDef(AliEMCALv1,8)//Implementation of EMCAL manager class to produce hits in a Central Calorimeter 
+  ClassDef(AliEMCALv1,9)//Implementation of EMCAL manager class to produce hits in a Central Calorimeter 
     
 };
 
diff --git a/EMCAL/AliEMCALv2.cxx b/EMCAL/AliEMCALv2.cxx
new file mode 100644
index 00000000000..7234e0cd586
--- /dev/null
+++ b/EMCAL/AliEMCALv2.cxx
@@ -0,0 +1,569 @@
+/**************************************************************************
+ * Copyright(c) 1998-2004, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+//_________________________________________________________________________
+//*-- Implementation version v2 of EMCAL Manager class 
+//*-- An object of this class does not produce digits
+//*-- It is the one to use if you do want to produce outputs in TREEH 
+//*--                  
+//*-- Author: Sahal Yacoob (LBL /UCT)
+//*--       : Jennifer Klay (LBL)
+// This Class not stores information on all particles prior to EMCAL entry - in order to facilitate analysis.
+// This is done by setting fIShunt =2, and flagging all parents of particles entering the EMCAL.
+
+// 15/02/2002 .... Yves Schutz
+//  1. fSamplingFraction and fLayerToPreshowerRatio have been removed
+//  2. Timing signal is collected and added to hit
+
+// --- ROOT system ---
+#include "TParticle.h"
+#include "TVirtualMC.h"
+#include "TBrowser.h"
+#include "TH1.h"
+#include "TH2.h"
+
+// --- Standard library ---
+
+// --- AliRoot header files ---
+#include "AliEMCALv2.h"
+#include "AliEMCALHit.h"
+#include "AliEMCALGeometry.h"
+#include "AliRun.h"
+#include "AliHeader.h"
+#include "AliMC.h"
+#include "AliPoints.h"
+#include "AliEMCALGetter.h"
+// for TRD1,2 case
+
+ClassImp(AliEMCALv2)
+
+  //extern "C" void gdaxis_(float &x0, float &y0, float &z0, float &axsiz);
+
+
+//______________________________________________________________________
+AliEMCALv2::AliEMCALv2():AliEMCALv1(), fGeometry(0){
+  // ctor
+
+}
+
+//______________________________________________________________________
+AliEMCALv2::AliEMCALv2(const char *name, const char *title): AliEMCALv1(name,title) {
+    // Standard Creator.
+
+  //    fGeant3 = (TGeant3*)gMC;
+    fHits= new TClonesArray("AliEMCALHit",1000);
+    gAlice->GetMCApp()->AddHitList(fHits);
+
+    fNhits    = 0;
+    fIshunt   = 2; // All hits are associated with particles entering the calorimeter
+    fTimeCut  = 30e-09;
+
+    fGeometry = GetGeometry(); 
+    fHDe = 0;
+    //    if (gDebug>0){
+    if (1){
+      TH1::AddDirectory(0);
+      fHDe = new TH1F("fHDe","De in EMCAL", 1000, 0., 1.);
+      fHistograms->Add(fHDe);
+      TH1::AddDirectory(1);
+    }
+}
+
+//______________________________________________________________________
+AliEMCALv2::~AliEMCALv2(){
+    // dtor
+
+    if ( fHits) {
+	fHits->Delete();
+	delete fHits;
+	fHits = 0;
+    }
+}
+
+//______________________________________________________________________
+void AliEMCALv2::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t iparent, Float_t ienergy, 
+			Int_t id, Float_t * hits,Float_t * p){
+    // Add a hit to the hit list.
+    // An EMCAL hit is the sum of all hits in a tower section
+    //   originating from the same entering particle 
+    static Int_t hitCounter;
+    static AliEMCALHit *newHit, *curHit;
+    static Bool_t deja;
+
+    deja = kFALSE;
+
+    newHit = new AliEMCALHit(shunt, primary, tracknumber, iparent, ienergy, id, hits, p);
+    for ( hitCounter = fNhits-1; hitCounter >= 0 && !deja; hitCounter-- ) {
+	curHit = (AliEMCALHit*) (*fHits)[hitCounter];
+	// We add hits with the same tracknumber, while GEANT treats
+	// primaries succesively
+	if(curHit->GetPrimary() != primary) 
+	  break;
+	if( *curHit == *newHit ) {
+	    *curHit = *curHit + *newHit;
+	    deja = kTRUE;
+	    //            break; // 30-aug-04 by PAI 
+	} // end if
+    } // end for hitCounter
+    
+    if ( !deja ) {
+	new((*fHits)[fNhits]) AliEMCALHit(*newHit);
+	fNhits++;
+    }
+    //    printf(" fNhits %i \n", fNhits); 
+    delete newHit;
+}
+//______________________________________________________________________
+void AliEMCALv2::StepManager(void){
+  // Accumulates hits as long as the track stays in a tower
+
+  // position wrt MRS and energy deposited
+  static Float_t        xyzte[5]={0.,0.,0.,0.,0.};// position wrt MRS, time and energy deposited
+  static Float_t        pmom[4]={0.,0.,0.,0.};
+  static TLorentzVector pos; // Lorentz vector of the track current position.
+  static TLorentzVector mom; // Lorentz vector of the track current momentum.
+  static Float_t ienergy = 0;
+  static TString curVolName;
+  static int supModuleNumber, moduleNumber, yNumber, xNumber, absid;
+  static int keyGeom=0;
+  static char *vn = "SX"; // 15-mar-05
+  static int nSMOP[7]={1,3,5,7,9,11}; // 30-mar-05
+  static int nSMON[7]={2,4,6,8,10,12};
+  static Float_t depositedEnergy=0.0; 
+
+  if(keyGeom == 0) {
+    keyGeom = 2;
+    if(gMC->VolId("PBMO")==0 || gMC->VolId("WSUC")==1) {
+      vn      = "SCMX";   // old TRD2(TRD1) or WSUC
+      keyGeom = 1;
+    }    
+    printf("AliEMCALv2::StepManager():  keyGeom %i : Sensetive volume %s \n", 
+    keyGeom, vn); 
+    if(gMC->VolId("WSUC")==1) printf(" WSUC - cosmic ray stand geometry \n");
+  }
+  Int_t tracknumber =  gAlice->GetMCApp()->GetCurrentTrackNumber();
+
+  curVolName = gMC->CurrentVolName();
+  if(curVolName.Contains(vn)) { // We are in a scintillator layer
+    //    printf(" keyGeom %i : Sensetive volume %s (%s) \n", keyGeom, curVolName.Data(), vn); 
+    
+    if( ((depositedEnergy = gMC->Edep()) > 0.)  && (gMC->TrackTime() < fTimeCut)){// Track is inside a scintillator and deposits some energy
+      //       Info("StepManager "," entry %i DE %f",++ientry, depositedEnergy); // for testing
+       if (fCurPrimary==-1) 
+	fCurPrimary=gAlice->GetMCApp()->GetPrimary(tracknumber);
+
+      if (fCurParent==-1 || tracknumber != fCurTrack) {
+	// Check parentage
+	Int_t parent=tracknumber;
+	if (fCurParent != -1) {
+	  while (parent != fCurParent && parent != -1) {
+	    TParticle *part=gAlice->GetMCApp()->Particle(parent);
+	    parent=part->GetFirstMother();
+	  }
+	}
+	if (fCurParent==-1 || parent==-1) {
+	  Int_t parent=tracknumber;
+	  TParticle *part=gAlice->GetMCApp()->Particle(parent);
+	  while (parent != -1 && fGeometry->IsInEMCAL(part->Vx(),part->Vy(),part->Vz())) {
+	    parent=part->GetFirstMother();
+	    if (parent!=-1) 
+	      part=gAlice->GetMCApp()->Particle(parent);
+	  } 
+	  fCurParent=parent;
+	  if (fCurParent==-1)
+	    Error("StepManager","Cannot find parent");
+	  else {
+	    TParticle *part=gAlice->GetMCApp()->Particle(fCurParent);
+	    ienergy = part->Energy(); 
+	  }
+	  while (parent != -1) {
+	    part=gAlice->GetMCApp()->Particle(parent);
+	    part->SetBit(kKeepBit);
+	    parent=part->GetFirstMother();
+	  }
+	}
+	fCurTrack=tracknumber;
+      }    
+      gMC->TrackPosition(pos);
+      xyzte[0] = pos[0];
+      xyzte[1] = pos[1];
+      xyzte[2] = pos[2];
+      xyzte[3] = gMC->TrackTime() ;       
+      
+      gMC->TrackMomentum(mom);
+      pmom[0] = mom[0];
+      pmom[1] = mom[1];
+      pmom[2] = mom[2];
+      pmom[3] = mom[3];
+      
+      //      if(ientry%200 > 0) return; // testing
+      supModuleNumber = moduleNumber = yNumber = xNumber = absid = 0;
+      if(keyGeom >= 1) { // old style
+        gMC->CurrentVolOffID(4, supModuleNumber);
+        gMC->CurrentVolOffID(3, moduleNumber);
+        gMC->CurrentVolOffID(1, yNumber);
+        gMC->CurrentVolOffID(0, xNumber); // really x number now
+      } else {
+        gMC->CurrentVolOffID(5, supModuleNumber);
+        gMC->CurrentVolOffID(4, moduleNumber);
+        gMC->CurrentVolOffID(1, yNumber);
+        gMC->CurrentVolOffID(0, xNumber); // really x number now
+        if    (strcmp(gMC->CurrentVolOffName(5),"SMOP")==0) supModuleNumber = nSMOP[supModuleNumber-1];
+        else if(strcmp(gMC->CurrentVolOffName(5),"SMON")==0) supModuleNumber = nSMON[supModuleNumber-1];
+        else   assert(0); // something wrong
+      }
+      absid = fGeometry->GetAbsCellId(supModuleNumber, moduleNumber, yNumber, xNumber);
+    
+      if (absid == -1) Fatal("StepManager()", "Wrong id ") ;
+
+      Float_t lightYield =  depositedEnergy ;
+      // Apply Birk's law (copied from G3BIRK)
+
+      if (gMC->TrackCharge()!=0) { // Check
+	  Float_t BirkC1_mod = 0;
+	if (fBirkC0==1){ // Apply correction for higher charge states
+	  if (TMath::Abs(gMC->TrackCharge())>=2) BirkC1_mod=fBirkC1*7.2/12.6;
+	  else                                    BirkC1_mod=fBirkC1;
+	}
+
+	Float_t dedxcm;
+	if (gMC->TrackStep()>0)  dedxcm=1000.*gMC->Edep()/gMC->TrackStep();
+	else                     dedxcm=0;
+	lightYield=lightYield/(1.+BirkC1_mod*dedxcm+fBirkC2*dedxcm*dedxcm);
+      } 
+
+      // use sampling fraction to get original energy --HG
+      //      xyzte[4] = lightYield * fGeometry->GetSampling();
+      xyzte[4] = lightYield; // 15-dec-04
+        
+      if (gDebug == -2) 
+      printf("#sm %2i #m %3i #x %1i #z %1i -> absid %i : xyzte[4] = %f\n",
+      supModuleNumber,moduleNumber,yNumber,xNumber,absid, xyzte[4]);
+
+      AddHit(fIshunt, fCurPrimary,tracknumber, fCurParent, ienergy, absid,  xyzte, pmom);
+    } // there is deposited energy
+  }
+}
+
+void AliEMCALv2::FinishEvent()
+{ // 26-may-05
+  static double de=0.;
+  de = GetDepositEnergy(0);
+  if(fHDe) fHDe->Fill(de);
+}
+
+Double_t AliEMCALv2::GetDepositEnergy(int print)
+{ // 23-mar-05 - for testing
+  if(fHits == 0) return 0.;
+  AliEMCALHit  *hit=0;
+  Double_t de=0.;
+  for(int ih=0; ihGetEntries(); 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(" 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; iUncheckedAt(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 (hit_it()) ) {
+    if (map[hit->GetIparent()]==-99)
+      cout << "Remapping, found -99 for parent id " << hit->GetIparent() << ", " << map[hit->GetIparent()] << ", i_hit " << i_hit << endl;
+    hit->SetIparent(map[hit->GetIparent()]);
+    if (map[hit->GetPrimary()]==-99)
+      cout << "Remapping, found -99 for primary id " << hit->GetPrimary() << ", " << map[hit->GetPrimary()] << ", i_hit " << i_hit << endl;
+    hit->SetPrimary(map[hit->GetPrimary()]);
+    i_hit++;
+  }
+}
+
+void AliEMCALv2::FinishPrimary() {
+  fCurPrimary=-1;
+  fCurParent=-1;
+  fCurTrack=-1;
+}
+*/
diff --git a/EMCAL/AliEMCALv2.h b/EMCAL/AliEMCALv2.h
new file mode 100644
index 00000000000..f3e299248f1
--- /dev/null
+++ b/EMCAL/AliEMCALv2.h
@@ -0,0 +1,75 @@
+#ifndef ALIEMCALV2_H
+#define ALIEMCALV2_H
+/* Copyright(c) 1998-2004, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                      
+         */
+/* $Id$ */
+
+//_________________________________________________________________________
+// Implementation version v2 of EMCAL Manager class for Shish-Kebab case 
+//*--                  
+//*-- Author:  Aleksei Pavlinov
+
+// --- ROOT system ---
+class TClonesArray;
+class TLorentzVector;
+class TFile;
+class TH1F;
+
+class TBrowser;
+class AliEMCALGeometry;
+
+// --- AliRoot header files ---
+#include "AliEMCALv1.h"
+
+// for TRD2 case
+//#include "TGeant3.h"
+
+class AliEMCALv2 : public AliEMCALv1 {
+  
+public:
+
+  AliEMCALv2(void) ; 
+  AliEMCALv2(const char *name, const char *title="") ;
+  // cpy ctor: no implementation yet
+  // requested by the Coding Convention
+  AliEMCALv2(const AliEMCALv1 & emcal):AliEMCALv1(emcal) {
+    Fatal("cpy ctor", "not implemented") ;  }
+  virtual ~AliEMCALv2(void) ;
+  virtual void  AddHit( Int_t shunt, Int_t primary, Int_t track, Int_t iparent, Float_t ienergy,
+			Int_t id, Float_t *hits, Float_t *p);
+
+  virtual void StepManager(void) ;
+  virtual void FinishEvent();
+
+  // Gives the version number 
+  virtual Int_t  IsVersion(void) const {return 2;}
+  virtual const TString Version(void)const {return TString("v2");}
+  //  virtual void RemapTrackHitIDs(Int_t *map);
+  //virtual void FinishPrimary();
+  // virtual void SetTimeCut(Float_t tc){ fTimeCut = tc;}
+  // virtual Float_t GetTimeCut(){return fTimeCut;}
+  // assignement operator requested by coding convention but not needed  
+  AliEMCALv2 & operator = (const AliEMCALv1 & /*rvalue*/){
+    Fatal("operator =", "not implemented") ;  
+    return *this;}
+  // 23-mar-05
+  Double_t GetDepositEnergy(int print=1); // *MENU*
+  // 30-aug-04
+  virtual void Browse(TBrowser* b);
+  // drawing
+  void DrawCalorimeterCut(const char *name="SMOD", int axis=3, double dcut=1.); // *MENU*
+  void DrawSuperModuleCut(const char *name="EMOD", int axis=2, double dcut=0.03, int fill = 6);//  *MENU*
+  void DrawTowerCut(const char *name="SCMY", int axis=2, double dcut=0., int fill=1, char *optShad="on");   //  *MENU*
+  void DrawAlicWithHits(int mode=1);                            // *MENU*
+  void SetVolumeAttributes(const char *name="SCM0",const int seen=1, const int color=1, const int fill=1); // *MENU*
+  void TestIndexTransition(int pri=0, int idmax=0); // *MENU*
+
+  AliEMCALGeometry* fGeometry; //!
+  TH1F*             fHDe;      //!
+
+  ClassDef(AliEMCALv2,1)    //Implementation of EMCAL manager class to produce hits in a Shish-Kebab
+    
+};
+
+#endif // AliEMCALV2_H
diff --git a/EMCAL/EMCALLinkDef.h b/EMCAL/EMCALLinkDef.h
index 14aa33f2c9a..f6a0b55ba9b 100644
--- a/EMCAL/EMCALLinkDef.h
+++ b/EMCAL/EMCALLinkDef.h
@@ -16,6 +16,7 @@
 #pragma link C++ class AliEMCALGeometry+; 
 #pragma link C++ class AliEMCALv0+;
 #pragma link C++ class AliEMCALv1+;
+#pragma link C++ class AliEMCALv2+;
 #pragma link C++ class AliEMCALHit+;
 #pragma link C++ class AliEMCALDigit+;
 #pragma link C++ class AliEMCALTick+;
@@ -49,9 +50,18 @@
 #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
diff --git a/EMCAL/libEMCAL.pkg b/EMCAL/libEMCAL.pkg
index a2022e47fa7..c10ee904efb 100644
--- a/EMCAL/libEMCAL.pkg
+++ b/EMCAL/libEMCAL.pkg
@@ -5,6 +5,7 @@ AliEMCAL.cxx \
 AliEMCALGeometry.cxx \
 AliEMCALv0.cxx \
 AliEMCALv1.cxx \
+AliEMCALv2.cxx \
 AliEMCALHit.cxx \
 AliEMCALDigit.cxx \
 AliEMCALSDigitizer.cxx \
@@ -38,14 +39,18 @@ AliEMCALRecParticle.cxx \
 AliEMCALFastRecParticle.cxx \
 AliEMCALPID.cxx \
 AliEMCALPIDv1.cxx \
-AliEMCALTracker.cxx \
 AliEMCALLoader.cxx \
 AliEMCALReconstructor.cxx \
-AliEMCALRawStream.cxx
+AliEMCALRawStream.cxx \
+AliEMCALGeneratorFactory.cxx \
+AliEMCALShishKebabModule.cxx\
+AliEMCALShishKebabTrd1Module.cxx\
+AliEMCALGeometryOfflineTrd1.cxx\
+AliEMCALWsuCosmicRaySetUp.cxx
 
 HDRS= $(SRCS:.cxx=.h) 
 
 DHDR:=EMCALLinkDef.h
 
-EINCLUDE:=PYTHIA6 RAW
+EINCLUDE:=PYTHIA6 RAW EVGEN THijing
 
-- 
2.39.3