Reconstruction in 0.1 eta over all sectors
authoralla <alla@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 29 Nov 2002 12:22:14 +0000 (12:22 +0000)
committeralla <alla@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 29 Nov 2002 12:22:14 +0000 (12:22 +0000)
FMD/AliFMD.cxx
FMD/AliFMDReconstParticles.cxx
FMD/AliFMDReconstParticles.h
FMD/AliFMDReconstruction.cxx
FMD/AliFMDReconstruction.h

index f1cd074..6094ab6 100644 (file)
@@ -36,6 +36,7 @@
 //                                                                           //
 ///////////////////////////////////////////////////////////////////////////////
 
+
 #define DEBUG
 #include <TMath.h>
 #include <TGeometry.h>
@@ -57,7 +58,6 @@
 #include "AliFMDReconstParticles.h"
 #include <stdlib.h>
 
-
 ClassImp (AliFMD)
   //_____________________________________________________________________________
 AliFMD::AliFMD ():AliDetector ()
@@ -231,6 +231,7 @@ void AliFMD::ResetDigits ()
 }
 
 //-------------------------------------------------------------------------
+
 void  AliFMD::Init ()
 {
   //
@@ -258,8 +259,10 @@ void  AliFMD::Init ()
     fIdSens2 = gMC->VolId ("GRN2");    //Si sensetive volume
     fIdSens3 = gMC->VolId ("GRN3");    //Si sensetive volume
     fIdSens4 = gMC->VolId ("GRN4");    //Si sensetive volume
+    fIdSens5 = gMC->VolId ("GRN5");    //Si sensetive volume
 
 }
+
 //---------------------------------------------------------------------
 void AliFMD::MakeBranch (Option_t * option, const char *file)
 {
@@ -337,7 +340,7 @@ void AliFMD::SetTreeAddress ()
 void AliFMD::SetRingsSi1(Int_t ringsSi1)
 {
   //  fRingsSi1=ringsSi1;
-  fRingsSi1=256;
+  fRingsSi1=768;
 }
 void AliFMD::SetSectorsSi1(Int_t sectorsSi1)
 {
@@ -345,7 +348,7 @@ void AliFMD::SetSectorsSi1(Int_t sectorsSi1)
 }
 void AliFMD::SetRingsSi2(Int_t ringsSi2)
 {
-  fRingsSi2=128;
+  fRingsSi2=384;
 }
 void AliFMD::SetSectorsSi2(Int_t sectorsSi2)
 {
@@ -434,7 +437,6 @@ void AliFMD::Digits2Reco()
   char * fileHeader=0;
   AliFMDReconstruction * reconstruction =
     new AliFMDReconstruction(fileHeader,fileReconParticles) ;
-  //  fReconParticles=new TClonesArray("AliFMDReconstParticles",1000);
   reconstruction->Exec("");
   delete  reconstruction;
 }
index 9a183f4..42bb482 100644 (file)
@@ -5,19 +5,25 @@ ClassImp(AliFMDReconstParticles)
 AliFMDReconstParticles::AliFMDReconstParticles()
 {
      fNumOfDet=0;
-     fNumOfSector=0;
-     fNumOfRing=0;
+     fNumOfMinSector=0;
+     fNumOfMaxSector=0;
+     fNumOfMinRing=0;
+     fNumOfMaxRing=0;
      fNumOfReconstParticles=0;
 }
 AliFMDReconstParticles::AliFMDReconstParticles(Int_t *RecParticles)
 {
      fNumOfDet=RecParticles[0];
-     fNumOfSector=RecParticles[1];
-     fNumOfRing=RecParticles[2];
-     fNumOfReconstParticles=RecParticles[3];
+     fNumOfMinSector=RecParticles[1];
+     fNumOfMaxSector=RecParticles[2];
+     fNumOfMinRing=RecParticles[3];
+     fNumOfMaxRing=RecParticles[4];
+     fNumOfReconstParticles=RecParticles[5];
 }
 Int_t AliFMDReconstParticles::GetVolume(){return fNumOfDet;}
-Int_t AliFMDReconstParticles::GetNumberOfSector() {return fNumOfSector;}
-Int_t AliFMDReconstParticles::GetNumberOfRing() {return fNumOfRing;}
+Int_t AliFMDReconstParticles::GetNumberOfMinSector() {return fNumOfMinSector;}
+Int_t AliFMDReconstParticles::GetNumberOfMaxSector() {return fNumOfMaxSector;}
+Int_t AliFMDReconstParticles::GetNumberOfMinRing() {return fNumOfMinRing;}
+Int_t AliFMDReconstParticles::GetNumberOfMaxRing() {return fNumOfMaxRing;}
 Int_t AliFMDReconstParticles::GetNumberOfReconstParticles() {return fNumOfReconstParticles;}
 
index c203614..1f1f4ec 100644 (file)
@@ -8,8 +8,10 @@ class AliFMDReconstParticles: public TObject
  //Reconstracted Particles Class
  private:
   Int_t fNumOfDet;                       //Number of FMD disk;
-  Int_t fNumOfSector;                    //Number of sector
-  Int_t fNumOfRing;                      //Number of ring
+  Int_t fNumOfMinSector;                    //Number of min sector 
+  Int_t fNumOfMaxSector;                    //Number of max sector
+  Int_t fNumOfMinRing;                   //Number of min ring
+  Int_t fNumOfMaxRing;                      //Number of max ring
   Int_t fNumOfReconstParticles;          //Number of reconstructed particles
 
  public:
@@ -17,9 +19,11 @@ class AliFMDReconstParticles: public TObject
   AliFMDReconstParticles (Int_t *RecParticles);
   virtual ~AliFMDReconstParticles(){};
   Int_t GetVolume();                    //Return the number of volume 
-  Int_t GetNumberOfSector();            //Return the number of sector
-  Int_t GetNumberOfRing();              //Return the number of ring 
-  Int_t GetNumberOfReconstParticles();  //Returnthe number of reconstructed particles
-  ClassDef(AliFMDReconstParticles,1)
+  Int_t GetNumberOfMinSector();           //Return the number of min sector
+  Int_t GetNumberOfMaxSector();           //Return the number of max sector
+  Int_t GetNumberOfMinRing();             //Return the number of min ring
+  Int_t GetNumberOfMaxRing();              //Return the number of max ring 
+  Int_t GetNumberOfReconstParticles();  //Returnthe the number of reconstructed particles
+  ClassDef(AliFMDReconstParticles,2)
 };
 #endif
index db64cc0..49480ab 100644 (file)
 #include "TFile.h"
 #include "TROOT.h"
 #include "TFolder.h"
+#include <Riostream.h>
 
 // --- Standard library ---
 #include <stdlib.h>
-#include <Riostream.h>
 
 // --- AliRoot header files ---
 
@@ -40,7 +40,6 @@
 #include "AliFMDv1.h"
 #include "AliFMDReconstruction.h"
 #include "AliRun.h"
-
 ClassImp(AliFMDReconstruction)
 
         
@@ -69,102 +68,236 @@ AliFMDReconstruction::AliFMDReconstruction(char* HeaderFile, char *ReconstPartic
 
 AliFMDReconstruction::~AliFMDReconstruction()
 {
 }
-
-//____________________________________________________________________________
+//----------------------------------------------------------------------------
 
 void AliFMDReconstruction::Exec(Option_t *option) 
 { 
+  printf (" AliFMDReconstruction starting \n");
  //Collects all digits in the same active volume into number of particles
-
-  AliFMD * FMD = (AliFMD *) gAlice->GetDetector("FMD");
+  
+  Int_t const NumVolums=5;
+  Int_t const NumSectors=40;
+  Int_t const NumRings=768;
+  Int_t PadADC[10][50][800];
+  Int_t ivol, iSector, iRing;
+  Int_t Ne1; 
+  //Int_t NumberOfRings[5]=  {256,128,256,128,256};
+  Int_t NumberOfSectors[5]=  {20,40,20,40,20};
+  // number of ring for boundary 0.1 eta
+   Int_t EtaIntervalInner []=
+      {0, 55, 110, 165, 221, 276, 331, 386, 442,
+       497, 552, 607, 663, 718, 767 };
+   /*
+     {0,  18, 36,   55,  73,
+      92, 110, 128, 147, 165,
+      184, 202, 221, 239, 255};*/
+   Int_t  EtaIntervalOuter []=  //{0, 21, 43, 65, 86, 108, 127};
+     {0, 65, 130, 195, 260, 325, 383};
+ Int_t EtaInterval[20];
+  
+
+  Int_t size_EtaIntervalInner=sizeof (EtaIntervalInner)/sizeof(EtaIntervalInner[0]);
+
+  Int_t size_EtaIntervalOuter=sizeof (EtaIntervalOuter)/sizeof(EtaIntervalOuter[0]);
+
+  AliFMD * FMD = (AliFMD *) gAlice->GetDetector("FMD");  
   TClonesArray *fReconParticles=FMD->ReconParticles();
-  if(fNevents == 0) fNevents=(Int_t)gAlice->TreeD()->GetEntries(); 
+  if(fNevents == 0) fNevents=(Int_t)gAlice->TreeE()->GetEntries(); 
   for(Int_t ievent=0;ievent<fNevents;ievent++)
     { 
+      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 ???
       gAlice->GetEvent(ievent) ;
       if(gAlice->TreeH()==0) return; 
       if(gAlice->TreeR()==0) gAlice->MakeTree("R");
       //Make branches
+      AliFMDdigit  *fmdDigit;
       FMD->MakeBranch("R");
       
-      Int_t threshold[]={ 0,     14,  28,    42,   57, 
-                         72,    89,  104,  124,  129, 
-                         164,  174,  179,  208,  228, 
-                         248,  268,   317,  337,  357, 
-                         392,  407,  416,  426,  436, 
-                         461,  468,  493,  506,  515, 
-                         541,  566,  592,  617,  642, 
-                         668,  693,  719,  744,  770, 
-                         795,  821,  846,  871,  897, 
-                         922,  948,  973,  999, 1024};
-      
+      Int_t zeroADC=1;
+      // Int_t  threshold_array_size=30;
 
-      /*
-       Int_t threshold[]={    0,  18,  37,  56,   76, 
-                             96, 119, 138, 165,  172, 
-                                    218, 231, 238, 277,  304, 
-                            330, 357, 423, 449,  476, 
-                            522, 542, 555, 568,  581, 
-                             614, 624, 657, 674,  687, 
-                            700, 713, 720, 727,  733, 
-                            740, 759, 778, 797,  816, 
-                            834, 853, 872, 891,  910, 
-                            929, 948, 967, 986,  1024};
-      */
-      
-      int threshold_array_size=sizeof(threshold)/sizeof(threshold[0]);
-      AliFMDdigit  *fmdDigit;
-      //    cout<<" AliFMDdigit "<<AliFMDdigit<<endl;
+      // cout<<" AliFMDdigit "<<AliFMDdigit<<endl;
       if (FMD)
        {
          gAlice->TreeD()->GetEvent(0); 
          TClonesArray *FMDdigits=FMD->Digits();
          Int_t nDigits=FMDdigits->GetEntries();
-          Int_t RecParticles[4];
+          Int_t RecParticles[6];
           Int_t nRecPart=0 ;
+          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);    
-             RecParticles[0] = fmdDigit->Volume();
-             RecParticles[1] = fmdDigit->NumberOfSector();
-             RecParticles[2] = fmdDigit->NumberOfRing();
-             Int_t ADC=fmdDigit->ADCsignal(); 
-             RecParticles[3]=0; //case when fmdDigit->ADCsignal()==0
-             for (int i=0;i<threshold_array_size-1;i++)
-                   {
-                     if(ADC>threshold[i]&&ADC<=threshold[i+1])
-                        RecParticles[3]=i;              
-                   }
-             new((*fReconParticles)[nRecPart++]) AliFMDReconstParticles(RecParticles); 
+              fmdDigit=(AliFMDdigit*)FMDdigits->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]=
+                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;
             } //digit loop
-        }//if FMD
-       gAlice->TreeR()->Reset();
-       gAlice->TreeR()->Fill(); 
-       gAlice->TreeR()->Write(0,TObject::kOverwrite);
+          Int_t Rmin=0; Int_t Rmax=0; //To avoid warning
+          Int_t Smin=0; Int_t Smax=0; //To avoid warning
+                  
+          for (ivol=1; ivol<=NumVolums; ivol++)
+            {
+              if (ivol==1||ivol==3||ivol==5)
+                {
+                  Ne1=size_EtaIntervalInner;
+                  for(Int_t ieta=0; ieta<Ne1; ieta++)
+                    EtaInterval[ieta]=EtaIntervalInner[ieta];
+                }
+              if (ivol==2||ivol==4)
+                {
+                  Ne1=size_EtaIntervalOuter;
+                  for( Int_t ieta=0; ieta<Ne1; ieta++)
+                    EtaInterval[ieta]=EtaIntervalOuter[ieta];
+                }
+
+
+                  for (Int_t e1=0;e1<Ne1-1;e1++) // vol<=NumVolums
+                    {
+                      Rmin=EtaInterval[e1];
+                  Rmax=EtaInterval[e1+1];
+                  ZeroPads=0;
+                  Smin=0;
+                  Smax=NumberOfSectors[ivol-1]; 
+                  for (Int_t iring=Rmin; iring<Rmax; iring++) 
+                    {
+                      NumberOfPads=(Rmax-Rmin)*(Smax-Smin);
+                      for 
+                        (Int_t isector=1;
+                         isector<=NumberOfSectors[ivol-1]; 
+                         isector++) 
+                        if(PadADC[ivol-1][isector-1][iring-1]<zeroADC)
+                          ZeroPads++;
+                    } //ring
+                  /*
+                  cout<<"\nRmin="<<Rmin;
+                  cout<<"\nRmax="<<Rmax;
+                  cout<<"\nSmin="<<Smin;
+                  cout<<"\nSmax="<<Smax;
+              
+                  cout<<"\nvolume "<<ivol<<" zero "<<ZeroPads<<
+                    " NumberOfPads "<<NumberOfPads<<endl;
+                  */
+                  Float_t zerosRatio= 
+                    (Float_t)ZeroPads/(Float_t)NumberOfPads;
+                  //  cout<<"\nzerosRatio="<<zerosRatio;                     
+                  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,Rmin,Rmax,Smin,Smax);
+                  else
+                    RecParticles[5]=
+                      Determination_by_thresholds
+                      (PadADC,ivol,Rmin,Rmax,Smin,Smax);
+                  //  cout<<" nDeterm "<<RecParticles[5]<<endl;
+                  new((*fReconParticles)[nRecPart++]) 
+                AliFMDReconstParticles(RecParticles);                 
+       } // eta
+    } // volume
+          
+          // if(zerosRatio<0.01) Determination_by_counting (ZeroPads) ;           
+       }//if FMD
+      gAlice->TreeR()->Reset();
+      gAlice->TreeR()->Fill(); 
+      gAlice->TreeR()->Write(0,TObject::kOverwrite);
     } //event loop
-  cout<<"\nAliFMDReconstruction::Exec finished"<<endl;
-}
-//__________________________________________________________________
+  // cout<<"\nAliFMDReconstruction::Exec finished"<<endl;
+};
 
-void AliFMDReconstruction::SetReconstParticlesFile(char * file )
+//------------------------------------------------------------
+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)
 {
-   if (!fReconstParticlesFile.IsNull()) 
-    cout<<"\nChanging reconstructed particles file from "<<
-      (char *) fReconstParticlesFile.Data()<< " to "<<file<<endl;
-    fReconstParticlesFile=file;
+  
+  Int_t thresholdInner[30]={
+    0,     14,  28,    42,   57,     
+    72,    89,  104,  124,  129, 
+    164,  174,  179,  208,  228, 
+    248,  268,   317,  337,  357, 
+    392,  407,  416,  426,  436, 
+    461,  468,  493,  506,  515}; 
+
+  Int_t thresholdOuter[30]={0, 18, 48, 77, 105,
+                           132, 165, 198, 231, 
+                           264, 286, 308, 334, 
+                           352, 374, 418, 440,
+                           462, 484, 506, 528,
+                           550, 572, 594, 616}; 
+  Int_t threshold[30];
+  for (Int_t it=0; it<30; it++) {
+    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;
+  //Inner Si
+  for (Int_t iring=Rmin; iring<Rmax; iring++) 
+    {
+      for (Int_t isector=Smin; isector<Smax; isector++) 
+       {
+         for (int i=0;i<threshold_array_size-1;i++)
+           {         
+             if(PadADC[volume-1][isector][iring]>threshold[i]
+                &&PadADC[volume-1][isector][iring]<=threshold[i+1])
+               {
+                 NumPart+=i;
+                 break;
+               }; // if in threshol interval
+           }; //threshold_array_size
+       }; //iring
+    }; //sector
+  // cout<<"\nEnd threshol method"<<endl;
+  return NumPart;
 }
-//__________________________________________________________________
-void AliFMDReconstruction::Print(Option_t* option)const
+
+
+ //__________________________________________________________________
+
+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)
 {
-  cout<<"------------------- "<<GetName() <<" -------------"<< endl ;
-  if(fReconstParticlesFile.IsNull())
-    cout<<"\nWriting reconstructed particles to file"<<endl ;
-  else
-    cout<<"\nWriting reconstructed particles to file  "<<
-      (char*) fReconstParticlesFile.Data() << endl ;
-}
+  Int_t threshold_adc=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++;
+       };
+    };
+  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
+  //  cout<<"\nEnd Poisson method"<<endl;
+  return Recon;
+};
+
+//------------------------------------------------------------------
+
+
+
+
+
 
 
 
index 42ca010..268b81c 100644 (file)
@@ -22,14 +22,21 @@ class AliFMDReconstruction: public TTask
   virtual void  Exec(Option_t *option); 
   void SetNEvents(Int_t Nevents){fNevents = Nevents;}
   Stat_t GetNEvents(){return fNevents;}
-  void SetReconstParticlesFile(char * file ) ;
-  virtual void Print(Option_t* option) const ;   
+   TClonesArray *Digits() const {return fDigits;}
+   Int_t Determination_by_thresholds(Int_t a[10][50][800], Int_t volume, Int_t Rmin, Int_t Rmax, 
+                                    Int_t Smin, Int_t Smax);
+   Int_t Determination_by_Poisson (Int_t PadADC[10][50][800], Int_t, Int_t, Int_t, Int_t, Int_t);
 
  private:
+  TClonesArray *fDigits;               // ! array with digits
   Int_t   fNevents ;                         // Number of events
   TString fReconstParticlesFile;             //output file 
   TString fHeadersFile ;                     //input file
+
+
   ClassDef(AliFMDReconstruction,2) 
+
+
 };
 #endif