Transition to NewIO
[u/mrichter/AliRoot.git] / TPC / AliTPCDigitizer.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2000, 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 /* $Id$ */
17
18 #include <stdlib.h>
19 #include <TTree.h> 
20 #include <TObjArray.h>
21 #include <TFile.h>
22 #include <TDirectory.h>
23 #include <Riostream.h>
24
25 #include "AliTPCDigitizer.h"
26
27 #include "AliTPC.h"
28 #include "AliTPCParam.h"
29 #include "AliTPCParamSR.h" 
30 #include "AliRun.h"
31 #include "AliPDG.h"
32 #include "AliRunDigitizer.h"
33 #include "AliSimDigits.h"
34
35
36 ClassImp(AliTPCDigitizer)
37
38 //___________________________________________
39 AliTPCDigitizer::AliTPCDigitizer() :AliDigitizer()
40 {
41 // Default ctor - don't use it
42   fDebug =0;
43 }
44
45 //___________________________________________
46 AliTPCDigitizer::AliTPCDigitizer(AliRunDigitizer* manager) 
47     :AliDigitizer(manager)
48 {
49 // ctor which should be used
50   fDebug =0;
51   if (GetDebug()>2) 
52     cerr<<"AliTPCDigitizer::AliTPCDigitizer"
53        <<"(AliRunDigitizer* manager) was processed"<<endl;
54 }
55
56 //------------------------------------------------------------------------
57 AliTPCDigitizer::~AliTPCDigitizer()
58 {
59 // Destructor
60 }
61
62
63
64 //------------------------------------------------------------------------
65 Bool_t AliTPCDigitizer::Init()
66 {
67 // Initialization 
68     
69  return kTRUE;
70 }
71
72
73 //------------------------------------------------------------------------
74 void AliTPCDigitizer::Exec(Option_t* option)
75 {
76   ExecFast(option);  
77 }
78 //------------------------------------------------------------------------
79 void AliTPCDigitizer::ExecFast(Option_t* option)
80 {
81   
82   // merge input tree's with summable digits
83   //output stored in TreeTPCD
84   char s[100]; 
85   char ss[100];
86   TString optionString = option;
87   if (optionString.Data() == "deb") {
88     cout<<"AliTPCDigitizer::Exec: called with option deb "<<endl;
89     fDebug = 3;
90   }
91   //get detector and geometry
92
93
94   AliRunLoader *rl, *orl;
95   AliLoader *gime, *ogime;
96   
97   if (gAlice == 0x0)
98    {
99      Warning("ExecFast","gAlice is NULL. Loading from input 0");
100      rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
101      if (rl == 0x0)
102       {
103         Error("ExecFast","Can not find Run Loader for input 0. Can not proceed.");
104         return;
105       }
106      rl->LoadgAlice();
107      rl->GetAliRun();
108    }
109   AliTPC *pTPC  = (AliTPC *) gAlice->GetModule("TPC");
110   AliTPCParam * param = pTPC->GetParam();
111   
112   sprintf(s,param->GetTitle());
113   sprintf(ss,"75x40_100x60");
114   if(strcmp(s,ss)==0){
115     printf("2 pad-length geom hits with 3 pad-lenght geom digits...\n");
116     delete param;
117     param=new AliTPCParamSR();
118   }
119   else{
120    sprintf(ss,"75x40_100x60_150x60");
121    if(strcmp(s,ss)!=0) {
122      printf("No TPC parameters found...\n");
123      exit(2); 
124    }
125   }
126   
127   pTPC->GenerNoise(500000); //create teble with noise
128   //
129   Int_t nInputs = fManager->GetNinputs();
130   Int_t * masks = new Int_t[nInputs];
131   for (Int_t i=0; i<nInputs;i++)
132     masks[i]= fManager->GetMask(i);
133   Short_t **pdig= new Short_t*[nInputs];   //pointers to the expanded digits array
134   Int_t **ptr=  new Int_t*[nInputs];       //pointers to teh expanded tracks array
135
136   
137   //create digits array for given sectors
138   // make indexes
139   Stat_t nentries = 0;//number of entries in TreeS
140   
141   AliSimDigits ** digarr = new AliSimDigits*[nInputs]; 
142   for (Int_t i1=0;i1<nInputs; i1++)
143     {
144       digarr[i1]=0;
145      //    intree[i1]
146       rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(i1));
147       gime = rl->GetLoader("TPCLoader");
148       gime->LoadSDigits("read");
149       TTree * treear =  gime->TreeS();
150      
151       if (!treear) 
152        {
153         cerr<<"AliTPCDigitizer: Input tree with SDigits not found in"
154             <<" input "<< i1<<endl;
155         return;
156        }
157       if (i1 == 0) nentries=treear->GetEntries();
158     
159       if (treear->GetIndex()==0)  
160         treear->BuildIndex("fSegmentID","fSegmentID");
161       treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
162     }
163
164
165   //create branch's in TPC treeD
166   AliSimDigits * digrow = new AliSimDigits;
167
168   orl = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
169   ogime = orl->GetLoader("TPCLoader");
170   
171   TTree * tree  = ogime->TreeD();
172   if (tree == 0x0)
173    {
174      ogime->MakeTree("D");
175      tree  = ogime->TreeD();
176    }
177   tree->Branch("Segment","AliSimDigits",&digrow);
178   //
179
180   param->SetZeroSup(2);
181
182   Int_t zerosup = param->GetZeroSup(); 
183   //
184   //Loop over segments of the TPC
185     
186   for (Int_t n=0; n<nentries; n++) 
187    {
188     rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
189     gime = rl->GetLoader("TPCLoader");
190     gime->TreeS()->GetEvent(n);
191     digarr[0]->ExpandBuffer();
192     digarr[0]->ExpandTrackBuffer();
193            
194     for (Int_t i=1;i<nInputs; i++) 
195      { 
196
197       rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(i));
198       gime = rl->GetLoader("TPCLoader");
199       
200       gime->TreeS()->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());
201       if ((digarr[0]->GetID()-digarr[i]->GetID())>0) 
202        printf("problem - not corresponding segment in background event\n");
203       
204       digarr[i]->ExpandBuffer();
205       digarr[i]->ExpandTrackBuffer();
206       
207      }   
208     Int_t sec, row;
209     if (!param->AdjustSectorRow(digarr[0]->GetID(),sec,row)) 
210      {
211       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr[0]->GetID()<<endl;
212       continue;
213      }
214
215     digrow->SetID(digarr[0]->GetID());
216
217     Int_t nrows = digarr[0]->GetNRows();
218     Int_t ncols = digarr[0]->GetNCols();
219     digrow->Allocate(nrows,ncols);
220     digrow->AllocateTrack(3);
221
222     Float_t q=0;
223     Int_t label[1000]; //stack for 300 events 
224     Int_t labptr = 0;
225
226     Int_t nElems = nrows*ncols;     
227  
228     for (Int_t i=0;i<nInputs; i++)
229      { 
230        pdig[i] = digarr[i]->GetDigits();
231        ptr[i]  = digarr[i]->GetTracks();
232      }
233      
234     Short_t *pdig1= digrow->GetDigits();
235     Int_t   *ptr1= digrow->GetTracks() ;
236
237     
238
239     for (Int_t elem=0;elem<nElems; elem++)
240      {    
241
242        q=0;
243        labptr=0;
244        // looop over digits 
245         for (Int_t i=0;i<nInputs; i++)
246          { 
247           //          q  += digarr[i]->GetDigitFast(rows,col);
248             q  += *(pdig[i]);
249          
250            for (Int_t tr=0;tr<3;tr++) 
251             {
252              //             Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
253              Int_t lab = ptr[i][tr*nElems];
254              if ( (lab > 1) && *(pdig[i])>zerosup) 
255               {
256                 label[labptr]=lab+masks[i];
257                 labptr++;
258               }          
259             }
260            pdig[i]++;
261            ptr[i]++;
262          }
263         q/=16.;  //conversion factor
264         //       Float_t noise  = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());  
265         Float_t noise  = pTPC->GetNoise();
266         q+=noise;
267         q=TMath::Nint(q);
268         if (q > zerosup)
269          { 
270           if(q > param->GetADCSat()) q = (Short_t)(param->GetADCSat());
271           //digrow->SetDigitFast((Short_t)q,rows,col);  
272           *pdig1 =Short_t(q);
273           for (Int_t tr=0;tr<3;tr++)
274            {
275             if (tr<labptr) 
276              ptr1[tr*nElems] = label[tr];
277            }
278           }
279         pdig1++;
280         ptr1++;
281      }
282     
283     digrow->CompresBuffer(1,zerosup);
284     digrow->CompresTrackBuffer(1);
285     tree->Fill();
286     if (fDebug>0) cerr<<sec<<"\t"<<row<<"\n";  
287    } //for (Int_t n=0; n<nentries; n++) 
288   
289
290   orl = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
291   ogime = orl->GetLoader("TPCLoader");
292   ogime->WriteDigits("OVERWRITE");
293   
294   //fManager->GetTreeDTPC()->Write(0,TObject::kOverwrite);
295   
296   delete digrow;     
297   for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
298   delete []masks;
299   delete digarr;  
300 }
301
302
303
304 //------------------------------------------------------------------------
305 void AliTPCDigitizer::ExecSave(Option_t* option)
306 {
307   
308   // merge input tree's with summable digits
309   //output stored in TreeTPCD
310
311   TString optionString = option;
312   if (optionString.Data() == "deb") {
313     cout<<"AliTPCDigitizer::Exec: called with option deb "<<endl;
314     fDebug = 3;
315   }
316   //get detector and geometry 
317   AliRunLoader *rl, *orl;
318   AliLoader *gime, *ogime;
319
320   
321   orl = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
322   ogime = orl->GetLoader("TPCLoader");
323   
324   rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
325   gime = rl->GetLoader("TPCLoader");
326   
327   rl->LoadgAlice();
328   AliRun* alirun = rl->GetAliRun();
329   
330   AliTPC *pTPC  = (AliTPC *) alirun->GetModule("TPC");
331   AliTPCParam * param = pTPC->GetParam();
332   pTPC->GenerNoise(500000); //create teble with noise
333   printf("noise %f \n",  param->GetNoise()*param->GetNoiseNormFac());
334   //
335   Int_t nInputs = fManager->GetNinputs();
336   Int_t * masks = new Int_t[nInputs];
337   for (Int_t i=0; i<nInputs;i++)
338     masks[i]= fManager->GetMask(i);
339
340   AliSimDigits ** digarr = new AliSimDigits*[nInputs]; 
341   for (Int_t i1=0;i1<nInputs; i1++)
342    {
343     digarr[i1]=0;
344     //    intree[i1]
345     rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(i1));
346     gime = rl->GetLoader("TPCLoader");
347
348     TTree * treear =  gime->TreeS();
349     TBranch * br = treear->GetBranch("fSegmentID");
350     if (br) br->GetFile()->cd();
351     if (!treear) {      
352       cerr<<" TPC -  not existing input = \n"<<i1<<" ";      
353     } 
354     treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
355   }
356   
357   rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
358   gime = rl->GetLoader("TPCLoader");
359   Stat_t nentries = gime->TreeS()->GetEntries();
360   
361
362   //create branch's in TPC treeD
363   AliSimDigits * digrow = new AliSimDigits;
364   TTree * tree  = ogime->TreeD();
365
366   tree->Branch("Segment","AliSimDigits",&digrow);
367   param->SetZeroSup(2);
368
369   Int_t zerosup = param->GetZeroSup();
370   //Loop over segments of the TPC
371     
372   for (Int_t n=0; n<nentries; n++) {
373     rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
374     gime = rl->GetLoader("TPCLoader");
375     gime->TreeS()->GetEvent(n);
376
377     digarr[0]->ExpandBuffer();
378     digarr[0]->ExpandTrackBuffer();
379            
380     for (Int_t i=1;i<nInputs; i++){ 
381 //      fManager->GetInputTreeTPCS(i)->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());      
382       rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(i));
383       gime = rl->GetLoader("TPCLoader");
384       gime->TreeS()->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());  
385       digarr[i]->ExpandBuffer();
386       digarr[i]->ExpandTrackBuffer();
387       if ((digarr[0]->GetID()-digarr[i]->GetID())>0) 
388        printf("problem\n");
389     
390     }   
391     
392     Int_t sec, row;
393     if (!param->AdjustSectorRow(digarr[0]->GetID(),sec,row)) {
394       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr[0]->GetID()<<endl;
395       continue;
396     }
397
398     digrow->SetID(digarr[0]->GetID());
399
400     Int_t nrows = digarr[0]->GetNRows();
401     Int_t ncols = digarr[0]->GetNCols();
402     digrow->Allocate(nrows,ncols);
403     digrow->AllocateTrack(3);
404
405     Float_t q=0;
406     Int_t label[1000]; //stack for 300 events 
407     Int_t labptr = 0;
408
409     
410
411     for (Int_t rows=0;rows<nrows; rows++){
412       for (Int_t col=0;col<ncols; col++){
413     
414        q=0;
415        labptr=0;
416        // looop over digits 
417         for (Int_t i=0;i<nInputs; i++){ 
418          q  += digarr[i]->GetDigitFast(rows,col);
419           //q  += *(pdig[i]);
420          
421           for (Int_t tr=0;tr<3;tr++) {
422            Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
423            //Int_t lab = ptr[i][tr*nElems];
424             if ( (lab > 1) ) {
425               label[labptr]=lab+masks[i];
426               labptr++;
427             }          
428           }
429          // pdig[i]++;
430          //ptr[i]++;
431          
432         }
433        q/=16.;  //conversion factor
434        //       Float_t noise  = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());  
435        Float_t noise  = pTPC->GetNoise();
436        q+=noise;
437         q=TMath::Nint(q);
438         if (q > zerosup){ 
439          
440          if(q > param->GetADCSat()) q = (Short_t)(param->GetADCSat());
441          digrow->SetDigitFast((Short_t)q,rows,col);  
442          // *pdig1 =Short_t(q);
443          for (Int_t tr=0;tr<3;tr++){
444            if (tr<labptr) 
445              ((AliSimDigits*)digrow)->SetTrackIDFast(label[tr],rows,col,tr);
446            //ptr1[tr*nElems] = label[tr];
447            //else
448              //           ((AliSimDigits*)digrow)->SetTrackIDFast(-1,rows,col,tr);          
449            //  ptr1[tr*nElems] = 1;
450          }
451        }
452        //pdig1++;
453        //ptr1++;
454     }
455     }
456     
457     digrow->CompresBuffer(1,zerosup);
458     digrow->CompresTrackBuffer(1);
459     tree->Fill();
460     if (fDebug>0) cerr<<sec<<"\t"<<row<<"\n";  
461   } 
462 //  printf("end TPC merging - end -Tree %s\t%p\n",fManager->GetInputTreeH(0)->GetName(),fManager->GetInputTreeH(0)->GetListOfBranches()->At(3));
463   //fManager->GetTreeDTPC()->Write(0,TObject::kOverwrite);
464     ogime->WriteDigits("OVERWRITE");
465
466     for (Int_t i=1;i<nInputs; i++) 
467      { 
468       rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(i));
469       gime = rl->GetLoader("TPCLoader");
470       gime->UnloadSDigits();
471      }
472     ogime->UnloadDigits();
473     
474   delete digrow;     
475   for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
476   delete []masks;
477   delete digarr;  
478 }