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