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