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 <Riostream.h>
  31
  32 // --- Standard library ---
  33 #include <stdlib.h>
  34
  35 // --- AliRoot header files ---
  36
  37 #include "AliFMDdigit.h"
  38 #include "AliFMDReconstParticles.h"
  39 #include "AliFMD.h"
  40 #include "AliFMDv1.h"
  41 #include "AliFMDReconstruction.h"
  42 #include "AliRun.h"
  43 ClassImp(AliFMDReconstruction)
  44
  45
  46 //____________________________________________________________________________
  47
  48 AliFMDReconstruction::AliFMDReconstruction():TTask("AliFMDReconstruction","")
  49 {
  50   fNevents = 0 ;  // Number of events to rreconnstraction, 0 means all events in current file
  51   // add Task to //root/Tasks folder
  52   TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ;
  53   roottasks->Add(this) ;
  54 }
  55 //____________________________________________________________________________
  56
  57 AliFMDReconstruction::AliFMDReconstruction(char* HeaderFile, char *ReconstParticlesFile):TTask("AliFMDReconstruction","")
  58 {
  59   fNevents = 0 ;    // Number of events to rreconnstraction, 0 means all events in current file
  60   fReconstParticlesFile=ReconstParticlesFile ;
  61   fHeadersFile=HeaderFile ;
  62   //add Task to //root/Tasks folder
  63   TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ;
  64   roottasks->Add(this) ;
  65 }
  66
  67 //____________________________________________________________________________
  68
  69 AliFMDReconstruction::~AliFMDReconstruction()
  70 {
  71 }
  72 //----------------------------------------------------------------------------
  73
  74 void AliFMDReconstruction::Exec(Option_t *option)
  75 {
  76   printf (" AliFMDReconstruction starting \n");
  77  //Collects all digits in the same active volume into number of particles
  78
  79   Int_t const NumVolums=5;
  80   Int_t const NumSectors=40;
  81   Int_t const NumRings=768;
  82   Int_t PadADC[10][50][800];
  83   Int_t ivol, iSector, iRing;
  84   Int_t Ne1;
  85   //Int_t NumberOfRings[5]=  {256,128,256,128,256};
  86   Int_t NumberOfSectors[5]=  {20,40,20,40,20};
  87   // number of ring for boundary 0.1 eta
  88    Int_t EtaIntervalInner []=
  89       {0, 55, 110, 165, 221, 276, 331, 386, 442,
  90         497, 552, 607, 663, 718, 767 };
  91    /*
  92      {0,  18, 36,   55,  73,
  93       92, 110, 128, 147, 165,
  94       184, 202, 221, 239, 255};*/
  95    Int_t  EtaIntervalOuter []=  //{0, 21, 43, 65, 86, 108, 127};
  96      {0, 65, 130, 195, 260, 325, 383};
  97  Int_t EtaInterval[20];
  98
  99
 100   Int_t size_EtaIntervalInner=sizeof (EtaIntervalInner)/sizeof(EtaIntervalInner[0]);
 101
 102   Int_t size_EtaIntervalOuter=sizeof (EtaIntervalOuter)/sizeof(EtaIntervalOuter[0]);
 103
 104   AliFMD * FMD = (AliFMD *) gAlice->GetDetector("FMD");
 105   TClonesArray *fReconParticles=FMD->ReconParticles();
 106   if(fNevents == 0) fNevents=(Int_t)gAlice->TreeE()->GetEntries();
 107   for(Int_t ievent=0;ievent<fNevents;ievent++)
 108     {
 109       for (Int_t i=0; i<NumVolums; i++)
 110         for(Int_t j=0; j<NumSectors; j++)
 111           for(Int_t ij=0; ij<NumRings; ij++)
 112             PadADC[i][j][ij]=0;                    //zhachem ???
 113       gAlice->GetEvent(ievent) ;
 114       if(gAlice->TreeH()==0) return;
 115       if(gAlice->TreeR()==0) gAlice->MakeTree("R");
 116       //Make branches
 117       AliFMDdigit  *fmdDigit;
 118       FMD->MakeBranch("R");
 119
 120       Int_t zeroADC=1;
 121       // Int_t  threshold_array_size=30;
 122
 123       // cout<<" AliFMDdigit "<<AliFMDdigit<<endl;
 124       if (FMD)
 125         {
 126           gAlice->TreeD()->GetEvent(0);
 127           TClonesArray *FMDdigits=FMD->Digits();
 128           Int_t nDigits=FMDdigits->GetEntries();
 129            Int_t RecParticles[6];
 130            Int_t nRecPart=0 ;
 131            Int_t ZeroPads=0;
 132            Int_t NumberOfPads=0; //To avoid warning
 133            Int_t pedestal;
 134            Float_t channelWidth=(22400*50)/1024;
 135            for (Int_t digit=0;digit<nDigits;digit++)
 136              {
 137                fmdDigit=(AliFMDdigit*)FMDdigits->UncheckedAt(digit);
 138                ivol=fmdDigit->Volume();
 139                iSector=fmdDigit->NumberOfSector();
 140                iRing=fmdDigit->NumberOfRing();
 141                pedestal=Int_t(gRandom->Gaus(500,250));
 142                PadADC[ivol-1][iSector-1][iRing-1]=
 143                  fmdDigit->ADCsignal()
 144                  -Int_t(Float_t(pedestal)/channelWidth);
 145                if (PadADC[ivol-1][iSector-1][iRing-1]<0)
 146                  PadADC[ivol-1][iSector-1][iRing-1]=0;
 147              } //digit loop
 148            Int_t Rmin=0; Int_t Rmax=0; //To avoid warning
 149            Int_t Smin=0; Int_t Smax=0; //To avoid warning
 150
 151            for (ivol=1; ivol<=NumVolums; ivol++)
 152              {
 153                if (ivol==1||ivol==3||ivol==5)
 154                  {
 155                    Ne1=size_EtaIntervalInner;
 156                    for(Int_t ieta=0; ieta<Ne1; ieta++)
 157                      EtaInterval[ieta]=EtaIntervalInner[ieta];
 158                  }
 159                if (ivol==2||ivol==4)
 160                  {
 161                    Ne1=size_EtaIntervalOuter;
 162                    for( Int_t ieta=0; ieta<Ne1; ieta++)
 163                      EtaInterval[ieta]=EtaIntervalOuter[ieta];
 164                  }
 165
 166
 167                    for (Int_t e1=0;e1<Ne1-1;e1++) // vol<=NumVolums
 168                      {
 169                        Rmin=EtaInterval[e1];
 170                    Rmax=EtaInterval[e1+1];
 171                    ZeroPads=0;
 172                    Smin=0;
 173                    Smax=NumberOfSectors[ivol-1];
 174                    for (Int_t iring=Rmin; iring<Rmax; iring++)
 175                      {
 176                        NumberOfPads=(Rmax-Rmin)*(Smax-Smin);
 177                        for
 178                          (Int_t isector=1;
 179                           isector<=NumberOfSectors[ivol-1];
 180                           isector++)
 181                          if(PadADC[ivol-1][isector-1][iring-1]<zeroADC)
 182                            ZeroPads++;
 183                      } //ring
 184                    /*
 185                    cout<<"\nRmin="<<Rmin;
 186                    cout<<"\nRmax="<<Rmax;
 187                    cout<<"\nSmin="<<Smin;
 188                    cout<<"\nSmax="<<Smax;
 189
 190                    cout<<"\nvolume "<<ivol<<" zero "<<ZeroPads<<
 191                      " NumberOfPads "<<NumberOfPads<<endl;
 192                    */
 193                    Float_t zerosRatio=
 194                      (Float_t)ZeroPads/(Float_t)NumberOfPads;
 195                    //  cout<<"\nzerosRatio="<<zerosRatio;
 196                    RecParticles[0]=ivol;
 197                    RecParticles[1]=Smin;
 198                    RecParticles[2]=Smax;
 199                    RecParticles[3]=Rmin;
 200                    RecParticles[4]=Rmax;
 201                    if (zerosRatio>0.1 ||(ivol==2||ivol==4))
 202                      RecParticles[5]=
 203                        Determination_by_Poisson
 204                        (PadADC,ivol,Rmin,Rmax,Smin,Smax);
 205                    else
 206                      RecParticles[5]=
 207                        Determination_by_thresholds
 208                        (PadADC,ivol,Rmin,Rmax,Smin,Smax);
 209                    //  cout<<" nDeterm "<<RecParticles[5]<<endl;
 210                    new((*fReconParticles)[nRecPart++])
 211                  AliFMDReconstParticles(RecParticles);
 212         } // eta
 213     } // volume
 214
 215            // if(zerosRatio<0.01) Determination_by_counting (ZeroPads) ;
 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::Determination_by_thresholds
 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   Int_t thresholdInner[30]={
 231     0,     14,  28,    42,   57,
 232     72,    89,  104,  124,  129,
 233     164,  174,  179,  208,  228,
 234     248,  268,   317,  337,  357,
 235     392,  407,  416,  426,  436,
 236     461,  468,  493,  506,  515};
 237
 238   Int_t thresholdOuter[30]={0, 18, 48, 77, 105,
 239                             132, 165, 198, 231,
 240                             264, 286, 308, 334,
 241                             352, 374, 418, 440,
 242                             462, 484, 506, 528,
 243                             550, 572, 594, 616};
 244   Int_t threshold[30];
 245   for (Int_t it=0; it<30; it++) {
 246     if(volume==1||volume==3||volume==5) threshold[it]=thresholdInner[it];
 247     if(volume==2||volume==4) threshold[it]=thresholdOuter[it];
 248   }
 249   Int_t threshold_array_size = 30;
 250   Int_t NumPart=0;
 251   //Inner Si
 252   for (Int_t iring=Rmin; iring<Rmax; iring++)
 253     {
 254       for (Int_t isector=Smin; isector<Smax; isector++)
 255         {
 256           for (int i=0;i<threshold_array_size-1;i++)
 257             {
 258               if(PadADC[volume-1][isector][iring]>threshold[i]
 259                  &&PadADC[volume-1][isector][iring]<=threshold[i+1])
 260                 {
 261                   NumPart+=i;
 262                   break;
 263                 }; // if in threshol interval
 264             }; //threshold_array_size
 265         }; //iring
 266     }; //sector
 267   // cout<<"\nEnd threshol method"<<endl;
 268   return NumPart;
 269 }
 270
 271
 272  //__________________________________________________________________
 273
 274 Int_t AliFMDReconstruction::Determination_by_Poisson (Int_t PadADC[][50][800],
 275                                            Int_t vol, Int_t rmin, Int_t rmax,
 276                                            Int_t secmin, Int_t secmax)
 277 {
 278   Int_t threshold_adc=1;
 279   Int_t zeropads=0;
 280   for (Int_t i=rmin;i<rmax;i++)
 281     {
 282       for (Int_t j=secmin;j<secmax;j++)
 283         {
 284           if (PadADC[vol-1][j][i]<threshold_adc) zeropads++;
 285         };
 286     };
 287   Float_t lambda=-TMath::Log(Float_t(zeropads)/
 288                              ( Float_t(secmax-secmin)*
 289                                Float_t(rmax-rmin))); //+1 zdes' ne nado
 290   Int_t Recon=(Int_t)(lambda*(secmax-secmin)*(rmax-rmin)); //+1 zdes' ne nado
 291   //  cout<<"\nEnd Poisson method"<<endl;
 292   return Recon;
 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
 318
 319
 320
 321
 322
 323