f98f6dc7f9def856c65ef2c14227f6d499599a73
[u/mrichter/AliRoot.git] / 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 "TH2F.h"
31
32 // --- Standard library ---
33 #include <stdlib.h>
34 #include <Riostream.h>
35
36 // --- AliRoot header files ---
37
38 #include "AliRunLoader.h"
39 #include "AliLoader.h"
40
41 #include "AliFMDdigit.h"
42 #include "AliFMDhit.h"
43 #include "AliFMDReconstParticles.h"
44 #include "AliFMD.h"
45 #include "AliFMDv1.h"
46 #include "AliFMDReconstruction.h"
47 #include "AliRun.h"
48 #include "AliConfig.h"
49 #include "AliHeader.h"
50 #include "AliGenEventHeader.h"
51
52 ClassImp(AliFMDReconstruction)
53
54         
55 //____________________________________________________________________________ 
56
57 AliFMDReconstruction::AliFMDReconstruction():TTask("AliFMDReconstruction","") 
58 {
59   fNevents = 0 ;  // Number of events to rreconnstraction, 0 means all events in current file
60   fRunLoader = 0x0;
61   
62 }
63 //____________________________________________________________________________ 
64
65 AliFMDReconstruction::AliFMDReconstruction(AliRunLoader* rl):TTask("AliFMDReconstruction","")
66 {
67
68   if (rl == 0x0)
69    {
70      Fatal("AliFMDReconstruction","Argument AliRunLoader* is null!");
71      return;
72    }
73    
74   fNevents = 0 ;    // Number of events to rreconnstraction, 0 means all events in current file
75
76   fRunLoader = rl;
77   AliLoader* gime = fRunLoader->GetLoader("FMDLoader");
78   if (gime == 0x0)
79    {
80      Fatal("AliFMDReconstruction","Can not find FMD (loader) in specified event");
81      return;//never reached
82    }
83   //add Task to //root/Tasks folder
84   gime->PostReconstructioner(this);
85 }
86
87 //____________________________________________________________________________ 
88
89 AliFMDReconstruction::~AliFMDReconstruction()
90 {
91   AliLoader* gime = fRunLoader->GetLoader("FMDLoader");
92   gime->CleanReconstructioner();
93 }
94
95 //____________________________________________________________________________
96
97 void AliFMDReconstruction::Exec() 
98
99  //Collects all digits in the same active volume into number of particles
100   /*
101     Reconstruct number of particles 
102     in given group of pads for given FMDvolume
103     determine by numberOfVolume , 
104     numberOfMinSector,numberOfMaxSector,
105     numberOfMinRing, numberOgMaxRing
106     Reconstruction method choose dependence on number of empty pads  
107   */
108
109
110   cout<<"\nStart AliFMDReconstruction::Exec(...)"<<endl;
111   Int_t const knumVolumes=5;
112   Int_t padADC;
113   Float_t eta, etain,etaout,rad,theta;
114   Int_t ivol, iSector, iRing;
115   Float_t rin[5]={4.2,15.4,4.2,15.4,4.2};
116   Float_t rout[5]={17.4,28.4,17.4,28.4,17.4};
117   Float_t z[5]={-62.8, -75.2, 83.4, 75.2, 340.};
118   Int_t numberOfRings[5]={512,256,512,256,512};
119   Int_t numberOfSectors[5]=  {20,40,20,40,20};
120   Int_t numberOfEtaIntervals[5];
121   // number of ring for boundary 0.1 eta
122
123   
124   if (fRunLoader == 0x0)
125    {
126     Error("Exec","Run Loader loader is NULL - Session not opened");
127     return;
128    }
129  
130  
131  AliLoader* gime = fRunLoader->GetLoader("FMDLoader");
132   if (gime == 0x0)
133    {
134      Fatal("AliFMDReconstruction","Can not find FMD (loader) in specified event");
135      return;//never reached
136    }
137    
138   if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
139   if (!fRunLoader->TreeE()) fRunLoader->LoadHeader();
140   Int_t retval;
141   TDirectory* cwd = gDirectory;
142   gDirectory = 0x0;
143     Text_t buf1[20];
144   TH2F* hTotal[10];
145   for (Int_t j=1; j<=5; j++){
146     sprintf(buf1,"hTotal%d",j);
147     
148     hTotal[j] = new TH2F(buf1," Number of primary particles ",
149                          numberOfSectors[j-1],1,numberOfSectors[j-1],
150                          numberOfRings[j-1],1,numberOfRings[j-1]);
151   }   
152   gDirectory = cwd;
153  
154        
155   if(fNevents == 0) fNevents=Int_t (fRunLoader->TreeE()->GetEntries()); 
156   cout<<" fNevents "<<fNevents<<endl;
157    for(Int_t ievent=0;ievent<fNevents;ievent++)
158     { 
159       fRunLoader->GetEvent(ievent) ;
160
161       cout<<" ievent "<<ievent<<endl;
162     for (Int_t i=0; i<5; i++)
163         hTotal[i+1]->Reset();
164       
165       retval = gime->LoadDigits("READ"); 
166       if (retval)
167         {
168           Error("Exec","Error occured while loading digits. Exiting.");
169           return;
170         }
171       
172       AliFMD * fFMD = (AliFMD *) gAlice->GetDetector("FMD");
173       TClonesArray *fReconParticles=fFMD->ReconParticles();
174       TClonesArray *fDigits=fFMD->Digits();
175  
176       TTree* treeD = gime->TreeD();
177       if (treeD == 0x0)
178         {
179           Error("Exec","Can not get Tree with Digits. Nothing to reconstruct - Exiting");
180           return;
181         }
182       
183       TBranch *brDigits=0;
184             
185       brDigits=treeD->GetBranch("FMD");
186       cout<<" brDigits "<<brDigits<<endl;
187       if (brDigits) {
188         brDigits->SetAddress(&fDigits);
189         //        fFMD->SetHitsAddressBranch(brHits);
190       }else{
191         cerr<<"EXEC Branch FMD digits not found"<<endl;
192         return;
193       } 
194           
195       if(gime->TreeR()==0) gime->MakeTree("R");
196
197       //Make branches
198       fFMD->MakeBranch("R");
199
200       
201       Int_t zeroADC=6;
202  
203       AliFMDdigit  *fmdDigit;
204        if (fFMD)
205         {
206           gime->TreeD()->GetEvent(0); 
207                   
208           Int_t nDigits=fDigits->GetEntries();
209           cout<<" nDigits "<<nDigits<<endl;
210           Int_t recParticles[6];
211           Int_t nRecPart=0 ;
212            Int_t zeroPads=0;
213            Int_t numberOfPads=0; //To avoid warning
214            Int_t pedestal;
215            Float_t channelWidth=(22400*50)/1024;
216            for (Int_t digit=0;digit<nDigits;digit++) 
217              {
218                fmdDigit=(AliFMDdigit*)fDigits->UncheckedAt(digit);    
219                ivol=fmdDigit->Volume();
220                iSector=fmdDigit->NumberOfSector();
221                iRing=fmdDigit->NumberOfRing();
222                pedestal=Int_t(gRandom->Gaus(500,250));
223                padADC= fmdDigit->ADCsignal()
224                  -Int_t(Float_t(pedestal)/channelWidth);
225                if (padADC<0) padADC=0;
226               hTotal[ivol]->Fill(iSector,iRing,padADC);
227              } //digit loop
228            Int_t rmin=0; Int_t rmax=0; //To avoid warning
229            Int_t smin=0; Int_t smax=0; //To avoid warning
230            AliHeader *header = fRunLoader->GetHeader();
231            AliGenEventHeader* genHeader = header->GenEventHeader();
232            TArrayF *o = new TArrayF(3); 
233            genHeader->PrimaryVertex(*o);
234            Float_t zVertex=o->At(2);
235            for (ivol=0; ivol<knumVolumes; ivol++)
236              {
237                Float_t ring2number=Float_t (numberOfRings[ivol])/
238                  (rout[ivol]-rin[ivol]);
239                Float_t realZ=zVertex+z[ivol];
240                theta=TMath::ATan(rout[ivol]/TMath::Abs(realZ));
241                etain = - TMath::Log( TMath::Tan(theta/2.));
242                theta=TMath::ATan(rin[ivol]/TMath::Abs(realZ));
243                etaout=- TMath::Log( TMath::Tan(theta/2.));
244                numberOfEtaIntervals[ivol]=Int_t((etaout-etain)*10)-1;
245                eta=etain;
246                for (Int_t e1=0;e1<=numberOfEtaIntervals[ivol];e1++) 
247                  {
248                    theta = 2.*TMath::ATan (TMath::Exp (-eta));
249                    Float_t radmin = TMath::Abs(realZ) * (TMath::Tan (theta));
250                    rmax= Int_t ( (radmin-rin[ivol])*ring2number+0.5);
251                    cout<<ivol<<" "<<" eta "<<eta<<" radmin "<<radmin<<
252
253                      " Rmax "<<rmax<<" "<<endl;;
254                    eta=eta+0.1;
255                    theta = 2.*TMath::ATan (TMath::Exp (-eta));
256                    rad = TMath::Abs(realZ) * (TMath::Tan (theta));
257                    rmin=Int_t( (rad-rin[ivol])*ring2number+0.5);
258                    cout<<"eta "<<eta<<" rad "<<rad<<" Rmin "<<rmin<<endl;
259                    //      if(ivol==2&&e1==13) rmin=0;
260                    zeroPads=0;
261                    smin=0;
262                    smax=numberOfSectors[ivol]; 
263                    numberOfPads=(rmax-rmin)*(smax-smin);
264                    for (Int_t iring=rmin; iring<rmax; iring++) 
265                      {
266                        for 
267                          (Int_t isector=0;
268                           isector<numberOfSectors[ivol]; 
269                           isector++) 
270                          {
271                            if(hTotal[ivol+1]->GetBinContent(isector+1,iring+1)
272                             <zeroADC) zeroPads++;}
273
274                      } //ring
275                    Float_t numberOfPads=Float_t(smax-smin)*Float_t(rmax-rmin);
276
277                    cout<<" zero "<<zeroPads++<<" pads "<<numberOfPads;
278                    Double_t lambda=-TMath::Log(Double_t(zeroPads)/numberOfPads);
279                    Int_t fRecon=Int_t (lambda*numberOfPads+0.5);
280                    
281                    Float_t zerosRatio= 
282                      (Float_t)zeroPads/(Float_t)numberOfPads;
283                    cout<<" zerosRatio "<<zerosRatio<<" recon "<<fRecon<<endl;
284                    recParticles[0]=ivol+1;
285                    recParticles[1]=smin;
286                    recParticles[2]=smax;
287                    recParticles[3]=rmin;
288                    recParticles[4]=rmax;
289                    recParticles[5]= fRecon;
290                    new((*fReconParticles)[nRecPart++]) 
291                      AliFMDReconstParticles(recParticles);             
292
293
294                  } // eta
295              } // volume
296            
297         }//if FMD
298        gime->TreeR()->Reset();
299        gime->TreeR()->Fill(); 
300        gime->WriteRecPoints("OVERWRITE");
301     } //event loop
302   cout<<"\nAliFMDReconstruction::Exec finished"<<endl;
303 }
304