Just something left from v2...now clean
[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   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); 
189   //
190   //Loop over segments of the TPC
191     
192   for (Int_t segmentID=0; segmentID<param->GetNRowsTotal(); segmentID++) 
193    {
194     Int_t sec, row;
195     if (!param->AdjustSectorRow(segmentID,sec,row)) 
196      {
197       cerr<<"AliTPC warning: invalid segment ID ! "<<segmentID<<endl;
198       continue;
199      }
200     AliTPCCalROC * gainROC = gainTPC->GetCalROC(sec);  // pad gains per given sector
201     AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sec);  // noise per given sector
202     digrow->SetID(segmentID);
203
204     Int_t nrows = 0;
205     Int_t ncols = 0;
206
207     Bool_t digitize = kFALSE;
208     for (Int_t i=0;i<nInputs; i++) 
209      { 
210
211       rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(i));
212       gime = rl->GetLoader("TPCLoader");
213       
214       if (gime->TreeS()->GetEntryWithIndex(segmentID,segmentID) >= 0) {
215         digarr[i]->ExpandBuffer();
216         digarr[i]->ExpandTrackBuffer();
217         nrows = digarr[i]->GetNRows();
218         ncols = digarr[i]->GetNCols();
219         active[i] = kTRUE;
220         if (!fRegionOfInterest || (i == 0)) digitize = kTRUE;
221       } else {
222         active[i] = kFALSE;
223       }
224       if (fRegionOfInterest && !digitize) break;
225      }   
226     if (!digitize) continue;
227
228     digrow->Allocate(nrows,ncols);
229     digrow->AllocateTrack(3);
230
231     Float_t q=0;
232     Int_t label[1000]; //stack for 300 events 
233     Int_t labptr = 0;
234
235     Int_t nElems = nrows*ncols;     
236  
237     for (Int_t i=0;i<nInputs; i++)
238      if (active[i]) { 
239        pdig[i] = digarr[i]->GetDigits();
240        ptr[i]  = digarr[i]->GetTracks();
241       }
242      
243     Short_t *pdig1= digrow->GetDigits();
244     Int_t   *ptr1= digrow->GetTracks() ;
245
246     
247
248     for (Int_t elem=0;elem<nElems; elem++)
249      {    
250
251        q=0;
252        labptr=0;
253        // looop over digits 
254         for (Int_t i=0;i<nInputs; i++) if (active[i]) 
255          { 
256           //          q  += digarr[i]->GetDigitFast(rows,col);
257             q  += *(pdig[i]);
258          
259            for (Int_t tr=0;tr<3;tr++) 
260             {
261              //             Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
262              Int_t lab = ptr[i][tr*nElems];
263              if ( (lab > 1) && *(pdig[i])>zerosup) 
264               {
265                 label[labptr]=lab+masks[i];
266                 labptr++;
267               }          
268             }
269            pdig[i]++;
270            ptr[i]++;
271          }
272         q/=16.;  //conversion factor
273         Float_t gain = gainROC->GetValue(row,elem/nrows);  // get gain for given - pad-row pad
274         if (gain<0.5){
275           printf("problem\n");
276         }
277         q*= gain;
278         Float_t noisePad = noiseROC->GetValue(row,elem/nrows);
279         //       Float_t noise  = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());  
280         Float_t noise  = pTPC->GetNoise();
281         q+=noise*noisePad;
282         q=TMath::Nint(q);
283         if (q > zerosup)
284          { 
285           if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
286           //digrow->SetDigitFast((Short_t)q,rows,col);  
287           *pdig1 =Short_t(q);
288           for (Int_t tr=0;tr<3;tr++)
289            {
290             if (tr<labptr) 
291              ptr1[tr*nElems] = label[tr];
292            }
293           }
294         pdig1++;
295         ptr1++;
296      }
297     
298     digrow->CompresBuffer(1,zerosup);
299     digrow->CompresTrackBuffer(1);
300     tree->Fill();
301     if (fDebug>0) cerr<<sec<<"\t"<<row<<"\n";  
302    } //for (Int_t n=0; n<param->GetNRowsTotal(); n++) 
303   
304
305   orl = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
306   ogime = orl->GetLoader("TPCLoader");
307   ogime->WriteDigits("OVERWRITE");
308   
309   //fManager->GetTreeDTPC()->Write(0,TObject::kOverwrite);
310   
311   delete digrow;     
312   for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
313   delete []masks;
314   delete []pdig;
315   delete []ptr;
316   delete []active;
317   delete []digarr;  
318 }
319
320
321
322 //------------------------------------------------------------------------
323 void AliTPCDigitizer::ExecSave(Option_t* option)
324 {
325   
326   // merge input tree's with summable digits
327   //output stored in TreeTPCD
328
329   TString optionString = option;
330   if (optionString.Data() == "deb") {
331     cout<<"AliTPCDigitizer::Exec: called with option deb "<<endl;
332     fDebug = 3;
333   }
334   //get detector and geometry 
335   AliRunLoader *rl, *orl;
336   AliLoader *gime, *ogime;
337
338   
339   orl = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
340   ogime = orl->GetLoader("TPCLoader");
341   
342   rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
343   gime = rl->GetLoader("TPCLoader");
344   
345   rl->LoadgAlice();
346   AliRun* alirun = rl->GetAliRun();
347   
348   AliTPC *pTPC  = (AliTPC *) alirun->GetModule("TPC");
349   AliTPCParam * param = pTPC->GetParam();
350   pTPC->GenerNoise(500000); //create teble with noise
351   printf("noise %f \n",  param->GetNoise()*param->GetNoiseNormFac());
352   //
353   Int_t nInputs = fManager->GetNinputs();
354   Int_t * masks = new Int_t[nInputs];
355   for (Int_t i=0; i<nInputs;i++)
356     masks[i]= fManager->GetMask(i);
357
358   AliSimDigits ** digarr = new AliSimDigits*[nInputs]; 
359   for (Int_t i1=0;i1<nInputs; i1++)
360    {
361     digarr[i1]=0;
362     //    intree[i1]
363     rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(i1));
364     gime = rl->GetLoader("TPCLoader");
365
366     TTree * treear =  gime->TreeS();
367     TBranch * br = treear->GetBranch("fSegmentID");
368     if (br) br->GetFile()->cd();
369     if (!treear) {      
370       cerr<<" TPC -  not existing input = \n"<<i1<<" ";      
371     } 
372     treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
373   }
374   
375   rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
376   gime = rl->GetLoader("TPCLoader");
377   Stat_t nentries = gime->TreeS()->GetEntries();
378   
379
380   //create branch's in TPC treeD
381   AliSimDigits * digrow = new AliSimDigits;
382   TTree * tree  = ogime->TreeD();
383
384   tree->Branch("Segment","AliSimDigits",&digrow);
385   param->SetZeroSup(2);
386
387   Int_t zerosup = param->GetZeroSup();
388   //Loop over segments of the TPC
389     
390   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
391   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
392   for (Int_t n=0; n<nentries; n++) {
393     rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
394     gime = rl->GetLoader("TPCLoader");
395     gime->TreeS()->GetEvent(n);
396
397     digarr[0]->ExpandBuffer();
398     digarr[0]->ExpandTrackBuffer();
399
400
401     for (Int_t i=1;i<nInputs; i++){ 
402 //      fManager->GetInputTreeTPCS(i)->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());      
403       rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(i));
404       gime = rl->GetLoader("TPCLoader");
405       gime->TreeS()->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());  
406       digarr[i]->ExpandBuffer();
407       digarr[i]->ExpandTrackBuffer();
408       if ((digarr[0]->GetID()-digarr[i]->GetID())>0) 
409        printf("problem\n");
410     
411     }   
412     
413     Int_t sec, row;
414     if (!param->AdjustSectorRow(digarr[0]->GetID(),sec,row)) {
415       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr[0]->GetID()<<endl;
416       continue;
417     }
418
419     AliTPCCalROC * gainROC = gainTPC->GetCalROC(sec);  // pad gains per given sector
420     AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sec);  // noise per given sector
421     digrow->SetID(digarr[0]->GetID());
422
423     Int_t nrows = digarr[0]->GetNRows();
424     Int_t ncols = digarr[0]->GetNCols();
425     digrow->Allocate(nrows,ncols);
426     digrow->AllocateTrack(3);
427
428     Float_t q=0;
429     Int_t label[1000]; //stack for 300 events 
430     Int_t labptr = 0;
431
432     
433
434     for (Int_t rows=0;rows<nrows; rows++){
435       for (Int_t col=0;col<ncols; col++){
436     
437        q=0;
438        labptr=0;
439        // looop over digits 
440         for (Int_t i=0;i<nInputs; i++){ 
441          q  += digarr[i]->GetDigitFast(rows,col);
442           //q  += *(pdig[i]);
443          
444           for (Int_t tr=0;tr<3;tr++) {
445            Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
446            //Int_t lab = ptr[i][tr*nElems];
447             if ( (lab > 1) ) {
448               label[labptr]=lab+masks[i];
449               labptr++;
450             }          
451           }
452          // pdig[i]++;
453          //ptr[i]++;
454          
455         }
456        q/=16.;  //conversion factor
457        //       Float_t noise  = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());  
458        Float_t gain = gainROC->GetValue(row,col);
459        q*= gain;
460        Float_t noisePad = noiseROC->GetValue(row, col);
461
462        Float_t noise  = pTPC->GetNoise();
463        q+=noise*noisePad;
464         q=TMath::Nint(q);
465         if (q > zerosup){ 
466          
467          if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
468          digrow->SetDigitFast((Short_t)q,rows,col);  
469          // *pdig1 =Short_t(q);
470          for (Int_t tr=0;tr<3;tr++){
471            if (tr<labptr) 
472              ((AliSimDigits*)digrow)->SetTrackIDFast(label[tr],rows,col,tr);
473            //ptr1[tr*nElems] = label[tr];
474            //else
475              //           ((AliSimDigits*)digrow)->SetTrackIDFast(-1,rows,col,tr);          
476            //  ptr1[tr*nElems] = 1;
477          }
478        }
479        //pdig1++;
480        //ptr1++;
481     }
482     }
483     
484     digrow->CompresBuffer(1,zerosup);
485     digrow->CompresTrackBuffer(1);
486     tree->Fill();
487     if (fDebug>0) cerr<<sec<<"\t"<<row<<"\n";  
488   } 
489 //  printf("end TPC merging - end -Tree %s\t%p\n",fManager->GetInputTreeH(0)->GetName(),fManager->GetInputTreeH(0)->GetListOfBranches()->At(3));
490   //fManager->GetTreeDTPC()->Write(0,TObject::kOverwrite);
491     ogime->WriteDigits("OVERWRITE");
492
493     for (Int_t i=1;i<nInputs; i++) 
494      { 
495       rl = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(i));
496       gime = rl->GetLoader("TPCLoader");
497       gime->UnloadSDigits();
498      }
499     ogime->UnloadDigits();
500     
501   delete digrow;     
502   for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
503   delete []masks;
504   delete digarr;  
505 }