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
  31 // --- Standard library ---
  32 #include <stdlib.h>
  33 #include "Riostream.h"
  34
  35 // --- AliRoot header files ---
  36
  37 #include "AliFMDdigit.h"
  38 #include "AliFMDhit.h"
  39 #include "AliFMDReconstParticles.h"
  40 #include "AliFMD.h"
  41 #include "AliFMDv1.h"
  42 #include "AliFMDReconstruction.h"
  43 #include "AliRun.h"
  44 #include "AliHeader.h"
  45 #include "AliGenEventHeader.h"
  46 ClassImp(AliFMDReconstruction)
  47
  48
  49 //____________________________________________________________________________
  50
  51 AliFMDReconstruction::AliFMDReconstruction():TTask("AliFMDReconstruction","")
  52 {
  53   fNevents = 0 ;  // Number of events to rreconnstraction, 0 means all events in current file
  54   // add Task to //root/Tasks folder
  55   TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ;
  56   roottasks->Add(this) ;
  57 }
  58 //____________________________________________________________________________
  59
  60 AliFMDReconstruction::AliFMDReconstruction(char* HeaderFile, char *ReconstParticlesFile):TTask("AliFMDReconstruction","")
  61 {
  62   fNevents = 0 ;    // Number of events to rreconnstraction, 0 means all events in current file
  63   fReconstParticlesFile=ReconstParticlesFile ;
  64   fHeadersFile=HeaderFile ;
  65   //add Task to //root/Tasks folder
  66   TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ;
  67   roottasks->Add(this) ;
  68 }
  69
  70 //____________________________________________________________________________
  71
  72 AliFMDReconstruction::~AliFMDReconstruction()
  73 {
  74 }
  75 //----------------------------------------------------------------------------
  76
  77 void AliFMDReconstruction::Exec(Option_t *option)
  78 {
  79   /*
  80     Reconstruct nember of particles
  81     in given group of pads for given FMDvolume
  82     determine by numberOfVolume ,
  83     numberOfMinSector,numberOfMaxSector,
  84     numberOfMinRing, numberOgMaxRing
  85     Reconstruction method choose dependence on number of empty pads
  86   */
  87
  88
  89   cout<<"\nStart AliFMDReconstruction::Exec(...)"<<endl;
  90   Int_t const knumVolumes=5;
  91   Int_t const knumSectors=40;
  92   Int_t const knumRings=768;
  93   Int_t padADC[10][50][800];
  94   Float_t eta, etain,etaout,rad,theta;
  95   Int_t ivol, iSector, iRing;
  96   Float_t rin[5]={4.2,15.4,4.2,15.4,4.2};
  97   Float_t rout[5]={17.4,28.4,17.4,28.4,17.4};
  98   Float_t z[5]={62.8, 75.2, -83.4, -75.2, -340.};
  99   Int_t numberOfRings[5]={768,384,768,384,768};
 100   Int_t numberOfSectors[5]=  {20,40,20,40,20};
 101   Int_t numberOfEtaIntervals[5];
 102   // number of ring for boundary 0.1 eta
 103   TBranch *brDigits=0;
 104
 105   AliFMD * fFMD = (AliFMD *) gAlice->GetDetector("FMD");
 106   TClonesArray *fReconParticles=fFMD->ReconParticles();
 107   if(fNevents == 0) fNevents=(Int_t)gAlice->TreeE()->GetEntries();
 108   cout<<" fNevents "<<fNevents<<endl;
 109   for(Int_t ievent=0;ievent<fNevents;ievent++)
 110     {
 111       cout<<" ievent "<<ievent<<endl;
 112       for (Int_t i=0; i<knumVolumes; i++)
 113         for(Int_t j=0; j<knumSectors; j++)
 114           for(Int_t ij=0; ij<knumRings; ij++)
 115             padADC[i][j][ij]=0;                    //zhachem ???
 116       gAlice->GetEvent(ievent) ;
 117       if(gAlice->TreeH()==0) return;
 118       if(gAlice->TreeD()==0) return;
 119         brDigits=gAlice->TreeD()->GetBranch("FMD");
 120         if (!brDigits){
 121           cerr<<"EXEC Branch FMD digits not found"<<endl;
 122         return;
 123       }
 124
 125       if(gAlice->TreeR()==0) gAlice->MakeTree("R");
 126       //Make branches
 127       fFMD->MakeBranch("R");
 128
 129
 130       Int_t zeroADC=1;
 131
 132       AliFMDdigit  *fmdDigit;
 133        if (fFMD)
 134         {
 135           gAlice->TreeD()->GetEvent(0);
 136           TClonesArray *fFMDdigits=fFMD->Digits();
 137           Int_t nDigits=fFMDdigits->GetEntries();
 138           cout<<" nDigits "<<nDigits<<endl;
 139            Int_t recParticles[6];
 140            Int_t nRecPart=0 ;
 141            Int_t zeroPads=0;
 142            Int_t numberOfPads=0; //To avoid warning
 143            Int_t pedestal;
 144            Float_t channelWidth=(22400*50)/1024;
 145            for (Int_t digit=0;digit<nDigits;digit++)
 146              {
 147                fmdDigit=(AliFMDdigit*)fFMDdigits->UncheckedAt(digit);
 148                ivol=fmdDigit->Volume();
 149                iSector=fmdDigit->NumberOfSector();
 150                iRing=fmdDigit->NumberOfRing();
 151                pedestal=Int_t(gRandom->Gaus(500,250));
 152                padADC[ivol-1][iSector-1][iRing-1]=
 153                  fmdDigit->ADCsignal()
 154                  -Int_t(Float_t(pedestal)/channelWidth);
 155                if (padADC[ivol-1][iSector-1][iRing-1]<0)
 156                  padADC[ivol-1][iSector-1][iRing-1]=0;
 157              } //digit loop
 158            Int_t rmin=0; Int_t rmax=0; //To avoid warning
 159            Int_t smin=0; Int_t smax=0; //To avoid warning
 160            AliHeader *header = gAlice->GetHeader();
 161            AliGenEventHeader* genHeader = header->GenEventHeader();
 162            TArrayF *o = new TArrayF(3);
 163            genHeader->PrimaryVertex(*o);
 164            Float_t zVertex=o->At(2);
 165            for (ivol=0; ivol<knumVolumes; ivol++)
 166              {
 167                Float_t realZ=zVertex+z[ivol];
 168                theta=TMath::ATan(rin[ivol]/TMath::Abs(realZ));
 169                etain = - TMath::Log( TMath::Tan(theta/2.));
 170                theta=TMath::ATan(rout[ivol]/TMath::Abs(realZ));
 171                etaout=- TMath::Log( TMath::Tan(theta/2.));
 172                numberOfEtaIntervals[ivol]=Int_t ((etain-etaout)*10);
 173                eta=etain;
 174                for (Int_t e1=0;e1<=numberOfEtaIntervals[ivol];e1++)
 175                  {
 176                    theta = 2.*TMath::ATan (TMath::Exp (-eta));
 177                    Float_t radmin = TMath::Abs(realZ) * (TMath::Tan (theta));
 178                    rmin= Int_t ( (radmin-rin[ivol])*numberOfRings[ivol]/(rout[ivol]-rin[ivol]));
 179                    eta=eta-0.1;
 180                    theta = 2.*TMath::ATan (TMath::Exp (-eta));
 181                    rad = TMath::Abs(realZ) * (TMath::Tan (theta));
 182                    rmax=Int_t( (rad-rin[ivol])*numberOfRings[ivol]/(rout[ivol]-rin[ivol]));
 183                    zeroPads=0;
 184                    smin=0;
 185                    smax=numberOfSectors[ivol];
 186                    for (Int_t iring=rmin; iring<rmax; iring++)
 187                      {
 188                        numberOfPads=(rmax-rmin)*(smax-smin);
 189                        for
 190                          (Int_t isector=1;
 191                           isector<=numberOfSectors[ivol];
 192                           isector++)
 193                          if(padADC[ivol][isector-1][iring-1]<zeroADC)
 194                            zeroPads++;
 195                      } //ring
 196                    Float_t zerosRatio=
 197                      (Float_t)zeroPads/(Float_t)numberOfPads;
 198                    recParticles[0]=ivol;
 199                    recParticles[1]=smin;
 200                    recParticles[2]=smax;
 201                    recParticles[3]=rmin;
 202                    recParticles[4]=rmax;
 203                    if (zerosRatio>0.1 )
 204                      recParticles[5]=
 205                        DeterminationByPoisson
 206                        (padADC,ivol+1,rmin,rmax,smin,smax);
 207                    else
 208                      recParticles[5]=
 209                        DeterminationByThresholds
 210                        (padADC,ivol+1,rmin,rmax,smin,smax);
 211                    new((*fReconParticles)[nRecPart++])
 212                  AliFMDReconstParticles(recParticles);
 213                  } // eta
 214              } // volume
 215
 216         }//if FMD
 217        gAlice->TreeR()->Reset();
 218        gAlice->TreeR()->Fill();
 219        gAlice->TreeR()->Write(0,TObject::kOverwrite);
 220     } //event loop
 221   cout<<"\nAliFMDReconstruction::Exec finished"<<endl;
 222 };
 223
 224 //------------------------------------------------------------
 225 Int_t AliFMDReconstruction::DeterminationByThresholds
 226 (Int_t padAdc[][50][800],Int_t volume, Int_t rmin, Int_t rmax,
 227  Int_t smin, Int_t smax)
 228 {
 229   /*
 230     reconstruct number of particles by
 231     energy deposited threshold method
 232     Using if number of empty pads less then 10%
 233
 234 */
 235   cout<<"\nStart threshold method\n";
 236
 237   Int_t thresholdInner[30]={
 238     0,     14,  28,    42,   57,
 239     72,    89,  104,  124,  129,
 240     164,  174,  179,  208,  228,
 241     248,  268,   317,  337,  357,
 242     392,  407,  416,  426,  436,
 243     461,  468,  493,  506,  515};
 244
 245   Int_t thresholdOuter[30]={0, 18, 48, 77, 105,
 246                             132, 165, 198, 231,
 247                             264, 286, 308, 334,
 248                             352, 374, 418, 440,
 249                             462, 484, 506, 528,
 250                             550, 572, 594, 616};
 251   Int_t threshold[30];
 252   for (Int_t it=0; it<30; it++) {
 253     if(volume==1||volume==3||volume==5) threshold[it]=thresholdInner[it];
 254     if(volume==2||volume==4) threshold[it]=thresholdOuter[it];
 255   }
 256   Int_t thresholdArraySize = 30;
 257   Int_t numPart=0;
 258   //Inner Si
 259   for (Int_t iring=rmin; iring<rmax; iring++)
 260     {
 261       for (Int_t isector=smin; isector<smax; isector++)
 262         {
 263           for (int i=0;i<thresholdArraySize-1;i++)
 264             {
 265               if(padAdc[volume-1][isector][iring]>threshold[i]
 266                  &&padAdc[volume-1][isector][iring]<=threshold[i+1])
 267                 {
 268                   numPart+=i;
 269                   break;
 270                 }; // if in threshol interval
 271             }; //threshold_array_size
 272         }; //iring
 273     }; //sector
 274   cout<<"\nEnd threshol method"<<endl;
 275   return numPart;
 276 }
 277
 278
 279  //__________________________________________________________________
 280
 281 Int_t AliFMDReconstruction::DeterminationByPoisson
 282 (Int_t padAdc[][50][800],
 283  Int_t vol, Int_t rmin, Int_t rmax,
 284  Int_t secmin, Int_t secmax)
 285 {
 286   /*
 287      reconstruct number of particles by Poisson statistic method
 288      according number of empty pad in chosen group of pads
 289      using if number of empty pads more then 10%
 290
 291   */
 292
 293   //  cout<<"\nStart Poisson method";
 294   Int_t thresholdADC=1;
 295   Int_t zeropads=0;
 296   for (Int_t i=rmin;i<rmax;i++)
 297     {
 298       for (Int_t j=secmin;j<secmax;j++)
 299         {
 300           if (padAdc[vol-1][j][i]<thresholdADC) zeropads++;
 301         };
 302     };
 303   Float_t lambda=-TMath::Log(Float_t(zeropads)/
 304                              ( Float_t(secmax-secmin)*
 305                                Float_t(rmax-rmin))); //+1 zdes' ne nado
 306   Int_t fRecon=(Int_t)(lambda*(secmax-secmin)*(rmax-rmin)); //+1 zdes' ne nado
 307   //  cout<<"\nEnd Poisson method"<<endl;
 308   return fRecon;
 309 };
 310
 311 //------------------------------------------------------------------
 312
 313
 314
 315
 316
 317
 318
 319
 320
 321
 322
 323
 324
 325
 326
 327
 328
 329
 330
 331
 332
 333
 334
 335
 336
 337
 338
 339
 340