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