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