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  //Collects all digits in the same active volume into number of particles
  80   cout<<"\nStart AliFMDReconstruction::Exec(...)"<<endl;
  81   Int_t const NumVolums=5;
  82   Int_t const NumSectors=40;
  83   Int_t const NumRings=768;
  84   Int_t PadADC[10][50][800];
  85   Float_t eta, etain,etaout,rad,theta;
  86   Int_t ivol, iSector, iRing;
  87   Float_t rin[5]={4.2,15.4,4.2,15.4,4.2};
  88   Float_t rout[5]={17.4,28.4,17.4,28.4,17.4};
  89   Float_t z[5]={62.8, 75.2, -83.4, -75.2, -340.};
  90   Int_t NumberOfRings[5]={768,384,768,384,768};
  91   Int_t NumberOfSectors[5]=  {20,40,20,40,20};
  92   Int_t NumberOfEtaIntervals[5];
  93   // number of ring for boundary 0.1 eta
  94   TBranch *brDigits=0;
  95
  96   AliFMD * FMD = (AliFMD *) gAlice->GetDetector("FMD");
  97   TClonesArray *fReconParticles=FMD->ReconParticles();
  98   if(fNevents == 0) fNevents=(Int_t)gAlice->TreeE()->GetEntries();
  99   cout<<" fNevents "<<fNevents<<endl;
 100   for(Int_t ievent=0;ievent<fNevents;ievent++)
 101     {
 102       cout<<" ievent "<<ievent<<endl;
 103       for (Int_t i=0; i<NumVolums; i++)
 104         for(Int_t j=0; j<NumSectors; j++)
 105           for(Int_t ij=0; ij<NumRings; ij++)
 106             PadADC[i][j][ij]=0;                    //zhachem ???
 107       gAlice->GetEvent(ievent) ;
 108       if(gAlice->TreeH()==0) return;
 109       if(gAlice->TreeD()==0) return;
 110         brDigits=gAlice->TreeD()->GetBranch("FMD");
 111         if (!brDigits){
 112           cerr<<"EXEC Branch FMD digits not found"<<endl;
 113         return;
 114       }
 115
 116       if(gAlice->TreeR()==0) gAlice->MakeTree("R");
 117       //Make branches
 118       FMD->MakeBranch("R");
 119
 120
 121       Int_t zeroADC=1;
 122
 123       AliFMDdigit  *fmdDigit;
 124        if (FMD)
 125         {
 126           gAlice->TreeD()->GetEvent(0);
 127           TClonesArray *FMDdigits=FMD->Digits();
 128           Int_t nDigits=FMDdigits->GetEntries();
 129           cout<<" nDigits "<<nDigits<<endl;
 130            Int_t RecParticles[6];
 131            Int_t nRecPart=0 ;
 132            Int_t ZeroPads=0;
 133            Int_t NumberOfPads=0; //To avoid warning
 134            Int_t pedestal;
 135            Float_t channelWidth=(22400*50)/1024;
 136            for (Int_t digit=0;digit<nDigits;digit++)
 137              {
 138                fmdDigit=(AliFMDdigit*)FMDdigits->UncheckedAt(digit);
 139                ivol=fmdDigit->Volume();
 140                iSector=fmdDigit->NumberOfSector();
 141                iRing=fmdDigit->NumberOfRing();
 142                pedestal=Int_t(gRandom->Gaus(500,250));
 143                PadADC[ivol-1][iSector-1][iRing-1]=
 144                  fmdDigit->ADCsignal()
 145                  -Int_t(Float_t(pedestal)/channelWidth);
 146                if (PadADC[ivol-1][iSector-1][iRing-1]<0)
 147                  PadADC[ivol-1][iSector-1][iRing-1]=0;
 148              } //digit loop
 149            Int_t Rmin=0; Int_t Rmax=0; //To avoid warning
 150            Int_t Smin=0; Int_t Smax=0; //To avoid warning
 151            AliHeader *header = gAlice->GetHeader();
 152            AliGenEventHeader* genHeader = header->GenEventHeader();
 153            TArrayF *o = new TArrayF(3);
 154            genHeader->PrimaryVertex(*o);
 155            Float_t zVertex=o->At(2);
 156            for (ivol=0; ivol<NumVolums; ivol++)
 157              {
 158                Float_t realZ=zVertex+z[ivol];
 159                theta=TMath::ATan(rin[ivol]/TMath::Abs(realZ));
 160                etain = - TMath::Log( TMath::Tan(theta/2.));
 161                theta=TMath::ATan(rout[ivol]/TMath::Abs(realZ));
 162                etaout=- TMath::Log( TMath::Tan(theta/2.));
 163                NumberOfEtaIntervals[ivol]=Int_t ((etain-etaout)*10);
 164                eta=etain;
 165                for (Int_t e1=0;e1<=NumberOfEtaIntervals[ivol];e1++)
 166                  {
 167                    theta = 2.*TMath::ATan (TMath::Exp (-eta));
 168                    Float_t radmin = TMath::Abs(realZ) * (TMath::Tan (theta));
 169                    Rmin= Int_t ( (radmin-rin[ivol])*NumberOfRings[ivol]/(rout[ivol]-rin[ivol]));
 170                    eta=eta-0.1;
 171                    theta = 2.*TMath::ATan (TMath::Exp (-eta));
 172                    rad = TMath::Abs(realZ) * (TMath::Tan (theta));
 173                    Rmax=Int_t( (rad-rin[ivol])*NumberOfRings[ivol]/(rout[ivol]-rin[ivol]));
 174                    ZeroPads=0;
 175                    Smin=0;
 176                    Smax=NumberOfSectors[ivol];
 177                    for (Int_t iring=Rmin; iring<Rmax; iring++)
 178                      {
 179                        NumberOfPads=(Rmax-Rmin)*(Smax-Smin);
 180                        for
 181                          (Int_t isector=1;
 182                           isector<=NumberOfSectors[ivol];
 183                           isector++)
 184                          if(PadADC[ivol][isector-1][iring-1]<zeroADC)
 185                            ZeroPads++;
 186                      } //ring
 187                    Float_t zerosRatio=
 188                      (Float_t)ZeroPads/(Float_t)NumberOfPads;
 189                    RecParticles[0]=ivol;
 190                    RecParticles[1]=Smin;
 191                    RecParticles[2]=Smax;
 192                    RecParticles[3]=Rmin;
 193                    RecParticles[4]=Rmax;
 194                    if (zerosRatio>0.1 ||(ivol==2||ivol==4))
 195                      RecParticles[5]=
 196                        Determination_by_Poisson
 197                        (PadADC,ivol+1,Rmin,Rmax,Smin,Smax);
 198                    else
 199                      RecParticles[5]=
 200                        Determination_by_thresholds
 201                        (PadADC,ivol+1,Rmin,Rmax,Smin,Smax);
 202                    new((*fReconParticles)[nRecPart++])
 203                  AliFMDReconstParticles(RecParticles);
 204                  } // eta
 205              } // volume
 206
 207         }//if FMD
 208        gAlice->TreeR()->Reset();
 209        gAlice->TreeR()->Fill();
 210        gAlice->TreeR()->Write(0,TObject::kOverwrite);
 211     } //event loop
 212   cout<<"\nAliFMDReconstruction::Exec finished"<<endl;
 213 };
 214
 215 //------------------------------------------------------------
 216 Int_t AliFMDReconstruction::Determination_by_thresholds
 217 (Int_t PadADC[][50][800],Int_t volume, Int_t Rmin, Int_t Rmax,
 218  Int_t Smin, Int_t Smax)
 219 {
 220   cout<<"\nStart threshold method\n";
 221
 222   Int_t thresholdInner[30]={
 223     0,     14,  28,    42,   57,
 224     72,    89,  104,  124,  129,
 225     164,  174,  179,  208,  228,
 226     248,  268,   317,  337,  357,
 227     392,  407,  416,  426,  436,
 228     461,  468,  493,  506,  515};
 229
 230   Int_t thresholdOuter[30]={0, 18, 48, 77, 105,
 231                             132, 165, 198, 231,
 232                             264, 286, 308, 334,
 233                             352, 374, 418, 440,
 234                             462, 484, 506, 528,
 235                             550, 572, 594, 616};
 236   Int_t threshold[30];
 237   for (Int_t it=0; it<30; it++) {
 238     if(volume==1||volume==3||volume==5) threshold[it]=thresholdInner[it];
 239     if(volume==2||volume==4) threshold[it]=thresholdOuter[it];
 240   }
 241   Int_t threshold_array_size = 30;
 242   Int_t NumPart=0;
 243   //Inner Si
 244   for (Int_t iring=Rmin; iring<Rmax; iring++)
 245     {
 246       for (Int_t isector=Smin; isector<Smax; isector++)
 247         {
 248           for (int i=0;i<threshold_array_size-1;i++)
 249             {
 250               if(PadADC[volume-1][isector][iring]>threshold[i]
 251                  &&PadADC[volume-1][isector][iring]<=threshold[i+1])
 252                 {
 253                   NumPart+=i;
 254                   break;
 255                 }; // if in threshol interval
 256             }; //threshold_array_size
 257         }; //iring
 258     }; //sector
 259   cout<<"\nEnd threshol method"<<endl;
 260   return NumPart;
 261 }
 262
 263
 264  //__________________________________________________________________
 265
 266 Int_t AliFMDReconstruction::Determination_by_Poisson (Int_t PadADC[][50][800],
 267                                            Int_t vol, Int_t rmin, Int_t rmax,
 268                                            Int_t secmin, Int_t secmax)
 269 {
 270   //  cout<<"\nStart Poisson method";
 271   Int_t threshold_adc=1;
 272   Int_t zeropads=0;
 273   for (Int_t i=rmin;i<rmax;i++)
 274     {
 275       for (Int_t j=secmin;j<secmax;j++)
 276         {
 277           if (PadADC[vol-1][j][i]<threshold_adc) zeropads++;
 278         };
 279     };
 280   Float_t lambda=-TMath::Log(Float_t(zeropads)/
 281                              ( Float_t(secmax-secmin)*
 282                                Float_t(rmax-rmin))); //+1 zdes' ne nado
 283   Int_t Recon=(Int_t)(lambda*(secmax-secmin)*(rmax-rmin)); //+1 zdes' ne nado
 284   //  cout<<"\nEnd Poisson method"<<endl;
 285   return Recon;
 286 };
 287
 288 //------------------------------------------------------------------
 289
 290
 291
 292
 293
 294
 295
 296
 297
 298
 299
 300
 301
 302
 303
 304
 305
 306
 307
 308
 309
 310
 311
 312
 313
 314
 315
 316
 317