EMCAL/AliEMCALReconstructor.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //_________________________________________________________________________
  19 //--
  20 //-- Yves Schutz (SUBATECH)
  21 // Reconstruction class. Redesigned from the old AliReconstructionner class and
  22 // derived from STEER/AliReconstructor.
  23 //
  24 // --- ROOT system ---
  25
  26 // --- Standard library ---
  27
  28 // --- AliRoot header files ---
  29 #include "AliEMCALReconstructor.h"
  30
  31 #include "AliESDEvent.h"
  32 #include "AliESDCaloCluster.h"
  33 #include "AliESDtrack.h"
  34 #include "AliEMCALLoader.h"
  35 #include "AliEMCALRawUtils.h"
  36 #include "AliEMCALClusterizerv1.h"
  37 #include "AliEMCALRecPoint.h"
  38 #include "AliEMCALPID.h"
  39 #include "AliEMCALTrigger.h"
  40 #include "AliRawReader.h"
  41 #include "AliCDBEntry.h"
  42 #include "AliCDBManager.h"
  43 #include "AliEMCALRecParam.h"
  44 // to be removed - it is here just because of geom
  45 #include "AliRun.h"
  46 #include "AliRunLoader.h"
  47
  48 ClassImp(AliEMCALReconstructor)
  49
  50 AliEMCALRecParam* AliEMCALReconstructor::fgkRecParam = 0;  // EMCAL rec. parameters
  51
  52 //____________________________________________________________________________
  53 AliEMCALReconstructor::AliEMCALReconstructor()
  54   : fDebug(kFALSE)
  55 {
  56   // ctor
  57 }
  58
  59 //____________________________________________________________________________
  60 AliEMCALReconstructor::AliEMCALReconstructor(const AliEMCALReconstructor & rec)
  61   : AliReconstructor(rec),
  62     fDebug(rec.fDebug)
  63 {
  64   //copy ctor
  65 }
  66
  67 //____________________________________________________________________________
  68 AliEMCALReconstructor::~AliEMCALReconstructor()
  69 {
  70   // dtor
  71 }
  72
  73 //____________________________________________________________________________
  74 void AliEMCALReconstructor::InitRecParam() const
  75 {
  76   // Check if the instance of AliEMCALRecParam exists,
  77   // if not, get it from OCDB if available, otherwise create a default one
  78
  79  if (!fgkRecParam  && (AliCDBManager::Instance()->IsDefaultStorageSet())) {
  80     AliCDBEntry *entry = (AliCDBEntry*)
  81       AliCDBManager::Instance()->Get("EMCAL/Config/RecParam");
  82     if (entry) fgkRecParam =  (AliEMCALRecParam*) entry->GetObject();
  83   }
  84
  85   if(!fgkRecParam){
  86     AliWarning("The Reconstruction parameters for EMCAL nonitialized - Used default one");
  87     fgkRecParam = new AliEMCALRecParam;
  88   }
  89 }
  90
  91 //____________________________________________________________________________
  92 void AliEMCALReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
  93 {
  94   // method called by AliReconstruction;
  95   // Only the clusterization is performed,; the rest of the reconstruction is done in FillESD because the track
  96   // segment maker needs access to the AliESD object to retrieve the tracks reconstructed by
  97   // the global tracking.
  98   // Works on the current event.
  99
 100   InitRecParam();
 101   AliEMCALClusterizerv1 clu;
 102   clu.SetInput(digitsTree);
 103   clu.SetOutput(clustersTree);
 104   if ( Debug() )
 105     clu.Digits2Clusters("deb all") ;
 106   else
 107     clu.Digits2Clusters("pseudo") ;
 108 }
 109
 110 //____________________________________________________________________________
 111 void AliEMCALReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const
 112
 113 {
 114   // Conversion from raw data to
 115   // EMCAL digits.
 116   // Works on a single-event basis
 117
 118   rawReader->Reset() ;
 119
 120   TClonesArray *digitsArr = new TClonesArray("AliEMCALDigit",100);
 121   Int_t bufsize = 32000;
 122   digitsTree->Branch("EMCAL", &digitsArr, bufsize);
 123
 124   static AliEMCALRawUtils rawUtils;
 125   rawUtils.Raw2Digits(rawReader,digitsArr);
 126
 127   digitsTree->Fill();
 128   digitsArr->Delete();
 129   delete digitsArr;
 130
 131 }
 132
 133 //____________________________________________________________________________
 134 void AliEMCALReconstructor::FillESD(TTree* /*digitsTree*/, TTree* clustersTree,
 135                                     AliESDEvent* esd) const
 136 {
 137   // Called by AliReconstruct after Reconstruct() and global tracking and vertexing
 138   // Works on the current event
 139   const double timeScale = 1.e+11; // transition constant from sec to 0.01 ns
 140
 141   // Creates AliESDCaloCluster from AliEMCALRecPoints
 142
 143   TObjArray *clusters = new TObjArray(100);
 144   TBranch *branch = clustersTree->GetBranch("EMCALECARP");
 145   branch->SetAddress(&clusters);
 146   clustersTree->GetEvent(0);
 147
 148   Int_t nClusters = clusters->GetEntries(), nClustersNew=0;
 149   AliDebug(1,Form("%d clusters",nClusters));
 150   //  Int_t nRP=0, nPC=0; // in input
 151   esd->SetFirstEMCALCluster(esd->GetNumberOfCaloClusters()); // Put after Phos clusters
 152
 153   //######################################################
 154   //#########Calculate trigger and set trigger info###########
 155   //######################################################
 156
 157   AliEMCALTrigger tr ;
 158   //   tr.SetPatchSize(1);//create 4x4 patches
 159   tr.Trigger();
 160
 161   Float_t maxAmp2x2  = tr.Get2x2MaxAmplitude();
 162   Float_t maxAmpnxn  = tr.GetnxnMaxAmplitude();
 163   Float_t ampOutOfPatch2x2  = tr.Get2x2AmpOutOfPatch() ;
 164   Float_t ampOutOfPatchnxn  = tr.GetnxnAmpOutOfPatch() ;
 165
 166   AliEMCALGeometry * geom = 0;
 167   AliRunLoader *runLoader = AliRunLoader::GetRunLoader();
 168   if (runLoader->GetAliRun() && runLoader->GetAliRun()->GetDetector("EMCAL"))
 169     geom = dynamic_cast<AliEMCAL*>(runLoader->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
 170   if (geom == 0)
 171     geom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaulGeometryName());
 172
 173   Int_t iSM2x2      = tr.Get2x2SuperModule();
 174   Int_t iSMnxn      = tr.GetnxnSuperModule();
 175   Int_t iCellPhi2x2 = tr.Get2x2CellPhi();
 176   Int_t iCellPhinxn = tr.GetnxnCellPhi();
 177   Int_t iCellEta2x2 = tr.Get2x2CellEta();
 178   Int_t iCellEtanxn = tr.GetnxnCellEta();
 179
 180   AliDebug(2, Form("Trigger 2x2 max amp %f, out amp %f, SM %d, iphi %d ieta %d",  maxAmp2x2, ampOutOfPatch2x2, iSM2x2,iCellPhi2x2, iCellEta2x2));
 181   AliDebug(2, Form("Trigger 4x4 max amp %f , out amp %f, SM %d, iphi %d, ieta %d",  maxAmpnxn, ampOutOfPatchnxn, iSMnxn,iCellPhinxn, iCellEtanxn));
 182
 183   TVector3    pos2x2(-1,-1,-1);
 184   TVector3    posnxn(-1,-1,-1);
 185
 186   Int_t iAbsId2x2 = geom->GetAbsCellIdFromCellIndexes( iSM2x2, iCellPhi2x2, iCellEta2x2) ;
 187   Int_t iAbsIdnxn = geom->GetAbsCellIdFromCellIndexes( iSMnxn, iCellPhinxn, iCellEtanxn) ;
 188   geom->GetGlobal(iAbsId2x2, pos2x2);
 189   geom->GetGlobal(iAbsIdnxn, posnxn);
 190
 191   TArrayF triggerPosition(6);
 192   triggerPosition[0] = pos2x2(0) ;
 193   triggerPosition[1] = pos2x2(1) ;
 194   triggerPosition[2] = pos2x2(2) ;
 195   triggerPosition[3] = posnxn(0) ;
 196   triggerPosition[4] = posnxn(1) ;
 197   triggerPosition[5] = posnxn(2) ;
 198
 199   TArrayF triggerAmplitudes(4);
 200   triggerAmplitudes[0] = maxAmp2x2 ;
 201   triggerAmplitudes[1] = ampOutOfPatch2x2 ;
 202   triggerAmplitudes[2] = maxAmpnxn ;
 203   triggerAmplitudes[3] = ampOutOfPatchnxn ;
 204
 205   esd->AddEMCALTriggerPosition(triggerPosition);
 206   esd->AddEMCALTriggerAmplitudes(triggerAmplitudes);
 207
 208   //######################################################
 209   //#######################TRACK MATCHING###############
 210   //######################################################
 211   //Fill list of integers, each one is index of track to which the cluster belongs.
 212
 213   // step 1 - initialize array of matched track indexes
 214   Int_t *matchedTrack = new Int_t[nClusters];
 215   for (Int_t iclus = 0; iclus < nClusters; iclus++)
 216     matchedTrack[iclus] = -1;  // neg. index --> no matched track
 217
 218   // step 2, change the flag for all matched clusters found in tracks
 219   Int_t iemcalMatch = -1;
 220   Int_t endtpc = esd->GetNumberOfTracks();
 221   for (Int_t itrack = 0; itrack < endtpc; itrack++) {
 222     AliESDtrack * track = esd->GetTrack(itrack) ; // retrieve track
 223     iemcalMatch = track->GetEMCALcluster();
 224     if(iemcalMatch >= 0) matchedTrack[iemcalMatch] = itrack;
 225   }
 226
 227   //########################################
 228   //##############Fill CaloClusters#############
 229   //########################################
 230
 231
 232   for (Int_t iClust = 0 ; iClust < nClusters ; iClust++) {
 233     const AliEMCALRecPoint * clust = (const AliEMCALRecPoint*)clusters->At(iClust);
 234     //if(clust->GetClusterType()== AliESDCaloCluster::kEMCALClusterv1) nRP++; else nPC++;
 235     if (Debug()) clust->Print();
 236     // Get information from EMCAL reconstruction points
 237     Float_t xyz[3];
 238     TVector3 gpos;
 239     clust->GetGlobalPosition(gpos);
 240     for (Int_t ixyz=0; ixyz<3; ixyz++)
 241       xyz[ixyz] = gpos[ixyz];
 242
 243     Int_t digitMult = clust->GetMultiplicity();
 244     TArrayS amplList(digitMult);
 245     TArrayS timeList(digitMult);
 246     TArrayS digiList(digitMult);
 247     Float_t *amplFloat = clust->GetEnergiesList();
 248     Float_t *timeFloat = clust->GetTimeList();
 249     Int_t   *digitInts = clust->GetAbsId();
 250     Float_t elipAxis[2];
 251     clust->GetElipsAxis(elipAxis);
 252
 253     // Convert Float_t* and Int_t* to Short_t* to save memory
 254     // Problem : we should recalculate a cluster characteristics when discard digit(s)
 255     Int_t newdigitMult = 0;
 256     for (Int_t iDigit=0; iDigit<digitMult; iDigit++) {
 257       if (amplFloat[iDigit] > 0) {
 258         amplList[newdigitMult] = (UShort_t)(amplFloat[iDigit]*500);
 259         // Time in units of 0.01 ns = 10 ps
 260         if(timeFloat[iDigit] < 65536./timeScale)
 261           timeList[newdigitMult] = (UShort_t)(timeFloat[iDigit]*timeScale);
 262         else
 263           timeList[newdigitMult] = 65535;
 264         digiList[newdigitMult] = (UShort_t)(digitInts[iDigit]);
 265         newdigitMult++;
 266       }
 267       else if (clust->GetClusterType() != AliESDCaloCluster::kEMCALPseudoCluster)
 268         Warning("FillESD()","Negative or 0 digit amplitude in cluster");
 269     }
 270
 271     if(newdigitMult > 0) { // accept cluster if it has some digit
 272       nClustersNew++;
 273
 274       amplList.Set(newdigitMult);
 275       timeList.Set(newdigitMult);
 276       digiList.Set(newdigitMult);
 277
 278       //Primaries
 279       Int_t  parentMult  = 0;
 280       Int_t *parentList =  clust->GetParents(parentMult);
 281
 282       // fills the ESDCaloCluster
 283       AliESDCaloCluster * ec = new AliESDCaloCluster() ;
 284       ec->SetClusterType(clust->GetClusterType());
 285       ec->SetPosition(xyz);
 286       ec->SetE(clust->GetEnergy());
 287       ec->AddDigitAmplitude(amplList);
 288       ec->AddDigitTime(timeList);
 289       ec->AddDigitIndex(digiList);
 290
 291       if(clust->GetClusterType()== AliESDCaloCluster::kEMCALClusterv1){
 292
 293         ec->SetClusterDisp(clust->GetDispersion());
 294         ec->SetClusterChi2(-1); //not yet implemented
 295         ec->SetM02(elipAxis[0]*elipAxis[0]) ;
 296         ec->SetM20(elipAxis[1]*elipAxis[1]) ;
 297         ec->SetM11(-1) ;        //not yet implemented
 298
 299        TArrayI arrayTrackMatched(1);// Only one track, temporal solution.
 300        arrayTrackMatched[0]= matchedTrack[iClust];
 301        ec->AddTracksMatched(arrayTrackMatched);
 302
 303        TArrayI arrayParents(parentMult,parentList);
 304        ec->AddLabels(arrayParents);
 305       }
 306
 307       // add the cluster to the esd object
 308       esd->AddCaloCluster(ec);
 309       delete ec;
 310     }
 311
 312   } // cycle on clusters
 313
 314   delete [] matchedTrack;
 315
 316   esd->SetNumberOfEMCALClusters(nClustersNew);
 317   //if(nClustersNew != nClusters)
 318   //printf(" ##### nClusters %i -> new %i ##### \n", nClusters, nClustersNew );
 319
 320   //Fill ESDCaloCluster with PID weights
 321   AliEMCALPID *pid = new AliEMCALPID;
 322   //pid->SetPrintInfo(kTRUE);
 323   pid->SetReconstructor(kTRUE);
 324   pid->RunPID(esd);
 325   delete pid;
 326 }
 327
 328