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