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