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* gime = fRunLoader->GetLoader("FMDLoader");
  78   if (gime == 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* gime = fRunLoader->GetLoader("FMDLoader");
  92   gime->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   cout<<"\nStart AliFMDReconstruction::Exec(...)"<<endl;
 111   Int_t const knumVolumes=5;
 112   Int_t padADC;
 113   Float_t eta, etain,etaout,rad,theta;
 114   Int_t ivol, iSector, iRing;
 115   Float_t rin[5]={4.2,15.4,4.2,15.4,4.2};
 116   Float_t rout[5]={17.4,28.4,17.4,28.4,17.4};
 117   Float_t z[5]={-62.8, -75.2, 83.4, 75.2, 340.};
 118   Int_t numberOfRings[5]={512,256,512,256,512};
 119   Int_t numberOfSectors[5]=  {20,40,20,40,20};
 120   Int_t numberOfEtaIntervals[5];
 121   // number of ring for boundary 0.1 eta
 122
 123
 124   if (fRunLoader == 0x0)
 125    {
 126     Error("Exec","Run Loader loader is NULL - Session not opened");
 127     return;
 128    }
 129
 130
 131  AliLoader* gime = fRunLoader->GetLoader("FMDLoader");
 132   if (gime == 0x0)
 133    {
 134      Fatal("AliFMDReconstruction","Can not find FMD (loader) in specified event");
 135      return;//never reached
 136    }
 137
 138   fRunLoader->LoadgAlice();
 139   fRunLoader->LoadHeader();
 140   Int_t retval;
 141   TDirectory* cwd = gDirectory;
 142   gDirectory = 0x0;
 143     Text_t buf1[20];
 144   TH2F* hTotal[10];
 145   for (Int_t j=1; j<=5; j++){
 146     sprintf(buf1,"hTotal%d",j);
 147
 148     hTotal[j] = new TH2F(buf1," Number of primary particles ",
 149                          numberOfSectors[j-1],1,numberOfSectors[j-1],
 150                          numberOfRings[j-1],1,numberOfRings[j-1]);
 151   }
 152   gDirectory = cwd;
 153
 154
 155   if(fNevents == 0) fNevents=Int_t (fRunLoader->TreeE()->GetEntries());
 156   cout<<" fNevents "<<fNevents<<endl;
 157    for(Int_t ievent=0;ievent<fNevents;ievent++)
 158     {
 159       fRunLoader->GetEvent(ievent) ;
 160
 161       cout<<" ievent "<<ievent<<endl;
 162     for (Int_t i=0; i<5; i++)
 163         hTotal[i+1]->Reset();
 164
 165       retval = gime->LoadDigits("READ");
 166       if (retval)
 167         {
 168           Error("Exec","Error occured while loading digits. Exiting.");
 169           return;
 170         }
 171
 172       AliFMD * fFMD = (AliFMD *) gAlice->GetDetector("FMD");
 173       TClonesArray *fReconParticles=fFMD->ReconParticles();
 174       TClonesArray *fDigits=fFMD->Digits();
 175
 176       TTree* treeD = gime->TreeD();
 177       if (treeD == 0x0)
 178         {
 179           Error("Exec","Can not get Tree with Digits. Nothing to reconstruct - Exiting");
 180           return;
 181         }
 182
 183       TBranch *brDigits=0;
 184
 185       brDigits=treeD->GetBranch("FMD");
 186       cout<<" brDigits "<<brDigits<<endl;
 187       if (brDigits) {
 188         brDigits->SetAddress(&fDigits);
 189         //        fFMD->SetHitsAddressBranch(brHits);
 190       }else{
 191         cerr<<"EXEC Branch FMD digits not found"<<endl;
 192         return;
 193       }
 194
 195       if(gime->TreeR()==0) gime->MakeTree("R");
 196
 197       //Make branches
 198       fFMD->MakeBranch("R");
 199
 200
 201       Int_t zeroADC=6;
 202
 203       AliFMDdigit  *fmdDigit;
 204        if (fFMD)
 205         {
 206           gime->TreeD()->GetEvent(0);
 207
 208           Int_t nDigits=fDigits->GetEntries();
 209           cout<<" nDigits "<<nDigits<<endl;
 210           Int_t recParticles[6];
 211           Int_t nRecPart=0 ;
 212            Int_t zeroPads=0;
 213            Int_t numberOfPads=0; //To avoid warning
 214            Int_t pedestal;
 215            Float_t channelWidth=(22400*50)/1024;
 216            for (Int_t digit=0;digit<nDigits;digit++)
 217              {
 218                fmdDigit=(AliFMDdigit*)fDigits->UncheckedAt(digit);
 219                ivol=fmdDigit->Volume();
 220                iSector=fmdDigit->NumberOfSector();
 221                iRing=fmdDigit->NumberOfRing();
 222                pedestal=Int_t(gRandom->Gaus(500,250));
 223                padADC= fmdDigit->ADCsignal()
 224                  -Int_t(Float_t(pedestal)/channelWidth);
 225                if (padADC<0) padADC=0;
 226               hTotal[ivol]->Fill(iSector,iRing,padADC);
 227              } //digit loop
 228            Int_t rmin=0; Int_t rmax=0; //To avoid warning
 229            Int_t smin=0; Int_t smax=0; //To avoid warning
 230            AliHeader *header = fRunLoader->GetHeader();
 231            AliGenEventHeader* genHeader = header->GenEventHeader();
 232            TArrayF *o = new TArrayF(3);
 233            genHeader->PrimaryVertex(*o);
 234            Float_t zVertex=o->At(2);
 235            for (ivol=0; ivol<knumVolumes; ivol++)
 236              {
 237                Float_t ring2number=Float_t (numberOfRings[ivol])/
 238                  (rout[ivol]-rin[ivol]);
 239                Float_t realZ=zVertex+z[ivol];
 240                theta=TMath::ATan(rout[ivol]/TMath::Abs(realZ));
 241                etain = - TMath::Log( TMath::Tan(theta/2.));
 242                theta=TMath::ATan(rin[ivol]/TMath::Abs(realZ));
 243                etaout=- TMath::Log( TMath::Tan(theta/2.));
 244                numberOfEtaIntervals[ivol]=Int_t((etaout-etain)*10)-1;
 245                eta=etain;
 246                for (Int_t e1=0;e1<=numberOfEtaIntervals[ivol];e1++)
 247                  {
 248                    theta = 2.*TMath::ATan (TMath::Exp (-eta));
 249                    Float_t radmin = TMath::Abs(realZ) * (TMath::Tan (theta));
 250                    rmax= Int_t ( (radmin-rin[ivol])*ring2number+0.5);
 251                    cout<<ivol<<" "<<" eta "<<eta<<" radmin "<<radmin<<
 252
 253                      " Rmax "<<rmax<<" "<<endl;;
 254                    eta=eta+0.1;
 255                    theta = 2.*TMath::ATan (TMath::Exp (-eta));
 256                    rad = TMath::Abs(realZ) * (TMath::Tan (theta));
 257                    rmin=Int_t( (rad-rin[ivol])*ring2number+0.5);
 258                    cout<<"eta "<<eta<<" rad "<<rad<<" Rmin "<<rmin<<endl;
 259                    //      if(ivol==2&&e1==13) rmin=0;
 260                    zeroPads=0;
 261                    smin=0;
 262                    smax=numberOfSectors[ivol];
 263                    numberOfPads=(rmax-rmin)*(smax-smin);
 264                    for (Int_t iring=rmin; iring<rmax; iring++)
 265                      {
 266                        for
 267                          (Int_t isector=0;
 268                           isector<numberOfSectors[ivol];
 269                           isector++)
 270                          {
 271                            if(hTotal[ivol+1]->GetBinContent(isector+1,iring+1)
 272                             <zeroADC) zeroPads++;}
 273
 274                      } //ring
 275                    Float_t numberOfPads=Float_t(smax-smin)*Float_t(rmax-rmin);
 276
 277                    cout<<" zero "<<zeroPads++<<" pads "<<numberOfPads;
 278                    Double_t lambda=-TMath::Log(Double_t(zeroPads)/numberOfPads);
 279                    Int_t fRecon=Int_t (lambda*numberOfPads+0.5);
 280
 281                    Float_t zerosRatio=
 282                      (Float_t)zeroPads/(Float_t)numberOfPads;
 283                    cout<<" zerosRatio "<<zerosRatio<<" recon "<<fRecon<<endl;
 284                    recParticles[0]=ivol+1;
 285                    recParticles[1]=smin;
 286                    recParticles[2]=smax;
 287                    recParticles[3]=rmin;
 288                    recParticles[4]=rmax;
 289                    recParticles[5]= fRecon;
 290                    new((*fReconParticles)[nRecPart++])
 291                      AliFMDReconstParticles(recParticles);
 292
 293
 294                  } // eta
 295              } // volume
 296
 297         }//if FMD
 298        gime->TreeR()->Reset();
 299        gime->TreeR()->Fill();
 300        gime->WriteRecPoints("OVERWRITE");
 301     } //event loop
 302   cout<<"\nAliFMDReconstruction::Exec finished"<<endl;
 303 }
 304