]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALv2.cxx
from Per Thomas: files for generation of PeakFinder OCDB object
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALv2.cxx
index ebc5284540db4ed51451122c44c9dcc01b46414a..0a87fcb44756dd31d4aa32c43643ee70074fe838 100644 (file)
@@ -20,7 +20,7 @@
 //*-- 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 : Aleksei Pavlinov (WSU)
+//*-- Author : Alexei Pavlinov (WSU)
 
 // 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.
@@ -44,7 +44,7 @@
 #include "AliHeader.h"
 #include "AliMC.h"
 #include "AliStack.h"
-#include "AliPoints.h"
+#include "AliTrackReference.h"
 // for TRD1 case only; May 31,2006
 
 ClassImp(AliEMCALv2)
@@ -62,7 +62,7 @@ AliEMCALv2::AliEMCALv2(const char *name, const char *title)
 {
     // Standard Creator.
 
-    fHits= new TClonesArray("AliEMCALHit",1000);
+    //fHits= new TClonesArray("AliEMCALHit",1000); //Already done in ctor of v1
     gAlice->GetMCApp()->AddHitList(fHits);
 
     fNhits    = 0;
@@ -76,11 +76,12 @@ AliEMCALv2::AliEMCALv2(const char *name, const char *title)
 AliEMCALv2::~AliEMCALv2(){
     // dtor
 
-    if ( fHits) {
-       fHits->Delete();
-       delete fHits;
-       fHits = 0;
-    }
+  //Already done in dtor of v1
+//  if ( fHits ) {
+//    fHits->Clear();
+//    delete fHits;
+//    fHits = 0;
+//  }
 }
 
 //______________________________________________________________________
@@ -91,8 +92,8 @@ void AliEMCALv2::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t ipa
     //   originating from the same entering particle 
     static Int_t hitCounter;
     static AliEMCALHit *newHit, *curHit;
-    static Bool_t deja;
-
+    static Bool_t deja ;
+  
     deja = kFALSE;
 
     newHit = new AliEMCALHit(shunt, primary, tracknumber, iparent, ienergy, id, hits, p);
@@ -127,10 +128,10 @@ void AliEMCALv2::StepManager(void){
   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; // part->Energy();
-  static TString curVolName;
-  static int supModuleNumber, moduleNumber, yNumber, xNumber, absid;
+  static TString curVolName="";
+  static int supModuleNumber=-1, moduleNumber=-1, yNumber=-1, xNumber=-1, absid=-1;
   static int keyGeom=1;  //real TRD1 geometry
-  static char *vn = "SCMX"; // Apr 13, 2006 - only TRD1 case now
+  static const char *vn = "SCMX"; // Apr 13, 2006 - only TRD1 case now
   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; 
@@ -146,8 +147,8 @@ void AliEMCALv2::StepManager(void){
     if(gMC->VolId("WSUC")==1) printf(" WSUC - cosmic ray stand geometry \n");
   }
   Int_t tracknumber =  gAlice->GetMCApp()->GetCurrentTrackNumber();
-  Int_t parent;
-  TParticle* part;
+  Int_t parent=0;
+  TParticle* part=0;
 
   curVolName = gMC->CurrentVolName();
   if(curVolName.Contains(vn) || curVolName.Contains("SCX")) { // We are in a scintillator layer; SCX for 3X3
@@ -159,8 +160,8 @@ void AliEMCALv2::StepManager(void){
 
       if (fCurParent==-1 || tracknumber != fCurTrack) {
        // Check parentage
-       //Int_t parent=tracknumber;
        parent=tracknumber;
+
        if (fCurParent != -1) {
          while (parent != fCurParent && parent != -1) {
            //TParticle *part=gAlice->GetMCApp()->Particle(parent);
@@ -185,6 +186,10 @@ void AliEMCALv2::StepManager(void){
            //TParticle *part=gAlice->GetMCApp()->Particle(fCurParent);
            part=gAlice->GetMCApp()->Particle(fCurParent);
            ienergy = part->Energy(); 
+
+           //Add reference to parent in TR tree.       
+           AddTrackReference(tracknumber, AliTrackReference::kEMCAL);
+
          }
          while (parent != -1) {
            part=gAlice->GetMCApp()->Particle(parent);
@@ -219,7 +224,7 @@ void AliEMCALv2::StepManager(void){
           if     (strcmp(gMC->CurrentVolOffName(0),"SCX1")==0) xNumber=1;
           else if(strcmp(gMC->CurrentVolOffName(0),"SCX2")==0) xNumber=2;
           else if(strcmp(gMC->CurrentVolOffName(0),"SCX3")==0) xNumber=3;
-          else Fatal("StepManager()", "Wrong name of sensetive volume in 3X3 case : %s ", gMC->CurrentVolOffName(0));
+          else Fatal("StepManager()", "Wrong name of sensitive volume in 3X3 case : %s ", gMC->CurrentVolOffName(0));
        }
       } else {
         gMC->CurrentVolOffID(5, supModuleNumber);
@@ -230,11 +235,30 @@ void AliEMCALv2::StepManager(void){
         else if(strcmp(gMC->CurrentVolOffName(5),"SMON")==0) supModuleNumber = nSMON[supModuleNumber-1];
         else   assert(0); // something wrong
       }
-      absid = fGeometry->GetAbsCellId(supModuleNumber-1, moduleNumber-1, yNumber-1, xNumber-1);
-    
+               
+      // Due to problem with index ordering conventions the calcultation of absid is no more like this:        
+      //absid = fGeometry->GetAbsCellId(smNumber, moduleNumber-1, yNumber-1, xNumber-1);
+      
+      //Swap A side in Phi and C side in Eta due to wrong indexing.
+      Int_t iphi = -1;
+      Int_t ieta = -1;
+      Int_t smNumber = supModuleNumber-1;
+      Int_t smType   = 1;
+      fGeometry->GetCellPhiEtaIndexInSModule(smNumber,moduleNumber-1,yNumber-1,xNumber-1, iphi, ieta);
+      if (smNumber%2 == 0) {
+       ieta = ((fGeometry->GetCentersOfCellsEtaDir()).GetSize()-1)-ieta;// 47-ieta, revert the ordering on A side in order to keep convention.
+      }
+      else {  
+       if(smNumber >= 10) smType = 2 ; //half supermodule
+       iphi= ((fGeometry->GetCentersOfCellsPhiDir()).GetSize()/smType-1)-iphi;//23-iphi, revert the ordering on C side in order to keep convention.
+      }
+      
+      //Once we know the indexes, calculate the absolute ID
+      absid = fGeometry->GetAbsCellIdFromCellIndexes(smNumber, iphi, ieta);
+      
       if (absid < 0) {
         printf(" supModuleNumber %i : moduleNumber %i : yNumber %i : xNumber %i \n",
-        supModuleNumber, moduleNumber, yNumber, xNumber); 
+              supModuleNumber, moduleNumber, yNumber, xNumber); 
        Fatal("StepManager()", "Wrong id : %i ", absid) ; 
       }
 
@@ -248,7 +272,7 @@ void AliEMCALv2::StepManager(void){
          else                                   birkC1Mod = fBirkC1;
        }
 
-       Float_t dedxcm;
+       Float_t dedxcm=0.;
        if (gMC->TrackStep()>0)  dedxcm=1000.*gMC->Edep()/gMC->TrackStep();
        else                     dedxcm=0;
        lightYield=lightYield/(1.+birkC1Mod*dedxcm+fBirkC2*dedxcm*dedxcm);
@@ -277,7 +301,7 @@ void AliEMCALv2::Browse(TBrowser* 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
-  TString g(fGeometry->GetName());
+  TString g(GetGeometry()->GetName());
   g.ToUpper();
   gMC->Gsatt("*", "seen", 0);
 
@@ -295,7 +319,7 @@ void AliEMCALv2::DrawCalorimeterCut(const char *name, int axis, double dcut)
   st += ", zcut, ";
   st += name;
 
-  char *optShad = "on", *optHide="on";
+  const char *optShad = "on", *optHide="on";
   double cxy=0.02;
   if     (axis==1) {
     dcut = 0.;
@@ -308,11 +332,12 @@ void AliEMCALv2::DrawCalorimeterCut(const char *name, int axis, double dcut)
   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);
+  const Int_t buffersize = 255;
+  char cmd[buffersize];
+  snprintf(cmd,buffersize,"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());
+  snprintf(cmd,buffersize,"gMC->Gdhead(1111, \"%s\")\n", st.Data());
   printf("%s\n",cmd); gROOT->ProcessLine(cmd);
 }
 
@@ -322,7 +347,7 @@ void AliEMCALv2::DrawSuperModuleCut(const char *name, int axis, double dcut, int
  // 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"};
+  const char *tit[3]={"xcut", "ycut", "zcut"};
   if(axis<1) axis=1; if(axis>3) axis=3;
 
   gMC->Gsatt("*", "seen", 0);
@@ -332,7 +357,7 @@ void AliEMCALv2::DrawSuperModuleCut(const char *name, int axis, double dcut, int
   SetVolumeAttributes(name, 1, 5, fill);    // yellow 
 
   double cxy=0.055, x0=10., y0=10.;
-  char *optShad = "on", *optHide="on";
+  const char *optShad = "on", *optHide="on";
   SetVolumeAttributes("STPL", 1, 3, fill);  // green 
   if     (axis==1) {
     gMC->Gsatt("STPL", "seen", 0);
@@ -348,11 +373,12 @@ void AliEMCALv2::DrawSuperModuleCut(const char *name, int axis, double dcut, int
   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);
+  const Int_t buffersize = 255;
+  char cmd[buffersize];
+  snprintf(cmd,buffersize,"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());
+  snprintf(cmd,buffersize,"gMC->Gdhead(1111, \"%s\")\n", st.Data());
   printf("%s\n",cmd); gROOT->ProcessLine(cmd);
   // hint for testing
   printf("Begin of super module\n");
@@ -364,7 +390,7 @@ void AliEMCALv2::DrawSuperModuleCut(const char *name, int axis, double dcut, int
 }
 
 //___________________________________________________________
-void AliEMCALv2::DrawTowerCut(const char *name, int axis, double dcut, int fill, char *optShad)
+void AliEMCALv2::DrawTowerCut(const char *name, int axis, double dcut, int fill, const 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;
@@ -400,15 +426,16 @@ void AliEMCALv2::DrawTowerCut(const char *name, int axis, double dcut, int fill,
   if (mn.Contains("EMOD")) {
     cx = cy = 0.5;
   }
-  char cmd[200];
+  const Int_t buffersize = 255;
+  char cmd[buffersize];
 
   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);
+  snprintf(cmd,buffersize,"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());
+  snprintf(cmd,buffersize,"gMC->Gdhead(1111, \"%s\")\n", st.Data());
   printf("%s\n",cmd); gROOT->ProcessLine(cmd);
 }
 
@@ -431,7 +458,7 @@ void AliEMCALv2::DrawAlicWithHits(int mode)
   TClonesArray *hits = Hits();
   Int_t nhits = hits->GetEntries(), absId, nSupMod, nModule, nIphi, nIeta, iphi, ieta;
   AliEMCALHit *hit = 0;
-  Double_t de, des=0.;
+  Double_t de=0., des=0.;
   for(Int_t i=0; i<nhits; i++) {
     hit   = (AliEMCALHit*)hits->UncheckedAt(i);
     absId = hit->GetId();