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