FMD/AliFMDReconstructor.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 //_________________________________________________________________________
  17 // This is a class that constructs ReconstParticles (reconstructed particles)
  18 // out of Digits
  19 //
  20 //-- Authors: Evgeny Karpechev(INR) and Alla Maevsksia
  21 //////////////////////////////////////////////////////////////////////////////
  22
  23 // --- ROOT system ---
  24 #include "TTask.h"
  25 #include "TTree.h"
  26 #include "TSystem.h"
  27 #include "TFile.h"
  28 #include "TROOT.h"
  29 #include "TFolder.h"
  30 #include "TH2F.h"
  31 #include "TRandom.h"
  32
  33 // --- Standard library ---
  34 #include <stdlib.h>
  35 #include <Riostream.h>
  36
  37 // --- AliRoot header files ---
  38
  39 #include "AliRunLoader.h"
  40 #include "AliLoader.h"
  41
  42 #include "AliFMDdigit.h"
  43 #include "AliFMDReconstParticles.h"
  44 #include "AliFMDReconstructor.h"
  45 #include "AliRun.h"
  46 #include "AliConfig.h"
  47 #include "AliHeader.h"
  48 #include "AliGenEventHeader.h"
  49
  50 ClassImp(AliFMDReconstructor)
  51
  52
  53 //____________________________________________________________________________
  54
  55 void AliFMDReconstructor::Reconstruct(AliRunLoader* runLoader) const
  56 {
  57  //Collects all digits in the same active volume into number of particles
  58   /*
  59     Reconstruct number of particles
  60     in given group of pads for given FMDvolume
  61     determine by numberOfVolume ,
  62     numberOfMinSector,numberOfMaxSector,
  63     numberOfMinRing, numberOgMaxRing
  64     Reconstructor method choose dependence on number of empty pads
  65   */
  66
  67
  68 #ifdef DEBUG
  69   Info("Exec ","Start");
  70 #endif
  71
  72   Int_t const knumVolumes=5;
  73   Int_t padADC;
  74   Float_t eta, etain,etaout,rad,theta;
  75   Int_t ivol, iSector, iRing;
  76   Float_t rin[5]={4.2,15.4,4.2,15.4,4.2};
  77   Float_t rout[5]={17.4,28.4,17.4,28.4,17.4};
  78   Float_t z[5]={-62.8, -75.2, 83.4, 75.2, 340.};
  79   Int_t numberOfRings[5]={512,256,512,256,512};
  80   Int_t numberOfSectors[5]=  {20,40,20,40,20};
  81   Int_t numberOfEtaIntervals[5];
  82   // number of ring for boundary 0.1 eta
  83
  84
  85   if (runLoader == 0x0)
  86    {
  87     Error("Exec","Run Loader loader is NULL - Session not opened");
  88     return;
  89    }
  90
  91   AliLoader* plFMD = runLoader->GetLoader("FMDLoader");
  92   if (plFMD == 0x0)
  93    {
  94      Fatal("AliFMDReconstructor","Can not find FMD (loader) in specified event");
  95      return;//never reached
  96    }
  97
  98   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
  99   if (!runLoader->TreeE()) runLoader->LoadHeader();
 100
 101   TDirectory* cwd = gDirectory;
 102   gDirectory = 0x0;
 103   Text_t buf1[20];
 104   TH2F* hTotal[10];
 105   for (Int_t j=1; j<=5; j++){
 106     sprintf(buf1,"hTotal%d",j);
 107
 108     hTotal[j] = new TH2F(buf1," Number of primary particles ",
 109                          numberOfSectors[j-1],1,numberOfSectors[j-1],
 110                          numberOfRings[j-1],1,numberOfRings[j-1]);
 111   }
 112   gDirectory = cwd;
 113
 114   plFMD->LoadRecPoints("RECREATE");
 115   TClonesArray* reconParticles = new TClonesArray("AliFMDReconstParticles");
 116   TClonesArray* digits = new TClonesArray("AliFMDdigit");
 117
 118   Int_t retval=0;
 119   Int_t nevents=Int_t (runLoader->TreeE()->GetEntries());
 120 #ifdef DEBUG
 121   cout<<" nevents "<<nevents<<endl;
 122 #endif
 123    for(Int_t ievent=0;ievent<nevents;ievent++)
 124     {
 125 #ifdef DEBUG
 126       cout<<" *** event "<<ievent<<endl;
 127 #endif
 128       runLoader->GetEvent(ievent) ;
 129       //event z-vertex for correction eta-rad dependence
 130       AliHeader *header = runLoader->GetHeader();
 131       AliGenEventHeader* genHeader = header->GenEventHeader();
 132       TArrayF *o = new TArrayF(3);
 133       if (genHeader) genHeader->PrimaryVertex(*o);
 134       Float_t zVertex=o->At(2);
 135
 136       for (Int_t i=0; i<5; i++)
 137         hTotal[i+1]->Reset();
 138
 139       retval = plFMD->LoadDigits("READ");
 140       if (retval)
 141         {
 142           Error("Exec","Error occured while loading digits. Exiting.");
 143           return;
 144         }
 145
 146       TTree* treeD = plFMD->TreeD();
 147       if (treeD == 0x0)
 148         {
 149           Error("Exec","Can not get Tree with Digits. Nothing to reconstruct - Exiting");
 150           return;
 151         }
 152
 153       TBranch *brDigits=0;
 154
 155       brDigits=treeD->GetBranch("FMD");
 156
 157       if (brDigits) {
 158         brDigits->SetAddress(&digits);
 159       }else{
 160         cerr<<"EXEC Branch FMD digits not found"<<endl;
 161         return;
 162       }
 163
 164       if(plFMD->TreeR()==0) plFMD->MakeTree("R");
 165
 166       //Make branches
 167       const Int_t kBufferSize = 16000;
 168       plFMD->TreeR()->Branch("FMD", &reconParticles, kBufferSize);
 169
 170
 171       Int_t zeroADC=6;
 172       // read Digits
 173       AliFMDdigit  *fmdDigit;
 174        if (plFMD->TreeD()->GetEvent(0))
 175         {
 176           Int_t nDigits=digits->GetEntries();
 177           Int_t recParticles[6];
 178           Int_t nRecPart=0 ;
 179           Int_t zeroPads=0;
 180           Int_t numberOfPads=0;
 181           Int_t pedestal;
 182           Float_t channelWidth=(22400*50)/1024;
 183           for (Int_t digit=0;digit<nDigits;digit++)
 184             {
 185               fmdDigit=(AliFMDdigit*)digits->UncheckedAt(digit);
 186               ivol=fmdDigit->Volume();
 187               iSector=fmdDigit->NumberOfSector();
 188               iRing=fmdDigit->NumberOfRing();
 189               pedestal=Int_t(gRandom->Gaus(500,250));
 190               padADC= fmdDigit->ADCsignal()
 191                 -Int_t(Float_t(pedestal)/channelWidth);
 192               if (padADC<0) padADC=0;
 193               hTotal[ivol]->Fill(iSector,iRing,padADC);
 194             } //digit loop
 195
 196           //reconstruct multiplicity in 0.1 eta according Poisson distribution
 197
 198           Int_t rmin=0; Int_t rmax=0;
 199           Int_t smin=0; Int_t smax=0;
 200           for (ivol=0; ivol<knumVolumes; ivol++)
 201             {
 202               Float_t ring2number=Float_t (numberOfRings[ivol])/
 203                 (rout[ivol]-rin[ivol]);
 204               Float_t realZ=zVertex+z[ivol];
 205               theta=TMath::ATan(rout[ivol]/TMath::Abs(realZ));
 206               etain = - TMath::Log( TMath::Tan(theta/2.));
 207               theta=TMath::ATan(rin[ivol]/TMath::Abs(realZ));
 208               etaout=- TMath::Log( TMath::Tan(theta/2.));
 209               numberOfEtaIntervals[ivol]=Int_t((etaout-etain)*10)-1;
 210               eta=etain;
 211               for (Int_t e1=0;e1<=numberOfEtaIntervals[ivol];e1++)
 212                 {
 213                   theta = 2.*TMath::ATan (TMath::Exp (-eta));
 214                   Float_t radmin = TMath::Abs(realZ) * (TMath::Tan (theta));
 215                   rmax= Int_t ( (radmin-rin[ivol])*ring2number+0.5);
 216                   eta=eta+0.1;
 217                   theta = 2.*TMath::ATan (TMath::Exp (-eta));
 218                   rad = TMath::Abs(realZ) * (TMath::Tan (theta));
 219                   rmin=Int_t( (rad-rin[ivol])*ring2number+0.5);
 220
 221                   zeroPads=0;
 222                   smin=0;
 223                   smax=numberOfSectors[ivol];
 224                   numberOfPads=(rmax-rmin)*(smax-smin);
 225                   for (Int_t iring=rmin; iring<rmax; iring++)
 226                     {
 227                       for
 228                         (Int_t isector=0;
 229                          isector<numberOfSectors[ivol];
 230                          isector++)
 231                         {
 232                           if(hTotal[ivol+1]->GetBinContent(isector+1,iring+1)
 233                              <zeroADC) zeroPads++;}
 234
 235                     } //ring
 236                   Float_t numberOfPads=Float_t(smax-smin)*Float_t(rmax-rmin);
 237
 238                   Double_t lambda=-TMath::Log(Double_t(zeroPads)/numberOfPads);
 239                   Int_t fRecon=Int_t (lambda*numberOfPads+0.5);
 240                   recParticles[0]=ivol+1;
 241                   recParticles[1]=smin;
 242                   recParticles[2]=smax;
 243                   recParticles[3]=rmin;
 244                   recParticles[4]=rmax;
 245                   recParticles[5]= fRecon;
 246                   new((*reconParticles)[nRecPart++])   AliFMDReconstParticles(recParticles);
 247
 248
 249                  } // eta
 250              } // volume
 251
 252         }//if (plFMD->TreeD()->GetEvent(0))
 253        plFMD->TreeR()->Reset();
 254        plFMD->TreeR()->Fill();
 255        plFMD->WriteRecPoints("OVERWRITE");
 256        plFMD->UnloadDigits();
 257     } //event loop
 258   plFMD->UnloadRecPoints();
 259 #ifdef DEBUG
 260   Info(" Exec"," finished");
 261 #endif
 262   //  delete hTotal[10];
 263
 264 }
 265
 266
 267 //_____________________________________________________________________________
 268 void AliFMDReconstructor::FillESD(AliRunLoader* /*runLoader*/,
 269                                   AliESD* /*esd*/) const
 270 {
 271 // nothing to be done
 272
 273 }
 274