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