]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - FMD/AliFMDReconstruction.cxx
volume overlap fixed
[u/mrichter/AliRoot.git] / FMD / AliFMDReconstruction.cxx
index a1ccba883da591c4242aa55836047ea85c72d294..dd5222ee3b80db105644bf8e721608e0cde9e5a3 100644 (file)
@@ -76,34 +76,43 @@ AliFMDReconstruction::~AliFMDReconstruction()
 
 void AliFMDReconstruction::Exec(Option_t *option) 
 { 
- //Collects all digits in the same active volume into number of particles
+  /*
+    Reconstruct nember of particles 
+    in given group of pads for given FMDvolume
+    determine by numberOfVolume , 
+    numberOfMinSector,numberOfMaxSector,
+    numberOfMinRing, numberOgMaxRing
+    Reconstruction method choose dependence on number of empty pads  
+  */
+
+
   cout<<"\nStart AliFMDReconstruction::Exec(...)"<<endl;
-  Int_t const NumVolums=5;
-  Int_t const NumSectors=40;
-  Int_t const NumRings=768;
-  Int_t PadADC[10][50][800];
+  Int_t const knumVolumes=5;
+  Int_t const knumSectors=40;
+  Int_t const knumRings=768;
+  Int_t padADC[10][50][800];
   Float_t eta, etain,etaout,rad,theta;
   Int_t ivol, iSector, iRing;
   Float_t rin[5]={4.2,15.4,4.2,15.4,4.2};
   Float_t rout[5]={17.4,28.4,17.4,28.4,17.4};
   Float_t z[5]={62.8, 75.2, -83.4, -75.2, -340.};
-  Int_t NumberOfRings[5]={768,384,768,384,768};
-  Int_t NumberOfSectors[5]=  {20,40,20,40,20};
-  Int_t NumberOfEtaIntervals[5];
+  Int_t numberOfRings[5]={768,384,768,384,768};
+  Int_t numberOfSectors[5]=  {20,40,20,40,20};
+  Int_t numberOfEtaIntervals[5];
   // number of ring for boundary 0.1 eta
   TBranch *brDigits=0;
 
-  AliFMD * FMD = (AliFMD *) gAlice->GetDetector("FMD");
-  TClonesArray *fReconParticles=FMD->ReconParticles();
+  AliFMD * fFMD = (AliFMD *) gAlice->GetDetector("FMD");
+  TClonesArray *fReconParticles=fFMD->ReconParticles();
   if(fNevents == 0) fNevents=(Int_t)gAlice->TreeE()->GetEntries(); 
   cout<<" fNevents "<<fNevents<<endl;
   for(Int_t ievent=0;ievent<fNevents;ievent++)
     { 
       cout<<" ievent "<<ievent<<endl;
-      for (Int_t i=0; i<NumVolums; i++)
-       for(Int_t j=0; j<NumSectors; j++)
-         for(Int_t ij=0; ij<NumRings; ij++)
-           PadADC[i][j][ij]=0;                    //zhachem ???
+      for (Int_t i=0; i<knumVolumes; i++)
+       for(Int_t j=0; j<knumSectors; j++)
+         for(Int_t ij=0; ij<knumRings; ij++)
+           padADC[i][j][ij]=0;                    //zhachem ???
       gAlice->GetEvent(ievent) ;
       if(gAlice->TreeH()==0) return; 
       if(gAlice->TreeD()==0) return;
@@ -115,92 +124,92 @@ void AliFMDReconstruction::Exec(Option_t *option)
  
       if(gAlice->TreeR()==0) gAlice->MakeTree("R");
       //Make branches
-      FMD->MakeBranch("R");
+      fFMD->MakeBranch("R");
 
       
       Int_t zeroADC=1;
  
       AliFMDdigit  *fmdDigit;
-       if (FMD)
+       if (fFMD)
        {
          gAlice->TreeD()->GetEvent(0); 
-         TClonesArray *FMDdigits=FMD->Digits();
-         Int_t nDigits=FMDdigits->GetEntries();
+         TClonesArray *fFMDdigits=fFMD->Digits();
+         Int_t nDigits=fFMDdigits->GetEntries();
          cout<<" nDigits "<<nDigits<<endl;
-          Int_t RecParticles[6];
+          Int_t recParticles[6];
           Int_t nRecPart=0 ;
-          Int_t ZeroPads=0;
-          Int_t NumberOfPads=0; //To avoid warning
+          Int_t zeroPads=0;
+          Int_t numberOfPads=0; //To avoid warning
           Int_t pedestal;
           Float_t channelWidth=(22400*50)/1024;
           for (Int_t digit=0;digit<nDigits;digit++) 
             {
-              fmdDigit=(AliFMDdigit*)FMDdigits->UncheckedAt(digit);    
+              fmdDigit=(AliFMDdigit*)fFMDdigits->UncheckedAt(digit);    
               ivol=fmdDigit->Volume();
               iSector=fmdDigit->NumberOfSector();
               iRing=fmdDigit->NumberOfRing();
               pedestal=Int_t(gRandom->Gaus(500,250));
-              PadADC[ivol-1][iSector-1][iRing-1]=
+              padADC[ivol-1][iSector-1][iRing-1]=
                 fmdDigit->ADCsignal()
                 -Int_t(Float_t(pedestal)/channelWidth);
-              if (PadADC[ivol-1][iSector-1][iRing-1]<0) 
-                PadADC[ivol-1][iSector-1][iRing-1]=0;
+              if (padADC[ivol-1][iSector-1][iRing-1]<0) 
+                padADC[ivol-1][iSector-1][iRing-1]=0;
             } //digit loop
-          Int_t Rmin=0; Int_t Rmax=0; //To avoid warning
-          Int_t Smin=0; Int_t Smax=0; //To avoid warning
+          Int_t rmin=0; Int_t rmax=0; //To avoid warning
+          Int_t smin=0; Int_t smax=0; //To avoid warning
           AliHeader *header = gAlice->GetHeader();
           AliGenEventHeader* genHeader = header->GenEventHeader();
           TArrayF *o = new TArrayF(3); 
           genHeader->PrimaryVertex(*o);
           Float_t zVertex=o->At(2);
-          for (ivol=0; ivol<NumVolums; ivol++)
+          for (ivol=0; ivol<knumVolumes; ivol++)
             {
               Float_t realZ=zVertex+z[ivol];
               theta=TMath::ATan(rin[ivol]/TMath::Abs(realZ));
               etain = - TMath::Log( TMath::Tan(theta/2.));
               theta=TMath::ATan(rout[ivol]/TMath::Abs(realZ));
               etaout=- TMath::Log( TMath::Tan(theta/2.));
-              NumberOfEtaIntervals[ivol]=Int_t ((etain-etaout)*10);
+              numberOfEtaIntervals[ivol]=Int_t ((etain-etaout)*10);
               eta=etain;
-              for (Int_t e1=0;e1<=NumberOfEtaIntervals[ivol];e1++) 
+              for (Int_t e1=0;e1<=numberOfEtaIntervals[ivol];e1++) 
                 {
                   theta = 2.*TMath::ATan (TMath::Exp (-eta));
                   Float_t radmin = TMath::Abs(realZ) * (TMath::Tan (theta));
-                  Rmin= Int_t ( (radmin-rin[ivol])*NumberOfRings[ivol]/(rout[ivol]-rin[ivol]));
+                  rmin= Int_t ( (radmin-rin[ivol])*numberOfRings[ivol]/(rout[ivol]-rin[ivol]));
                   eta=eta-0.1;
                   theta = 2.*TMath::ATan (TMath::Exp (-eta));
                   rad = TMath::Abs(realZ) * (TMath::Tan (theta));
-                  Rmax=Int_t( (rad-rin[ivol])*NumberOfRings[ivol]/(rout[ivol]-rin[ivol]));
-                  ZeroPads=0;
-                  Smin=0;
-                  Smax=NumberOfSectors[ivol]; 
-                  for (Int_t iring=Rmin; iring<Rmax; iring++) 
+                  rmax=Int_t( (rad-rin[ivol])*numberOfRings[ivol]/(rout[ivol]-rin[ivol]));
+                  zeroPads=0;
+                  smin=0;
+                  smax=numberOfSectors[ivol]; 
+                  for (Int_t iring=rmin; iring<rmax; iring++) 
                     {
-                      NumberOfPads=(Rmax-Rmin)*(Smax-Smin);
+                      numberOfPads=(rmax-rmin)*(smax-smin);
                       for 
                         (Int_t isector=1;
-                         isector<=NumberOfSectors[ivol]; 
+                         isector<=numberOfSectors[ivol]; 
                          isector++) 
-                        if(PadADC[ivol][isector-1][iring-1]<zeroADC)
-                          ZeroPads++;
+                        if(padADC[ivol][isector-1][iring-1]<zeroADC)
+                          zeroPads++;
                     } //ring
                   Float_t zerosRatio= 
-                    (Float_t)ZeroPads/(Float_t)NumberOfPads;
-                  RecParticles[0]=ivol;
-                  RecParticles[1]=Smin;
-                  RecParticles[2]=Smax;
-                  RecParticles[3]=Rmin;
-                  RecParticles[4]=Rmax;
-                  if (zerosRatio>0.1 ||(ivol==2||ivol==4))
-                    RecParticles[5]=
-                      Determination_by_Poisson
-                      (PadADC,ivol+1,Rmin,Rmax,Smin,Smax);
+                    (Float_t)zeroPads/(Float_t)numberOfPads;
+                  recParticles[0]=ivol;
+                  recParticles[1]=smin;
+                  recParticles[2]=smax;
+                  recParticles[3]=rmin;
+                  recParticles[4]=rmax;
+                  if (zerosRatio>0.1 )
+                    recParticles[5]=
+                      DeterminationByPoisson
+                      (padADC,ivol+1,rmin,rmax,smin,smax);
                   else
-                    RecParticles[5]=
-                      Determination_by_thresholds
-                      (PadADC,ivol+1,Rmin,Rmax,Smin,Smax);
+                    recParticles[5]=
+                      DeterminationByThresholds
+                      (padADC,ivol+1,rmin,rmax,smin,smax);
                   new((*fReconParticles)[nRecPart++]) 
-                AliFMDReconstParticles(RecParticles);                 
+                AliFMDReconstParticles(recParticles);                 
                 } // eta
             } // volume
           
@@ -213,10 +222,16 @@ void AliFMDReconstruction::Exec(Option_t *option)
 };
 
 //------------------------------------------------------------
-Int_t AliFMDReconstruction::Determination_by_thresholds
-(Int_t PadADC[][50][800],Int_t volume, Int_t Rmin, Int_t Rmax, 
- Int_t Smin, Int_t Smax)
+Int_t AliFMDReconstruction::DeterminationByThresholds
+(Int_t padAdc[][50][800],Int_t volume, Int_t rmin, Int_t rmax, 
+ Int_t smin, Int_t smax)
 {
+  /*
+    reconstruct number of particles by 
+    energy deposited threshold method 
+    Using if number of empty pads less then 10%
+
+*/
   cout<<"\nStart threshold method\n";
   
   Int_t thresholdInner[30]={
@@ -238,51 +253,59 @@ Int_t AliFMDReconstruction::Determination_by_thresholds
     if(volume==1||volume==3||volume==5) threshold[it]=thresholdInner[it];
     if(volume==2||volume==4) threshold[it]=thresholdOuter[it];
   }
-  Int_t threshold_array_size = 30;     
-  Int_t NumPart=0;
+  Int_t thresholdArraySize = 30;     
+  Int_t numPart=0;
   //Inner Si
-  for (Int_t iring=Rmin; iring<Rmax; iring++) 
+  for (Int_t iring=rmin; iring<rmax; iring++) 
     {
-      for (Int_t isector=Smin; isector<Smax; isector++) 
+      for (Int_t isector=smin; isector<smax; isector++) 
        {
-         for (int i=0;i<threshold_array_size-1;i++)
+         for (int i=0;i<thresholdArraySize-1;i++)
            {         
-             if(PadADC[volume-1][isector][iring]>threshold[i]
-                &&PadADC[volume-1][isector][iring]<=threshold[i+1])
+             if(padAdc[volume-1][isector][iring]>threshold[i]
+                &&padAdc[volume-1][isector][iring]<=threshold[i+1])
                {
-                 NumPart+=i;
+                 numPart+=i;
                  break;
                }; // if in threshol interval
            }; //threshold_array_size
        }; //iring
     }; //sector
   cout<<"\nEnd threshol method"<<endl;
-  return NumPart;
+  return numPart;
 }
 
 
  //__________________________________________________________________
 
-Int_t AliFMDReconstruction::Determination_by_Poisson (Int_t PadADC[][50][800],
-                                          Int_t vol, Int_t rmin, Int_t rmax, 
-                                          Int_t secmin, Int_t secmax)
+Int_t AliFMDReconstruction::DeterminationByPoisson 
+(Int_t padAdc[][50][800],
+ Int_t vol, Int_t rmin, Int_t rmax, 
+ Int_t secmin, Int_t secmax)
 {
+  /* 
+     reconstruct number of particles by Poisson statistic method 
+     according number of empty pad in chosen group of pads 
+     using if number of empty pads more then 10% 
+
+  */
+
   //  cout<<"\nStart Poisson method";
-  Int_t threshold_adc=1;   
+  Int_t thresholdADC=1;   
   Int_t zeropads=0;
   for (Int_t i=rmin;i<rmax;i++)
     {
       for (Int_t j=secmin;j<secmax;j++)
        {
-         if (PadADC[vol-1][j][i]<threshold_adc) zeropads++;
+         if (padAdc[vol-1][j][i]<thresholdADC) zeropads++;
        };
     };
   Float_t lambda=-TMath::Log(Float_t(zeropads)/
                             ( Float_t(secmax-secmin)*
                               Float_t(rmax-rmin))); //+1 zdes' ne nado
-  Int_t Recon=(Int_t)(lambda*(secmax-secmin)*(rmax-rmin)); //+1 zdes' ne nado
+  Int_t fRecon=(Int_t)(lambda*(secmax-secmin)*(rmax-rmin)); //+1 zdes' ne nado
   //  cout<<"\nEnd Poisson method"<<endl;
-  return Recon;
+  return fRecon;
 };
 
 //------------------------------------------------------------------