Last minute changes and new code for event mixing and reconstruction (A.Maevskaia)
[u/mrichter/AliRoot.git] / FMD / AliFMD.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 //  Forward Multiplicity Detector based on Silicon plates                    //
18 //  This class contains the base procedures for the Forward Multiplicity     //
19 //  detector                                                                 //
20 //  Detector consists of 6 Si volumes covered pseudorapidity interval         //
21 //  from 1.6 to 6.0.                                                         //
22 //                                                                           //
23 //Begin_Html
24 /*
25 <img src="gif/AliFMDClass.gif">
26 </pre>
27 <br clear=left>
28 <font size=+2 color=red>
29 <p>The responsible person for this module is
30 <a href="mailto:Alla.Maevskaia@cern.ch">Alla Maevskaia</a>.
31 </font>
32 <pre>
33 */
34 //End_Html
35 //                                                                           //
36 //                                                                           //
37 ///////////////////////////////////////////////////////////////////////////////
38
39 #define DEBUG
40 #include <TMath.h>
41 #include <TGeometry.h>
42 #include <TTUBE.h>
43 #include <TTree.h>
44 #include <TNode.h>
45 #include <TFile.h>
46
47 #include <TClonesArray.h>
48 #include <TLorentzVector.h>
49 #include "AliFMDv1.h"
50 #include "AliRun.h"
51 #include "AliMC.h"
52 #include "AliDetector.h"
53 #include <iostream.h>
54 #include <fstream.h>
55 #include "AliMagF.h"
56 #include "AliFMDhit.h"
57 #include "AliFMDdigit.h"
58 #include "AliFMDReconstruction.h"
59 #include <stdlib.h>
60
61
62 ClassImp (AliFMD)
63   //_____________________________________________________________________________
64 AliFMD::AliFMD ():AliDetector ()
65 {
66   //
67   // Default constructor for class AliFMD
68   //
69   fIshunt = 0;
70   fHits     = 0;
71   fDigits   = 0;
72   fSDigits  = 0;
73   fReconParticles=0; 
74 }
75
76 //_____________________________________________________________________________
77 AliFMD::AliFMD (const char *name, const char *title):
78 AliDetector (name, title)
79 {
80   //
81   // Standard constructor for Forward Multiplicity Detector
82   //
83
84   //
85   // Initialise Hit array
86   fHits = new TClonesArray ("AliFMDhit", 1000);
87   // Digits for each Si disk
88   fDigits = new TClonesArray ("AliFMDdigit", 1000);
89   fSDigits = new TClonesArray ("AliFMDdigit", 1000);
90   gAlice->AddHitList (fHits);
91
92   fIshunt = 0;
93   fIdSens1 = 0;
94
95   SetMarkerColor (kRed);
96 }
97
98 //-----------------------------------------------------------------------------
99 AliFMD::~AliFMD ()
100 {
101   if (fHits)
102     {
103       fHits->Delete ();
104       delete fHits;
105       fHits = 0;
106     }
107   if (fDigits)
108     {
109       fDigits->Delete ();
110       delete fDigits;
111       fDigits = 0;
112     }
113   if (fSDigits)
114     {
115       fSDigits->Delete ();
116       delete fSDigits;
117       fSDigits = 0;
118     }
119   if (fReconParticles)
120     {
121       fReconParticles->Delete ();
122       delete fReconParticles;
123       fReconParticles = 0;
124     }
125
126 }
127
128 //_____________________________________________________________________________
129 void AliFMD::AddHit (Int_t track, Int_t * vol, Float_t * hits)
130 {
131   //
132   // Add a hit to the list
133   //
134   TClonesArray & lhits = *fHits;
135   new (lhits[fNhits++]) AliFMDhit (fIshunt, track, vol, hits);
136 }
137
138 //_____________________________________________________________________________
139 void AliFMD::AddDigit (Int_t * digits)
140 {
141   // add a real digit - as coming from data
142
143   // printf("AddDigit\n");
144
145   TClonesArray & ldigits = *fDigits;
146   new (ldigits[fNdigits++]) AliFMDdigit (digits);
147
148 }
149
150 //_____________________________________________________________________________
151 void AliFMD::BuildGeometry ()
152 {
153   //
154   // Build simple ROOT TNode geometry for event display
155   //
156   TNode *node, *top;
157   const int kColorFMD = 7;
158   //
159   top = gAlice->GetGeometry ()->GetNode ("alice");
160
161   // FMD define the different volumes
162   new TRotMatrix ("rot901", "rot901", 90, 0, 90, 90, 180, 0);
163
164   new TTUBE ("S_FMD0", "FMD  volume 0", "void", 4.73, 17.7, 1.5);
165   top->cd ();
166   node = new TNode ("FMD0", "FMD0", "S_FMD0", 0, 0, 64, "");
167   node->SetLineColor (kColorFMD);
168   fNodes->Add (node);
169
170   new TTUBE ("S_FMD1", "FMD  volume 1", "void", 23.4, 36., 1.5);
171   top->cd ();
172   node = new TNode ("FMD1", "FMD1", "S_FMD1", 0, 0, 85, "");
173   node->SetLineColor (kColorFMD);
174   fNodes->Add (node);
175
176   new TTUBE ("S_FMD2", "FMD  volume 2", "void", 4.73, 17.7, 1.5);
177   top->cd ();
178   node = new TNode ("FMD2", "FMD2", "S_FMD2", 0, 0, -64, "");
179   node->SetLineColor (kColorFMD);
180   fNodes->Add (node);
181
182   new TTUBE ("S_FMD3", "FMD  volume 3", "void", 23.4, 36., 1.5);
183   top->cd ();
184   node = new TNode ("FMD3", "FMD3", "S_FMD3", 0, 0, -85, "");
185   node->SetLineColor (kColorFMD);
186   fNodes->Add (node);
187
188   new TTUBE ("S_FMD4", "FMD  volume 4", "void", 5, 15, 0.015);
189   top->cd ();
190   //  node = new TNode("FMD4","FMD4","S_FMD4",0,0,-270,"");
191   node = new TNode ("FMD4", "FMD4", "S_FMD4", 0, 0, -270, "");
192   node->SetLineColor (kColorFMD);
193   fNodes->Add (node);
194
195   /* 
196      new TTUBE("S_FMD5","FMD  volume 5","void",5,14,0.015);
197      top->cd();
198      node = new TNode("FMD5","FMD5","S_FMD5",0,0,-630,"");
199      node->SetLineColor(kColorFMD);
200      fNodes->Add(node);
201    */
202 }
203
204 //_____________________________________________________________________________
205 Int_t AliFMD::DistanceToPrimitive (Int_t px, Int_t py)
206 {
207   //
208   // Calculate the distance from the mouse to the FMD on the screen
209   // Dummy routine
210   //
211   return 9999;
212 }
213
214 //___________________________________________
215 void AliFMD::ResetHits ()
216 {
217   // Reset number of clusters and the cluster array for this detector
218   AliDetector::ResetHits ();
219 }
220
221 //____________________________________________
222 void AliFMD::ResetDigits ()
223 {
224   //
225   // Reset number of digits and the digits array for this detector
226   AliDetector::ResetHits ();
227   //
228 }
229
230 //-------------------------------------------------------------------------
231 void  AliFMD::Init ()
232 {
233   //
234   // Initialis the FMD after it has been built
235   Int_t i;
236   AliMC *pMC = AliMC::GetMC ();
237   //
238   if (fDebug)
239     {
240       printf ("\n%s: ", ClassName ());
241       for (i = 0; i < 35; i++)
242         printf ("*");
243       printf (" FMD_INIT ");
244       for (i = 0; i < 35; i++)
245         printf ("*");
246       printf ("\n%s: ", ClassName ());
247       //
248       // Here the FMD initialisation code (if any!)
249       for (i = 0; i < 80; i++)
250         printf ("*");
251       printf ("\n");
252     }
253   //
254   //
255   if (IsVersion () != 0)
256     fIdSens1 = pMC->VolId ("GRIN");     //Si sensetive volume
257   else
258     fIdSens1 = pMC->VolId ("GFSI");     //Si sensetive volume
259
260 }
261
262 //---------------------------------------------------------------------
263 void AliFMD::MakeBranch (Option_t * option, const char *file)
264 {
265   // Create Tree branches for the FMD.
266   char branchname[10];
267   const Int_t kBufferSize = 16000;
268   sprintf (branchname, "%s", GetName ());
269   
270   AliDetector::MakeBranch (option, file);
271   const char *cD = strstr(option,"D");
272   const char *cR = strstr(option,"R");
273   const char *cS = strstr(option,"S");
274   
275   if (cS){
276
277     MakeBranchInTree(gAlice->TreeS(), 
278                      branchname,&fSDigits, 
279                      kBufferSize, file);
280   }
281   if (cD){
282
283     MakeBranchInTree(gAlice->TreeD(), 
284                      branchname,&fDigits, 
285                      kBufferSize, file);
286   }
287   if (cR){
288     MakeBranchInTree(gAlice->TreeR(), 
289                      branchname,&fReconParticles,
290                      kBufferSize, file);
291   }
292   
293 }
294
295 //_____________________________________________________________________________
296 void AliFMD::SetTreeAddress ()
297 {
298   // Set branch address for the Hits and Digits Tree.
299   char branchname[30];
300   AliDetector::SetTreeAddress ();
301
302   TBranch *branch;
303   TTree *treeD = gAlice->TreeD ();
304
305
306   if (treeD)
307     {
308       if (fDigits)
309         {
310           branch = treeD->GetBranch (branchname);
311           if (branch)
312             branch->SetAddress (&fDigits);
313         }
314
315     }
316   if (fSDigits)
317     fSDigits->Clear ();
318
319   if (gAlice->TreeS () && fSDigits)
320     {
321       branch = gAlice->TreeS ()->GetBranch ("FMD");
322       if (branch)
323         branch->SetAddress (&fSDigits);
324     }
325
326   if(fReconParticles)
327     fReconParticles->Clear();
328   if (gAlice->TreeR()) 
329     {
330       branch = gAlice->TreeR()->GetBranch("FMD"); 
331       if (branch) branch->SetAddress(&fReconParticles) ;
332     } 
333
334 }
335
336 //---------------------------------------------------------------------
337
338 void AliFMD::SDigits2Digits() 
339 {
340   cout<<"AliFMD::SDigits2Digits"<<endl; 
341     if (fMerger) {
342       cout<<"AliFMD::SDigits2Digits fMerger"<<fMerger<<endl;
343       fMerger->Init();
344       cout<<"AliFMD::SDigits2Digits Init"<<endl; 
345   
346       fMerger->Digitise();
347     }
348
349 }
350 //---------------------------------------------------------------------
351 void   AliFMD::SetMerger(AliFMDMerger* merger)
352 {
353 // Set pointer to merger
354     fMerger = merger;
355 }
356
357 AliFMDMerger*  AliFMD::Merger()
358 {
359 // Return pointer to merger
360     return fMerger;
361 }
362
363 //---------------------------------------------------------------------
364
365
366
367 void
368 AliFMD::Eta2Radius (Float_t eta, Float_t zDisk, Float_t * radius)
369 {
370   Float_t expEta = TMath::Exp (-eta);
371   Float_t theta = TMath::ATan (expEta);
372   theta = 2. * theta;
373   Float_t rad = zDisk * (TMath::Tan (theta));
374   *radius = rad;
375
376   if (fDebug)
377     printf ("%s: eta %f radius %f\n", ClassName (), eta, rad);
378 }
379
380 //---------------------------------------------------------------------
381
382 void AliFMD::Hits2SDigits ()
383 {
384
385   AliFMD *FMD = (AliFMD *) gAlice->GetDetector ("FMD");
386
387   if (fNevents == 0)
388     fNevents = (Int_t) gAlice->TreeE ()->GetEntries ();
389
390   for (Int_t ievent = 0; ievent < fNevents; ievent++)
391     {
392       gAlice->GetEvent (ievent);
393       if (gAlice->TreeH () == 0)
394         return;
395       if (gAlice->TreeS () == 0)
396         gAlice->MakeTree ("S");
397
398
399
400       Int_t nSdigits = 0;
401       
402             //Make branches
403       char branchname[20];
404        sprintf (branchname, "%s", FMD->GetName ());
405       //Make branch for digits
406       FMD->MakeBranch ("S");
407           
408     
409        //Now made SDigits from hits, for PHOS it is the same
410       Int_t volume, sector, ring, charge;
411       Float_t e;
412       Float_t de[10][20][150];
413       Int_t ivol, isec, iring;
414       Int_t hit, nbytes;
415       TParticle *particle;
416       AliFMDhit *fmdHit;
417       TClonesArray *FMDhits = FMD->Hits ();
418
419       // Event ------------------------- LOOP  
420
421       for (ivol = 1; ivol <= 5; ivol++)
422         for (isec = 1; isec <= 16; isec++)
423           for (iring = 1; iring <= 128; iring++)
424             de[ivol][isec][iring] = 0;
425
426       if (FMD)
427         {
428           FMDhits = FMD->Hits ();
429           TTree *TH = gAlice->TreeH ();
430           Stat_t ntracks = TH->GetEntries ();
431           for (Int_t track = 0; track < ntracks; track++)
432             {
433               gAlice->ResetHits ();
434               nbytes += TH->GetEvent (track);
435               particle = gAlice->Particle (track);
436               Int_t nhits = FMDhits->GetEntriesFast ();
437
438               for (hit = 0; hit < nhits; hit++)
439                 {
440                   fmdHit = (AliFMDhit *) FMDhits->UncheckedAt (hit);
441
442                   volume = fmdHit->Volume ();
443                   sector = fmdHit->NumberOfSector ();
444                   ring = fmdHit->NumberOfRing ();
445                   e = fmdHit->Edep ();
446                   de[volume][sector][ring] = de[volume][sector][ring] + e;
447                 }               //hit loop
448             }                   //track loop
449         }                       //if FMD
450
451
452       Int_t digit[5];
453       Float_t I = 1.664 * 0.04 * 2.33 / 22400;  // = 0.69e-6;
454       for (ivol = 1; ivol <= 5; ivol++)
455         {
456           for (isec = 1; isec <= 16; isec++)
457             {
458               for (iring = 1; iring <= 128; iring++)
459                 {
460                       digit[0] = ivol;
461                       digit[1] = isec;
462                       digit[2] = iring;
463                       charge = Int_t (de[ivol][isec][iring] / I);
464                       digit[3] = charge;
465                       //                      if (charge!=0) cout<<" charge "<<charge<<endl;
466                       //dinamic diapason from MIP(0.155MeV) to 30MIP(4.65MeV)
467                       //1024 ADC channels 
468                       Float_t channelWidth = (22400 * 30) / 1024;
469                       digit[4] = Int_t (digit[3] / channelWidth);
470
471                       new ((*fSDigits)[nSdigits++]) AliFMDdigit (digit);
472
473                 }               // iring loop
474             }                   //sector loop
475         }                       // volume loop
476       
477       gAlice->TreeS ()->Fill ();
478       gAlice->TreeS ()->Print ();
479
480     }                           //event loop
481
482 }
483 //-----------------------------------------------------------------------
484
485 void AliFMD::Digits2Reco()
486 {
487 #ifdef DEBUG
488   cout<<"ALiFMD::Digits2Reco> start...";
489 #endif
490   char * fileReconParticles=0;
491   char * fileHeader=0;
492   AliFMDReconstruction * reconstruction =
493     new AliFMDReconstruction(fileHeader,fileReconParticles) ;
494   fReconParticles=new TClonesArray("AliFMDReconstParticles",1000);
495   reconstruction->Exec(fReconParticles,"");
496   delete  reconstruction;
497 }
498