PHOS/AliPHOSReconstructor.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 "AliLog.h"
  30 #include "AliAltroMapping.h"
  31 #include "AliESDEvent.h"
  32 #include "AliESDCaloCluster.h"
  33 #include "AliPHOSReconstructor.h"
  34 #include "AliPHOSClusterizerv1.h"
  35 #include "AliPHOSTrackSegmentMakerv1.h"
  36 #include "AliPHOSPIDv1.h"
  37 #include "AliPHOSTracker.h"
  38 #include "AliRawReader.h"
  39 #include "AliPHOSTrigger.h"
  40 #include "AliPHOSGeometry.h"
  41 #include "AliPHOSRecoParam.h"
  42 #include "AliPHOSRecoParamEmc.h"
  43 #include "AliPHOSRecoParamCpv.h"
  44 #include "AliPHOSDigit.h"
  45 #include "AliPHOSTrackSegment.h"
  46 #include "AliPHOSEmcRecPoint.h"
  47 #include "AliPHOSRecParticle.h"
  48 #include "AliPHOSRawDecoder.h"
  49 #include "AliPHOSRawDecoderv1.h"
  50 #include "AliPHOSRawDigiProducer.h"
  51 #include "AliPHOSPulseGenerator.h"
  52
  53 ClassImp(AliPHOSReconstructor)
  54
  55 Bool_t AliPHOSReconstructor::fgDebug = kFALSE ;
  56 AliPHOSRecoParam* AliPHOSReconstructor::fgkRecoParamEmc =0;  // EMC rec. parameters
  57 AliPHOSRecoParam* AliPHOSReconstructor::fgkRecoParamCpv =0;  // CPV rec. parameters
  58
  59 //____________________________________________________________________________
  60 AliPHOSReconstructor::AliPHOSReconstructor() :
  61   fGeom(NULL)
  62 {
  63   // ctor
  64
  65   if (!fgkRecoParamEmc) {
  66     AliWarning("The Reconstruction parameters for EMC nonitialized - Used default one");
  67     fgkRecoParamEmc = AliPHOSRecoParamEmc::GetEmcDefaultParameters();
  68   }
  69
  70   if (!fgkRecoParamCpv) {
  71     AliWarning("The Reconstruction parameters for CPV nonitialized - Used default one");
  72     fgkRecoParamCpv = AliPHOSRecoParamCpv::GetCpvDefaultParameters();
  73   }
  74
  75   fGeom = AliPHOSGeometry::GetInstance("IHEP","");
  76 }
  77
  78 //____________________________________________________________________________
  79   AliPHOSReconstructor::~AliPHOSReconstructor()
  80 {
  81   // dtor
  82   delete fGeom;
  83 }
  84
  85 //____________________________________________________________________________
  86 void AliPHOSReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
  87 {
  88   // 'single-event' local reco method called by AliReconstruction;
  89   // Only the clusterization is performed,; the rest of the reconstruction is done in FillESD because the track
  90   // segment maker needs access to the AliESDEvent object to retrieve the tracks reconstructed by
  91   // the global tracking.
  92
  93   AliPHOSClusterizerv1 clu(fGeom);
  94   clu.SetInput(digitsTree);
  95   clu.SetOutput(clustersTree);
  96   if ( Debug() )
  97     clu.Digits2Clusters("deb all") ;
  98   else
  99     clu.Digits2Clusters("") ;
 100 }
 101
 102 //____________________________________________________________________________
 103 void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
 104                                    AliESDEvent* esd) const
 105 {
 106   // This method produces PHOS rec-particles,
 107   // then it creates AliESDtracks out of them and
 108   // write tracks to the ESD
 109
 110   AliPHOSTrackSegmentMaker *tsm = new AliPHOSTrackSegmentMakerv1(fGeom);
 111   AliPHOSPID               *pid = new AliPHOSPIDv1              (fGeom);
 112
 113   // do current event; the loop over events is done by AliReconstruction::Run()
 114   tsm->SetESD(esd) ;
 115   tsm->SetInput(clustersTree);
 116   if ( Debug() )
 117     tsm->Clusters2TrackSegments("deb all") ;
 118   else
 119     tsm->Clusters2TrackSegments("") ;
 120
 121   pid->SetInput(clustersTree, tsm->GetTrackSegments()) ;
 122   pid->SetESD(esd) ;
 123   if ( Debug() )
 124     pid->TrackSegments2RecParticles("deb all") ;
 125   else
 126     pid->TrackSegments2RecParticles("") ;
 127
 128
 129   // This function creates AliESDtracks from AliPHOSRecParticles
 130   //         and
 131   // writes them to the ESD
 132
 133   TClonesArray *recParticles  = pid->GetRecParticles();
 134   Int_t nOfRecParticles = recParticles->GetEntries();
 135
 136   esd->SetNumberOfPHOSClusters(nOfRecParticles) ;
 137   esd->SetFirstPHOSCluster(esd->GetNumberOfCaloClusters()) ;
 138
 139   AliDebug(2,Form("%d rec. particles, option %s",nOfRecParticles,GetOption()));
 140
 141   //#########Calculate trigger and set trigger info###########
 142
 143   AliPHOSTrigger tr ;
 144   //   tr.SetPatchSize(1);//create 4x4 patches
 145   tr.Trigger();
 146
 147   Float_t maxAmp2x2  = tr.Get2x2MaxAmplitude();
 148   Float_t maxAmpnxn  = tr.GetnxnMaxAmplitude();
 149   Float_t ampOutOfPatch2x2  = tr.Get2x2AmpOutOfPatch() ;
 150   Float_t ampOutOfPatchnxn  = tr.GetnxnAmpOutOfPatch() ;
 151
 152   Int_t iSM2x2      = tr.Get2x2SuperModule();
 153   Int_t iSMnxn      = tr.GetnxnSuperModule();
 154   Int_t iCrystalPhi2x2 = tr.Get2x2CrystalPhi();
 155   Int_t iCrystalPhinxn = tr.GetnxnCrystalPhi();
 156   Int_t iCrystalEta2x2 = tr.Get2x2CrystalEta();
 157   Int_t iCrystalEtanxn = tr.GetnxnCrystalEta();
 158
 159   AliDebug(2, Form("Trigger 2x2 max amp %f, out amp %f, SM %d, iphi %d ieta %d",  maxAmp2x2, ampOutOfPatch2x2, iSM2x2,iCrystalPhi2x2, iCrystalEta2x2));
 160   AliDebug(2, Form("Trigger 4x4 max amp %f , out amp %f, SM %d, iphi %d, ieta %d",  maxAmpnxn, ampOutOfPatchnxn, iSMnxn,iCrystalPhinxn, iCrystalEtanxn));
 161
 162   Int_t iRelId2x2 []= {iSM2x2+1,0,iCrystalPhi2x2,iCrystalEta2x2};// PHOS modules in order to calculate AbsId need to be 1-5 not 0-4 as returns trigger.
 163   Int_t iAbsId2x2 =-1;
 164   Int_t iRelIdnxn []= {iSMnxn+1,0,iCrystalPhinxn,iCrystalEtanxn};// PHOS modules in order to calculate AbsId need to be 1-5 not 0-4 as returns trigger.
 165   Int_t iAbsIdnxn =-1;
 166   TVector3    pos2x2(-1,-1,-1);
 167   TVector3    posnxn(-1,-1,-1);
 168   fGeom->RelToAbsNumbering(iRelId2x2, iAbsId2x2);
 169   fGeom->RelToAbsNumbering(iRelIdnxn, iAbsIdnxn);
 170   fGeom->RelPosInAlice(iAbsId2x2, pos2x2);
 171   fGeom->RelPosInAlice(iAbsIdnxn, posnxn);
 172
 173   TArrayF triggerPosition(6);
 174   triggerPosition[0] = pos2x2(0) ;
 175   triggerPosition[1] = pos2x2(1) ;
 176   triggerPosition[2] = pos2x2(2) ;
 177   triggerPosition[3] = posnxn(0) ;
 178   triggerPosition[4] = posnxn(1) ;
 179   triggerPosition[5] = posnxn(2) ;
 180
 181   TArrayF triggerAmplitudes(4);
 182   triggerAmplitudes[0] = maxAmp2x2 ;
 183   triggerAmplitudes[1] = ampOutOfPatch2x2 ;
 184   triggerAmplitudes[2] = maxAmpnxn ;
 185   triggerAmplitudes[3] = ampOutOfPatchnxn ;
 186
 187   //esd->SetPHOSTriggerCells(triggerPosition);
 188   esd->AddPHOSTriggerPosition(triggerPosition);
 189   esd->AddPHOSTriggerAmplitudes(triggerAmplitudes);
 190
 191   //######################################
 192
 193   // Read digits array
 194   TBranch *branch = digitsTree->GetBranch("PHOS");
 195   if (!branch) {
 196     AliError("can't get the branch with the PHOS digits !");
 197     return;
 198   }
 199   TClonesArray *fDigitsArr    = new TClonesArray("AliPHOSDigit",100);
 200   branch->SetAddress(&fDigitsArr);
 201   branch->GetEntry(0);
 202
 203   // Get the clusters array
 204   TBranch *emcbranch = clustersTree->GetBranch("PHOSEmcRP");
 205   if (!emcbranch) {
 206     AliError("can't get the branch with the PHOS EMC clusters !");
 207     return;
 208   }
 209
 210   TObjArray *fEmcRecPoints = new TObjArray(100) ;
 211   emcbranch->SetAddress(&fEmcRecPoints);
 212   emcbranch->GetEntry(0);
 213
 214   //Fill CaloClusters
 215   const Float_t kBigShort = std::numeric_limits<short int>::max() - 1;
 216   const Float_t nsec100   = 1e9*100.; // units of 0.01 ns
 217   const Float_t gev500    = 500.;     // units of GeV/500
 218
 219   for (Int_t recpart = 0 ; recpart < nOfRecParticles ; recpart++) {
 220     AliPHOSRecParticle * rp = dynamic_cast<AliPHOSRecParticle*>(recParticles->At(recpart));
 221     if (Debug())
 222       rp->Print();
 223     // Get track segment and EMC rec.point associated with this rec.particle
 224     AliPHOSTrackSegment *ts    = static_cast<AliPHOSTrackSegment *>(tsm->GetTrackSegments()->At(rp->GetPHOSTSIndex()));
 225
 226     AliPHOSEmcRecPoint  *emcRP = static_cast<AliPHOSEmcRecPoint *>(fEmcRecPoints->At(ts->GetEmcIndex()));
 227     AliESDCaloCluster   *ec    = new AliESDCaloCluster() ;
 228
 229     Float_t xyz[3];
 230     for (Int_t ixyz=0; ixyz<3; ixyz++)
 231       xyz[ixyz] = rp->GetPos()[ixyz];
 232
 233     AliDebug(2,Form("Global position xyz=(%f,%f,%f)",xyz[0],xyz[1],xyz[2]));
 234
 235     //Create digits lists
 236     Int_t  digitMult  = emcRP->GetDigitsMultiplicity();
 237     Int_t *digitsList = emcRP->GetDigitsList();
 238     Float_t *rpElist   = emcRP->GetEnergiesList() ;
 239     Short_t *amplList  = new Short_t[digitMult];
 240     Short_t *timeList  = new Short_t[digitMult];
 241     Short_t *digiList  = new Short_t[digitMult];
 242
 243
 244     // Convert Float_t* and Int_t* to Short_t* to save memory
 245     for (Int_t iDigit=0; iDigit<digitMult; iDigit++) {
 246
 247       AliPHOSDigit *digit = static_cast<AliPHOSDigit *>(fDigitsArr->At(digitsList[iDigit]));
 248       amplList[iDigit] =
 249         (Short_t)(TMath::Min(rpElist[iDigit]*gev500,kBigShort)); // Energy in units of GeV/500
 250 // We should add here not full energy of digit, but unfolded one, stored in RecPoint
 251 //       amplList[iDigit] =
 252 //      (Short_t)(TMath::Min(digit->GetEnergy()*gev500,kBigShort)); // Energy in units of GeV/500
 253      timeList[iDigit] =
 254         (Short_t)(TMath::Max(-kBigShort,TMath::Min(digit->GetTime()*nsec100,kBigShort))); // time in units of 0.01 ns
 255       digiList[iDigit] = (Short_t)(digit->GetId());
 256     }
 257
 258
 259
 260     //Primaries
 261     Int_t  primMult  = 0;
 262     Int_t *primList =  emcRP->GetPrimaries(primMult);
 263
 264     // fills the ESDCaloCluster
 265
 266     ec->SetClusterType(AliESDCaloCluster::kPHOSCluster);
 267     ec->SetPosition(xyz);                 //rec.point position in MARS
 268     ec->SetE(rp->Energy());         //total particle energy
 269     ec->SetClusterDisp(emcRP->GetDispersion()); //cluster dispersion
 270     ec->SetPid(rp->GetPID()) ;                  //array of particle identification
 271     ec->SetM02(emcRP->GetM2x()) ;               //second moment M2x
 272     ec->SetM20(emcRP->GetM2z()) ;               //second moment M2z
 273     ec->SetNExMax(emcRP->GetNExMax());          //number of local maxima
 274     ec->SetEmcCpvDistance(ts->GetCpvDistance("r"));                  //Only radius, what about separate x,z????
 275     ec->SetClusterChi2(-1);                     //not yet implemented
 276     ec->SetM11(-1) ;                            //not yet implemented
 277
 278     //Digits Lists
 279     TArrayS arrayAmpList(digitMult,amplList);
 280     TArrayS arrayTimeList(digitMult,timeList);
 281     TArrayS arrayIndexList(digitMult,digiList);
 282     ec->AddDigitAmplitude(arrayAmpList);
 283     ec->AddDigitTime(arrayTimeList);
 284     ec->AddDigitIndex(arrayIndexList);
 285
 286     //Distance to the nearest bad crystal
 287     ec->SetDistanceToBadChannel(emcRP->GetDistanceToBadCrystal());
 288
 289     //Array of MC indeces
 290     TArrayI arrayPrim(primMult,primList);
 291     ec->AddLabels(arrayPrim);
 292
 293     //Array of tracks uncomment when available in future
 294     //TArrayS arrayTrackMatched(1);// Only one track, temporal solution.
 295     //arrayTrackMatched[0]= (Short_t)(matchedTrack[iClust]);
 296     //ec->AddTracksMatched(arrayTrackMatched);
 297
 298     // add the track to the esd object
 299     esd->AddCaloCluster(ec);
 300     delete ec;
 301     delete [] amplList;
 302     delete [] timeList;
 303     delete [] digiList;
 304   }
 305   fDigitsArr   ->Delete();
 306   delete fDigitsArr;
 307   fEmcRecPoints->Delete();
 308   delete fEmcRecPoints;
 309   delete tsm;
 310   delete pid;
 311 }
 312
 313 //____________________________________________________________________________
 314 AliTracker* AliPHOSReconstructor::CreateTracker() const
 315 {
 316   // creates the PHOS tracker
 317   return new AliPHOSTracker();
 318 }
 319
 320 //____________________________________________________________________________
 321 void  AliPHOSReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const
 322 {
 323   // Converts raw data to
 324   // PHOS digits
 325   // Works on a single-event basis
 326
 327   rawReader->Reset() ;
 328
 329   AliPHOSRawDecoder * dc ;
 330
 331   const TObjArray* maps = AliPHOSRecoParamEmc::GetMappings();
 332   if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
 333
 334   AliAltroMapping *mapping[4];
 335   for(Int_t i = 0; i < 4; i++) {
 336     mapping[i] = (AliAltroMapping*)maps->At(i);
 337   }
 338
 339   if(strcmp(fgkRecoParamEmc->DecoderVersion(),"v1")==0)
 340     dc=new AliPHOSRawDecoderv1(rawReader,mapping);
 341   else
 342     dc=new AliPHOSRawDecoder(rawReader,mapping);
 343
 344   TString option = GetOption();
 345   if (option.Contains("OldRCUFormat"))
 346     dc->SetOldRCUFormat(kTRUE);
 347   else
 348     dc->SetOldRCUFormat(kFALSE);
 349
 350   dc->SubtractPedestals(fgkRecoParamEmc->SubtractPedestals());
 351
 352   TClonesArray *digits = new TClonesArray("AliPHOSDigit",1);
 353   digits->SetName("DIGITS");
 354   Int_t bufsize = 32000;
 355   digitsTree->Branch("PHOS", &digits, bufsize);
 356
 357   AliPHOSRawDigiProducer pr;
 358   pr.MakeDigits(digits,dc);
 359
 360   delete dc ;
 361
 362   //ADC counts -> GeV
 363   if(strcmp(fgkRecoParamEmc->DecoderVersion(),"v1")==0){ //"Energy" calculated as fit
 364     for(Int_t i=0; i<digits->GetEntries(); i++) {
 365       AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(i);
 366       digit->SetEnergy(digit->GetEnergy()*0.005); //We assume here 5 MeV/ADC channel
 367       digit->SetTime(digit->GetTime()*1.e-7) ;    //Here we assume sample step==100 ns TO BE FIXED!!!!!!!!!!!!!
 368     }
 369   }
 370   else{ //Digits energy calculated as maximal energy
 371     for(Int_t i=0; i<digits->GetEntries(); i++) {
 372       AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(i);
 373       digit->SetEnergy(digit->GetEnergy()/AliPHOSPulseGenerator::GeV2ADC());
 374     }
 375   }
 376
 377   // Clean up digits below the noise threshold
 378   // Assuming the digit noise to be 4 MeV, we suppress digits within
 379   // 3-sigma of the noise.
 380   // This parameter should be passed via AliPHOSRecoParamEmc later
 381
 382   const Double_t emcDigitThreshold = 0.012;
 383   for(Int_t i=0; i<digits->GetEntries(); i++) {
 384     AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(i);
 385     if(digit->GetEnergy() < emcDigitThreshold)
 386       digits->RemoveAt(i) ;
 387   }
 388   digits->Compress() ;
 389
 390   //!!!!for debug!!!
 391   Int_t modMax=-111;
 392   Int_t colMax=-111;
 393   Int_t rowMax=-111;
 394   Float_t eMax=-333;
 395   //!!!for debug!!!
 396
 397   Int_t relId[4];
 398   for(Int_t iDigit=0; iDigit<digits->GetEntries(); iDigit++) {
 399     AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(iDigit);
 400     if(digit->GetEnergy()>eMax) {
 401       fGeom->AbsToRelNumbering(digit->GetId(),relId);
 402       eMax=digit->GetEnergy();
 403       modMax=relId[0];
 404       rowMax=relId[2];
 405       colMax=relId[3];
 406     }
 407   }
 408
 409   AliDebug(1,Form("Digit with max. energy:  modMax %d colMax %d rowMax %d  eMax %f\n\n",
 410                   modMax,colMax,rowMax,eMax));
 411
 412   digitsTree->Fill();
 413   digits->Delete();
 414   delete digits;
 415 }