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