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