690bb0e2f166af9c3a349dc011232688925a4a49
[u/mrichter/AliRoot.git] / TPC / TPCsim / 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 "AliDigitizationInput.h"
52 #include "AliSimDigits.h"
53 #include "AliLog.h"
54
55 #include "AliTPCcalibDB.h"
56 #include "AliTPCCalPad.h"
57 #include "AliTPCCalROC.h"
58 #include "TTreeStream.h"
59 #include "AliTPCReconstructor.h"
60 #include <TGraphErrors.h>
61
62 using std::cout;
63 using std::cerr;
64 using std::endl;
65 ClassImp(AliTPCDigitizer)
66
67 //___________________________________________
68   AliTPCDigitizer::AliTPCDigitizer() :AliDigitizer(),fDebug(0), fDebugStreamer(0)
69 {
70   //
71 // Default ctor - don't use it
72 //
73   
74 }
75
76 //___________________________________________
77 AliTPCDigitizer::AliTPCDigitizer(AliDigitizationInput* digInput) 
78   :AliDigitizer(digInput),fDebug(0), fDebugStreamer(0)
79 {
80   //
81 // ctor which should be used
82 //  
83   AliDebug(2,"(AliDigitizationInput* digInput) was processed");
84   if (AliTPCReconstructor::StreamLevel()>0)  fDebugStreamer = new TTreeSRedirector("TPCDigitDebug.root","recreate");
85
86 }
87
88 //------------------------------------------------------------------------
89 AliTPCDigitizer::~AliTPCDigitizer()
90 {
91 // Destructor
92   if (fDebugStreamer) delete fDebugStreamer;
93 }
94
95
96
97 //------------------------------------------------------------------------
98 Bool_t AliTPCDigitizer::Init()
99 {
100 // Initialization 
101     
102  return kTRUE;
103 }
104
105
106 //------------------------------------------------------------------------
107 void AliTPCDigitizer::Digitize(Option_t* option)
108 {
109   //DigitizeFast(option);  
110   DigitizeWithTailAndCrossTalk(option);
111   
112 }
113 //------------------------------------------------------------------------
114 void AliTPCDigitizer::DigitizeFast(Option_t* option)
115 {
116   
117   // merge input tree's with summable digits
118   //output stored in TreeTPCD
119   char s[100]; 
120   char ss[100];
121   TString optionString = option;
122   if (!strcmp(optionString.Data(),"deb")) {
123     cout<<"AliTPCDigitizer:::DigitizeFast called with option deb "<<endl;
124     fDebug = 3;
125   }
126   //get detector and geometry
127
128
129   AliRunLoader *rl, *orl;
130   AliLoader *gime, *ogime;
131   
132   if (gAlice == 0x0)
133    {
134      Warning("DigitizeFast","gAlice is NULL. Loading from input 0");
135      rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
136      if (rl == 0x0)
137       {
138         Error("DigitizeFast","Can not find Run Loader for input 0. Can not proceed.");
139         return;
140       }
141      rl->LoadgAlice();
142      rl->GetAliRun();
143    }
144   AliTPC *pTPC  = (AliTPC *) gAlice->GetModule("TPC");
145   AliTPCParam * param = pTPC->GetParam();
146   
147   //sprintf(s,param->GetTitle());
148   snprintf(s,100,"%s",param->GetTitle());
149   //sprintf(ss,"75x40_100x60");
150   snprintf(ss,100,"75x40_100x60");
151   if(strcmp(s,ss)==0){
152     printf("2 pad-length geom hits with 3 pad-lenght geom digits...\n");
153     delete param;
154     param=new AliTPCParamSR();
155   }
156   else{
157     //sprintf(ss,"75x40_100x60_150x60");
158     snprintf(ss,100,"75x40_100x60_150x60");
159    if(strcmp(s,ss)!=0) {
160      printf("No TPC parameters found...\n");
161      exit(2); 
162    }
163   }
164   
165   pTPC->GenerNoise(500000); //create table with noise
166   //
167   Int_t nInputs = fDigInput->GetNinputs();
168   Int_t * masks = new Int_t[nInputs];
169   for (Int_t i=0; i<nInputs;i++)
170     masks[i]= fDigInput->GetMask(i);
171   Short_t **pdig= new Short_t*[nInputs];   //pointers to the expanded digits array
172   Int_t **ptr=  new Int_t*[nInputs];       //pointers to the expanded tracks array
173   Bool_t *active=  new Bool_t[nInputs];    //flag for active input segments
174   Char_t phname[100];
175   
176   //create digits array for given sectors
177   // make indexes
178   //
179   //create branch's in TPC treeD
180   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
181   ogime = orl->GetLoader("TPCLoader");
182   TTree * tree  = ogime->TreeD();
183   AliSimDigits * digrow = new AliSimDigits;  
184
185   if (tree == 0x0)
186    {
187      ogime->MakeTree("D");
188      tree  = ogime->TreeD();
189    }
190   tree->Branch("Segment","AliSimDigits",&digrow);
191   //  
192   AliSimDigits ** digarr = new AliSimDigits*[nInputs]; 
193   for (Int_t i1=0;i1<nInputs; i1++)
194     {
195       digarr[i1]=0;
196      //    intree[i1]
197       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
198       gime = rl->GetLoader("TPCLoader");
199       gime->LoadSDigits("read");
200       TTree * treear =  gime->TreeS();
201      
202       if (!treear) 
203        {
204         cerr<<"AliTPCDigitizer: Input tree with SDigits not found in"
205             <<" input "<< i1<<endl;
206         for (Int_t i2=0;i2<i1+1; i2++){ 
207           
208           if(digarr[i2])  delete digarr[i2];
209         }
210         delete [] digarr;
211         delete [] active;
212         delete []masks;
213         delete []pdig;
214         delete []ptr;
215         return;
216        }
217
218       //sprintf(phname,"lhcphase%d",i1);
219       snprintf(phname,100,"lhcphase%d",i1);
220       TParameter<float> *ph = (TParameter<float>*)treear->GetUserInfo()
221                                ->FindObject("lhcphase0");
222       if(!ph){
223         cerr<<"AliTPCDigitizer: LHC phase  not found in"
224             <<" input "<< i1<<endl;
225         for (Int_t i2=0;i2<i1+1; i2++){ 
226           if(digarr[i2])  delete digarr[i2];
227         }
228         delete [] digarr;
229         delete [] active;
230         delete []masks;
231         delete []pdig;
232         delete []ptr;
233         return;
234       }
235       tree->GetUserInfo()->Add(new TParameter<float>(phname,ph->GetVal()));
236               //
237       if (treear->GetIndex()==0)  
238         treear->BuildIndex("fSegmentID","fSegmentID");
239       treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
240     }
241
242
243
244
245   //
246
247   param->SetZeroSup(2);
248
249   Int_t zerosup = param->GetZeroSup(); 
250   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor(); 
251   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); 
252   //
253   //Loop over segments of the TPC
254     
255   for (Int_t segmentID=0; segmentID<param->GetNRowsTotal(); segmentID++) 
256    {
257     Int_t sector, padRow;
258     if (!param->AdjustSectorRow(segmentID,sector,padRow)) 
259      {
260       cerr<<"AliTPC warning: invalid segment ID ! "<<segmentID<<endl;
261       continue;
262      }
263     AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector);  // pad gains per given sector
264     AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector);  // noise per given sector
265     digrow->SetID(segmentID);
266
267     Int_t nTimeBins = 0;
268     Int_t nPads = 0;
269
270     Bool_t digitize = kFALSE;
271     for (Int_t i=0;i<nInputs; i++) 
272      { 
273
274       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
275       gime = rl->GetLoader("TPCLoader");
276       
277       if (gime->TreeS()->GetEntryWithIndex(segmentID,segmentID) >= 0) {
278         digarr[i]->ExpandBuffer();
279         digarr[i]->ExpandTrackBuffer();
280         nTimeBins = digarr[i]->GetNRows();
281         nPads = digarr[i]->GetNCols();
282         active[i] = kTRUE;
283         if (!GetRegionOfInterest() || (i == 0)) digitize = kTRUE;
284       } else {
285         active[i] = kFALSE;
286       }
287       if (GetRegionOfInterest() && !digitize) break;
288      }   
289     if (!digitize) continue;
290
291     digrow->Allocate(nTimeBins,nPads);
292     digrow->AllocateTrack(3);
293
294     Float_t q=0;
295     Int_t label[1000]; //stack for 300 events 
296     Int_t labptr = 0;
297
298     Int_t nElems = nTimeBins*nPads;     
299  
300     for (Int_t i=0;i<nInputs; i++)
301      if (active[i]) { 
302        pdig[i] = digarr[i]->GetDigits();
303        ptr[i]  = digarr[i]->GetTracks();
304       }
305      
306     Short_t *pdig1= digrow->GetDigits();
307     Int_t   *ptr1= digrow->GetTracks() ;
308
309     
310
311     for (Int_t elem=0;elem<nElems; elem++)
312      {    
313
314        q=0;
315        labptr=0;
316        // looop over digits 
317         for (Int_t i=0;i<nInputs; i++) if (active[i]) 
318          { 
319           //          q  += digarr[i]->GetDigitFast(rows,col);
320             q  += *(pdig[i]);
321          
322            for (Int_t tr=0;tr<3;tr++) 
323             {
324              //             Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
325              Int_t lab = ptr[i][tr*nElems];
326              if ( (lab > 1) && *(pdig[i])>zerosup) 
327               {
328                 label[labptr]=lab+masks[i];
329                 labptr++;
330               }          
331             }
332            pdig[i]++;
333            ptr[i]++;
334          }
335         q/=16.;  //conversion factor
336         Float_t gain = gainROC->GetValue(padRow,elem/nTimeBins);  // get gain for given - pad-row pad
337         //if (gain<0.5){
338           //printf("problem\n");
339         //}
340         q*= gain;
341         Float_t noisePad = noiseROC->GetValue(padRow,elem/nTimeBins);
342         //       Float_t noise  = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());  
343         Float_t noise  = pTPC->GetNoise();
344         q+=noise*noisePad;
345         q=TMath::Nint(q);
346         if (q > zerosup)
347          { 
348           if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
349           //digrow->SetDigitFast((Short_t)q,rows,col);  
350           *pdig1 =Short_t(q);
351           for (Int_t tr=0;tr<3;tr++)
352            {
353             if (tr<labptr) 
354              ptr1[tr*nElems] = label[tr];
355            }
356           }
357         pdig1++;
358         ptr1++;
359      }
360     //
361     //  glitch filter
362     //
363     digrow->GlitchFilter();
364     //
365     digrow->CompresBuffer(1,zerosup);
366     digrow->CompresTrackBuffer(1);
367     tree->Fill();
368     if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";  
369    } //for (Int_t n=0; n<param->GetNRowsTotal(); n++) 
370   
371
372   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
373   ogime = orl->GetLoader("TPCLoader");
374   ogime->WriteDigits("OVERWRITE");
375   
376   //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
377   
378   delete digrow;     
379   for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
380   delete []masks;
381   delete []pdig;
382   delete []ptr;
383   delete []active;
384   delete []digarr;  
385 }
386
387
388
389 //------------------------------------------------------------------------
390 void AliTPCDigitizer::DigitizeSave(Option_t* option)
391 {
392   //
393   // Merge input tree's with summable digits
394   // Output digits stored in TreeTPCD
395   // 
396   // Not active for long time.
397   // Before adding modification (for ion tail calucation and for the crorsstalk) it should be 
398   //  checked one by one with currenlty used AliTPCDigitizer::DigitizeFast
399   //
400   TString optionString = option;
401   if (!strcmp(optionString.Data(),"deb")) {
402     cout<<"AliTPCDigitizer::Digitize: called with option deb "<<endl;
403     fDebug = 3;
404   }
405   //get detector and geometry 
406   AliRunLoader *rl, *orl;
407   AliLoader *gime, *ogime;
408
409   
410   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
411   ogime = orl->GetLoader("TPCLoader");
412   
413   rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
414   //gime = rl->GetLoader("TPCLoader");
415   rl->GetLoader("TPCLoader");
416   rl->LoadgAlice();
417   AliRun* alirun = rl->GetAliRun();
418   
419   AliTPC *pTPC  = (AliTPC *) alirun->GetModule("TPC");
420   AliTPCParam * param = pTPC->GetParam();
421   pTPC->GenerNoise(500000); //create teble with noise
422   printf("noise %f \n",  param->GetNoise()*param->GetNoiseNormFac());
423   //
424   Int_t nInputs = fDigInput->GetNinputs();
425   // stupid protection...
426   if (nInputs <= 0) return;
427   //
428   Int_t * masks = new Int_t[nInputs];
429   for (Int_t i=0; i<nInputs;i++)
430     masks[i]= fDigInput->GetMask(i);
431
432   AliSimDigits ** digarr = new AliSimDigits*[nInputs]; 
433   for(Int_t ii=0;ii<nInputs;ii++) digarr[ii]=0;
434
435   for (Int_t i1=0;i1<nInputs; i1++)
436    {
437      //digarr[i1]=0;
438     //    intree[i1]
439     rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
440     gime = rl->GetLoader("TPCLoader");
441
442     TTree * treear =  gime->TreeS();
443     //
444     if (!treear) {      
445       cerr<<" TPC -  not existing input = \n"<<i1<<" "; 
446       delete [] masks;  
447       for(Int_t i=0; i<nInputs; i++) delete digarr[i];
448       delete [] digarr;
449       return;   
450     } 
451     //
452     TBranch * br = treear->GetBranch("fSegmentID");
453     if (br) br->GetFile()->cd();
454     treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
455    }
456   
457   rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
458   gime = rl->GetLoader("TPCLoader");
459   Stat_t nentries = gime->TreeS()->GetEntries();
460   
461
462   //create branch's in TPC treeD
463   AliSimDigits * digrow = new AliSimDigits;
464   TTree * tree  = ogime->TreeD();
465
466   tree->Branch("Segment","AliSimDigits",&digrow);
467   param->SetZeroSup(2);
468
469   Int_t zerosup = param->GetZeroSup();
470   //Loop over segments of the TPC
471     
472   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
473   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
474   for (Int_t n=0; n<nentries; n++) {
475     rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
476     gime = rl->GetLoader("TPCLoader");
477     gime->TreeS()->GetEvent(n);
478
479     digarr[0]->ExpandBuffer();
480     digarr[0]->ExpandTrackBuffer();
481
482
483     for (Int_t i=1;i<nInputs; i++){ 
484 //      fDigInput->GetInputTreeTPCS(i)->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());      
485       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
486       gime = rl->GetLoader("TPCLoader");
487       gime->TreeS()->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());  
488       digarr[i]->ExpandBuffer();
489       digarr[i]->ExpandTrackBuffer();
490       if ((digarr[0]->GetID()-digarr[i]->GetID())>0) 
491        printf("problem\n");
492     
493     }   
494     
495     Int_t sector, padRow;
496     if (!param->AdjustSectorRow(digarr[0]->GetID(),sector,padRow)) {
497       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr[0]->GetID()<<endl;
498       continue;
499     }
500
501     AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector);  // pad gains per given sector
502     AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector);  // noise per given sector
503     digrow->SetID(digarr[0]->GetID());
504
505     Int_t nTimeBins = digarr[0]->GetNRows();
506     Int_t nPads = digarr[0]->GetNCols();
507     digrow->Allocate(nTimeBins,nPads);
508     digrow->AllocateTrack(3);
509
510     Float_t q=0;
511     Int_t label[1000]; //stack for 300 events 
512     Int_t labptr = 0;
513
514     
515
516     for (Int_t iTimeBin=0;iTimeBin<nTimeBins; iTimeBin++){   // iTimeBin
517       for (Int_t iPad=0;iPad<nPads; iPad++){    // pad
518     
519        q=0;
520        labptr=0;
521        // looop over digits 
522         for (Int_t i=0;i<nInputs; i++){ 
523          q  += digarr[i]->GetDigitFast(iTimeBin,iPad);
524           //q  += *(pdig[i]);
525          
526           for (Int_t tr=0;tr<3;tr++) {
527            Int_t lab = digarr[i]->GetTrackIDFast(iTimeBin,iPad,tr);
528            //Int_t lab = ptr[i][tr*nElems];
529             if ( (lab > 1) ) {
530               label[labptr]=lab+masks[i];
531               labptr++;
532             }          
533           }
534          // pdig[i]++;
535          //ptr[i]++;
536          
537         }
538        q/=16.;  //conversion factor
539        //       Float_t noise  = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());  
540        Float_t gain = gainROC->GetValue(padRow,iPad);
541        q*= gain;
542        Float_t noisePad = noiseROC->GetValue(padRow, iPad);
543
544        Float_t noise  = pTPC->GetNoise();
545        q+=noise*noisePad;
546        //
547        // here we can get digits from past and add signal
548        //
549        //
550        //for (Int_t jTimeBin=0; jTimeBin<iTimeBin; jTimeBin++)
551        //  q+=ionTail
552        //
553
554         q=TMath::Nint(q);
555         if (q > zerosup){ 
556          
557          if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
558          digrow->SetDigitFast((Short_t)q,iTimeBin,iPad);  
559          // *pdig1 =Short_t(q);
560          for (Int_t tr=0;tr<3;tr++){
561            if (tr<labptr) 
562              ((AliSimDigits*)digrow)->SetTrackIDFast(label[tr],iTimeBin,iPad,tr);
563            //ptr1[tr*nElems] = label[tr];
564            //else
565              //           ((AliSimDigits*)digrow)->SetTrackIDFast(-1,iTimeBin,iPad,tr);          
566            //  ptr1[tr*nElems] = 1;
567          }
568        }
569        //pdig1++;
570        //ptr1++;
571     }
572     }
573     
574     digrow->CompresBuffer(1,zerosup);
575     digrow->CompresTrackBuffer(1);
576     tree->Fill();
577     if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";  
578   } 
579 //  printf("end TPC merging - end -Tree %s\t%p\n",fDigInput->GetInputTreeH(0)->GetName(),fDigInput->GetInputTreeH(0)->GetListOfBranches()->At(3));
580   //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
581     ogime->WriteDigits("OVERWRITE");
582
583     for (Int_t i=1;i<nInputs; i++) 
584      { 
585       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
586       gime = rl->GetLoader("TPCLoader");
587       gime->UnloadSDigits();
588      }
589     ogime->UnloadDigits();
590     
591   delete digrow;     
592   for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
593   delete [] masks;
594   delete [] digarr;  
595 }
596
597
598
599
600
601
602
603 //------------------------------------------------------------------------
604 void AliTPCDigitizer::DigitizeWithTailAndCrossTalk(Option_t* option) 
605 {
606   // Modified version of the digitization function
607   // Modification: adding the ion tail and crosstalk:
608   //
609   // pcstream used in order to visually inspect data
610   //
611   //
612   // Crosstalk simulation:
613   //   1.) Calculate per time bin mean charge (per pad)  within anode wire segment 
614   //   2.) Subsract for the clusters at given time bin fraction of (mean charge) normalized by add hoc constant
615   //       AliTPCRecoParam::GetCrosstalkCorrection() (0 if not crosstalk, 1 if ideal crosstalk)
616   //       for simplicity we are assuming that wire segents are related to pad-rows
617   //       Wire segmentationn is obtatined from the      
618   //       AliTPCParam::GetWireSegment(Int_t sector, Int_t row); // to be implemented
619   //       AliTPCParam::GetNPadsPerSegment(Int_t segmentID); // to be implemented
620   //   3.) Substract form the signal contribution from the previous tracks - Ion tail in case speified in the AliTPCRecoParam
621   //        AliTPCRecoParam::GetUseIonTailCorrection()
622   //
623   // Ion tail simulation:
624   //    1.) Needs signal from pad+-1, taking signal from history
625   // merge input tree's with summable digits
626   // output stored in TreeTPCD
627   //  
628   AliTPCcalibDB* const calib=AliTPCcalibDB::Instance();
629   AliTPCRecoParam *recoParam = calib->GetRecoParam(0); 
630   AliDebug(1, Form(" recoParam->GetCrosstalkCorrection()  =   %f", recoParam->GetCrosstalkCorrection())); 
631   AliDebug(1,Form(" recoParam->GetUseIonTailCorrection() =  %d ", recoParam->GetUseIonTailCorrection()));
632   Int_t nROCs = 72;
633   char s[100]; 
634   char ss[100];
635   TString optionString = option;
636   if (!strcmp(optionString.Data(),"deb")) {
637     cout<<"AliTPCDigitizer:::DigitizeFast called with option deb "<<endl;
638     fDebug = 3;
639   }
640
641   // ======== get detector and geometry =======
642
643   AliRunLoader *rl, *orl;
644   AliLoader *gime, *ogime;
645
646   if (gAlice == 0x0)
647   {
648     Warning("DigitizeFast","gAlice is NULL. Loading from input 0");
649     rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
650     if (rl == 0x0)
651     {
652       Error("DigitizeFast","Can not find Run Loader for input 0. Can not proceed.");
653       return;
654     }
655     rl->LoadgAlice();
656     rl->GetAliRun();
657   }
658   AliTPC *pTPC  = (AliTPC *) gAlice->GetModule("TPC");
659   AliTPCParam * param = pTPC->GetParam();
660
661   //sprintf(s,param->GetTitle());
662   snprintf(s,100,"%s",param->GetTitle());
663   //sprintf(ss,"75x40_100x60");
664   snprintf(ss,100,"75x40_100x60");
665   if(strcmp(s,ss)==0){
666     printf("2 pad-length geom hits with 3 pad-lenght geom digits...\n");
667     delete param;
668     param=new AliTPCParamSR();
669   } else {
670     //sprintf(ss,"75x40_100x60_150x60");
671     snprintf(ss,100,"75x40_100x60_150x60");
672     if(strcmp(s,ss)!=0) {
673       printf("No TPC parameters found...\n");
674       exit(2); 
675     }
676   }
677
678   pTPC->GenerNoise(500000); //create table with noise
679   //
680   Int_t nInputs = fDigInput->GetNinputs();
681   Int_t * masks = new Int_t[nInputs];
682   for (Int_t i=0; i<nInputs;i++)
683     masks[i]= fDigInput->GetMask(i);
684   Short_t **pdig= new Short_t*[nInputs];   //pointers to the expanded digits array
685   Int_t **ptr=  new Int_t*[nInputs];       //pointers to the expanded tracks array
686   Bool_t *active=  new Bool_t[nInputs];    //flag for active input segments
687   Char_t phname[100];
688
689   //create digits array for given sectors
690   // make indexes
691   //
692   //create branch's in TPC treeD
693   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
694   ogime = orl->GetLoader("TPCLoader");
695   TTree * tree  = ogime->TreeD();
696   AliSimDigits * digrow = new AliSimDigits;  
697
698   if (tree == 0x0)
699   {
700     ogime->MakeTree("D");
701     tree  = ogime->TreeD();
702   }
703   tree->Branch("Segment","AliSimDigits",&digrow);
704   //  
705   AliSimDigits ** digarr = new AliSimDigits*[nInputs]; 
706   for (Int_t i1=0;i1<nInputs; i1++)
707   {
708     digarr[i1]=0;
709     //    intree[i1]
710     rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
711     gime = rl->GetLoader("TPCLoader");
712     gime->LoadSDigits("read");
713     TTree * treear =  gime->TreeS();
714
715     if (!treear) 
716     {
717       cerr<<"AliTPCDigitizer: Input tree with SDigits not found in"
718         <<" input "<< i1<<endl;
719       for (Int_t i2=0;i2<i1+1; i2++){ 
720
721         if(digarr[i2])  delete digarr[i2];
722       }
723       delete [] digarr;
724       delete [] active;
725       delete []masks;
726       delete []pdig;
727       delete []ptr;
728       return;
729     }
730
731     //sprintf(phname,"lhcphase%d",i1);
732     snprintf(phname,100,"lhcphase%d",i1);
733     TParameter<float> *ph = (TParameter<float>*)treear->GetUserInfo()
734       ->FindObject("lhcphase0");
735     if(!ph)
736     {
737       cerr<<"AliTPCDigitizer: LHC phase  not found in"
738         <<" input "<< i1<<endl;
739       for (Int_t i2=0;i2<i1+1; i2++){ 
740         if(digarr[i2])  delete digarr[i2];
741       }
742       delete [] digarr;
743       delete [] active;
744       delete []masks;
745       delete []pdig;
746       delete []ptr;
747       return;
748     }
749     tree->GetUserInfo()->Add(new TParameter<float>(phname,ph->GetVal()));
750     //
751     if (treear->GetIndex()==0)  
752       treear->BuildIndex("fSegmentID","fSegmentID");
753     treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
754   }
755
756
757
758
759   //
760   // zero supp, take gain and noise map of TPC from OCDB 
761   param->SetZeroSup(2);
762   Int_t zerosup = param->GetZeroSup(); 
763   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor(); 
764   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); 
765
766  
767   // 
768   // Cache the ion tail objects form OCDB
769   //      
770  
771   TObjArray *ionTailArr = (TObjArray*)AliTPCcalibDB::Instance()->GetIonTailArray();
772   if (!ionTailArr) {AliFatal("TPC - Missing IonTail OCDB object");}
773  //  TObject *rocFactorIROC  = ionTailArr->FindObject("factorIROC");
774 //   TObject *rocFactorOROC  = ionTailArr->FindObject("factorOROC");
775 //   Float_t factorIROC      = (atof(rocFactorIROC->GetTitle()));
776 //   Float_t factorOROC      = (atof(rocFactorOROC->GetTitle()));
777   Int_t nIonTailBins =0;
778   TObjArray timeResFunc(nROCs); 
779   for (Int_t isec = 0;isec<nROCs;isec++){        //loop overs sectors
780     // Array of TGraphErrors for a given sector
781     TGraphErrors ** graphRes   = new TGraphErrors *[20];
782     Float_t * trfIndexArr    = new Float_t[20];
783     for (Int_t icache=0; icache<20; icache++)
784     {
785       graphRes[icache]       = NULL;
786       trfIndexArr[icache] = 0;
787     }
788     if (!AliTPCcalibDB::Instance()->GetTailcancelationGraphs(isec,graphRes,trfIndexArr)) continue;
789     
790     // fill all TGraphErrors of trfs (time response function) of a given sector to a TObjArray
791     TObjArray *timeResArr = new TObjArray(20);  timeResArr -> SetOwner(kTRUE); 
792     for (Int_t ires = 0;ires<20;ires++) timeResArr->AddAt(graphRes[ires],ires);
793     timeResFunc.AddAt(timeResArr,isec); // Fill all trfs into a single TObjArray 
794     nIonTailBins = graphRes[3]->GetN();
795   }
796
797   //
798   // 1.) Make first loop to calculate mean amplitude per pad per segment for cross talk 
799   //
800   
801   TObjArray   crossTalkSignalArray(nROCs);  // for 36 sectors 
802   TVectorD  * qTotSectorOld  = new TVectorD(nROCs);
803   Float_t qTotTPC = 0.;
804   Float_t qTotPerSector = 0.;
805   Int_t nTimeBinsAll = 1100;
806   Int_t nWireSegments=11;
807   // 1.a) crorstalk matrix initialization
808   for (Int_t sector=0; sector<nROCs; sector++){
809     TMatrixD *pcrossTalkSignal = new TMatrixD(nWireSegments,nTimeBinsAll);
810     for (Int_t imatrix = 0; imatrix<11; imatrix++)
811       for (Int_t jmatrix = 0; jmatrix<nTimeBinsAll; jmatrix++){
812         (*pcrossTalkSignal)[imatrix][jmatrix]=0.;
813       }
814     crossTalkSignalArray.AddAt(pcrossTalkSignal,sector);
815   }
816   //  
817   // main loop over rows of whole TPC
818   for (Int_t globalRowID=0; globalRowID<param->GetNRowsTotal(); globalRowID++) {
819     Int_t sector, padRow;
820     if (!param->AdjustSectorRow(globalRowID,sector,padRow)) {
821       cerr<<"AliTPC warning: invalid segment ID ! "<<globalRowID<<endl;
822       continue;
823     }
824     // Calculate number of pads in a anode wire segment for normalization
825     Int_t wireSegmentID    = param->GetWireSegment(sector,padRow);
826     Float_t nPadsPerSegment = (Float_t)(param->GetNPadsPerSegment(wireSegmentID));
827     // structure with mean signal per pad to be filled for each timebin in first loop (11 anodeWireSegment and 1100 timebin)
828     TMatrixD &crossTalkSignal =  *((TMatrixD*)crossTalkSignalArray.At(sector));
829     AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector);  // pad gains per given sector
830     //  digrow->SetID(globalRowID);
831     Int_t nTimeBins = 0;
832     Int_t nPads = 0;
833     Bool_t digitize = kFALSE;
834     for (Int_t i=0;i<nInputs; i++){   //here we can have more than one input  - merging of separate events, signal1+signal2+background 
835       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
836       gime = rl->GetLoader("TPCLoader");
837       if (gime->TreeS()->GetEntryWithIndex(globalRowID,globalRowID) >= 0) {
838         digarr[i]->ExpandBuffer();
839         digarr[i]->ExpandTrackBuffer();
840         nTimeBins = digarr[i]->GetNRows();
841         nPads = digarr[i]->GetNCols();
842         active[i] = kTRUE;
843         if (!GetRegionOfInterest() || (i == 0)) digitize = kTRUE;
844       } else {
845         active[i] = kFALSE;
846       }
847       if (GetRegionOfInterest() && !digitize) break;
848     }   
849     if (!digitize) continue;
850     //digrow->Allocate(nTimeBins,nPads);
851     Float_t q    = 0;
852     Int_t labptr = 0;
853     Int_t nElems = nTimeBins*nPads; // element is a unit of a given row's pad-timebin space        
854     for (Int_t i=0;i<nInputs; i++)
855       if (active[i]) { 
856         pdig[i] = digarr[i]->GetDigits();
857       }
858     //    
859     // loop over elements i.e pad-timebin space of a row, "padNumber=elem/nTimeBins", "timeBin=elem%nTimeBins"
860     // 
861     for (Int_t elem=0;elem<nElems; elem++) {    
862       q=0;
863       labptr=0;
864       // looop over digits 
865       for (Int_t i=0;i<nInputs; i++) if (active[i])       { 
866         q  += *(pdig[i]);
867         pdig[i]++;
868       }
869       if (q<=0) continue;   
870       Int_t padNumber = elem/nTimeBins;
871       Int_t timeBin   = elem%nTimeBins;      
872       Float_t gain = gainROC->GetValue(padRow,padNumber);  // get gain for given - pad-row pad
873       q*= gain;
874       crossTalkSignal[wireSegmentID][timeBin]+= q/nPadsPerSegment;        // Qtot per segment for a given timebin
875       qTotSectorOld -> GetMatrixArray()[sector] += q;                      // Qtot for each sector
876       qTotTPC += q;                                                        // Qtot for whole TPC       
877     } // end of q loop
878   } // end of global row loop
879
880   //
881   // 1.b) Dump the content of the crossTalk signal to the debug stremer - to be corealted later with the same crosstalk correction
882   //      assumed during reconstruction 
883   //
884   if (AliTPCReconstructor::StreamLevel()==1) {
885     //
886     // dump the crosstalk matrices to tree for further investigation
887     //     a.) to estimate fluctuation of pedestal in indiviula wire segments
888     //     b.) to check correlation between regions
889     //     c.) to check relative conribution of signal below threshold to crosstalk
890     for (Int_t isector=0; isector<nROCs; isector++){  //set all ellemts of crosstalk matrix to 0
891       TMatrixD * crossTalkMatrix = (TMatrixD*)crossTalkSignalArray.At(isector);
892       //TMatrixD * crossTalkMatrixBelow = (TMatrixD*)fCrossTalkSignalArray->At(isector+nROCs);
893       TVectorD vecAll(crossTalkMatrix->GetNrows());
894       //TVectorD vecBelow(crossTalkMatrix->GetNrows());
895       //
896       for (Int_t itime=0; itime<crossTalkMatrix->GetNcols(); itime++){
897         for (Int_t iwire=0; iwire<crossTalkMatrix->GetNrows(); iwire++){
898           vecAll[iwire]=(*crossTalkMatrix)(iwire,itime);
899           //vecBelow[iwire]=(*crossTalkMatrixBelow)(iwire,itime);
900         }
901         (*fDebugStreamer)<<"crosstalkMatrix"<<
902           "sector="<<isector<<
903           "itime="<<itime<<
904           "vecAll.="<<&vecAll<<                // crosstalk charge + charge below threshold
905           //"vecBelow.="<<&vecBelow<<            // crosstalk contribution from signal below threshold
906           "\n";
907       }
908     }
909   }
910
911
912
913
914   //
915   // 2.) Loop over segments (padrows) of the TPC 
916   //  
917   for (Int_t globalRowID=0; globalRowID<param->GetNRowsTotal(); globalRowID++) {
918     Int_t sector, padRow;
919     if (!param->AdjustSectorRow(globalRowID,sector,padRow)) {
920       cerr<<"AliTPC warning: invalid segment ID ! "<<globalRowID<<endl;
921       continue;
922     }
923     TObjArray *arrTRF = (TObjArray*)timeResFunc.At(sector);  
924     //    TGraphErrors  *graphTRF = (TGraphErrors*)arrTRF->At(1);
925     Int_t wireSegmentID    = param->GetWireSegment(sector,padRow);
926     Float_t nPadsPerSegment = (Float_t)(param->GetNPadsPerSegment(wireSegmentID));
927     //    const Float_t ampfactor = (sector<36)?factorIROC:factorOROC;      // factor for the iontail which is ROC type dependent
928     AliTPCCalROC * gainROC  = gainTPC->GetCalROC(sector);  // pad gains per given sector
929     AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector);  // noise per given sector
930     digrow->SetID(globalRowID);
931     Int_t nTimeBins = 0;
932     Int_t nPads = 0;
933     Bool_t digitize = kFALSE;
934     for (Int_t i=0;i<nInputs; i++) {   //here we can have more than one input  - merging of separate events, signal1+signal2+background
935       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
936       gime = rl->GetLoader("TPCLoader");
937       if (gime->TreeS()->GetEntryWithIndex(globalRowID,globalRowID) >= 0) {
938         digarr[i]->ExpandBuffer();
939         digarr[i]->ExpandTrackBuffer();
940         nTimeBins = digarr[i]->GetNRows();
941         nPads = digarr[i]->GetNCols();
942         active[i] = kTRUE;
943         if (!GetRegionOfInterest() || (i == 0)) digitize = kTRUE;
944       } else {
945         active[i] = kFALSE;
946       }
947       if (GetRegionOfInterest() && !digitize) break;
948     }   
949     if (!digitize) continue;
950     
951     digrow->Allocate(nTimeBins,nPads);
952     digrow->AllocateTrack(3);
953     
954     Int_t localPad = 0;
955     Float_t q    = 0.;
956     Float_t qXtalk   = 0.;
957     Float_t qIonTail = 0.;
958     Float_t qOrig = 0.;
959     Int_t label[1000]; //stack for 300 events 
960     Int_t labptr = 0;
961     Int_t nElems = nTimeBins*nPads; // element is a unit of a given row's pad-timebin space    
962     for (Int_t i=0;i<nInputs; i++)
963       if (active[i]) { 
964         pdig[i] = digarr[i]->GetDigits();
965         ptr[i]  = digarr[i]->GetTracks();
966       }
967     Short_t *pdig1= digrow->GetDigits();
968     Int_t   *ptr1= digrow->GetTracks() ;
969     // loop over elements i.e pad-timebin space of a row
970     for (Int_t elem=0;elem<nElems; elem++)     {     
971       q=0; 
972       labptr=0;
973       // looop over digits 
974       for (Int_t i=0;i<nInputs; i++) if (active[i]){ 
975         q  += *(pdig[i]);
976         for (Int_t tr=0;tr<3;tr++)         {
977           Int_t lab = ptr[i][tr*nElems];
978           if ( (lab > 1) && *(pdig[i])>zerosup) {
979             label[labptr]=lab+masks[i];
980             labptr++;
981           }          
982         }
983         pdig[i]++;
984         ptr[i]++;
985       }
986       Int_t padNumber = elem/nTimeBins;
987       Int_t timeBin   = elem%nTimeBins;
988       localPad = padNumber-nPads/2;
989       
990       Float_t gain = gainROC->GetValue(padRow,padNumber);  // get gain for given - pad-row pad
991       //if (gain<0.5){
992       //printf("problem\n");
993       //}
994       q*= gain;
995       qOrig = q;
996       Float_t noisePad = noiseROC->GetValue(padRow,padNumber);
997       Float_t noise  = pTPC->GetNoise()*noisePad;
998       if ( (q/16.+noise)> zerosup){
999         // Crosstalk correction 
1000         qXtalk = (*(TMatrixD*)crossTalkSignalArray.At(sector))[wireSegmentID][timeBin];
1001         qTotPerSector = qTotSectorOld -> GetMatrixArray()[sector];    
1002         
1003         // Ion tail correction: being elem=padNumber*nTimeBins+timeBin;
1004         Int_t lowerElem=elem-nIonTailBins;    
1005         Int_t zeroElem =(elem/nTimeBins)*nTimeBins;
1006         if (zeroElem<0) zeroElem=0;
1007         if (lowerElem<zeroElem) lowerElem=zeroElem;
1008         // 
1009         qIonTail=0;
1010         if (q>0 && recoParam->GetUseIonTailCorrection()){
1011           for (Int_t i=0;i<nInputs; i++) if (active[i]){ 
1012             Short_t *pdigC= digarr[i]->GetDigits();
1013             if (padNumber==0) continue;
1014             if (padNumber>=nPads-1) continue;
1015             for (Int_t dpad=-1; dpad<=1; dpad++){          // calculate iontail due signals from neigborhood pads
1016               for (Int_t celem=elem-1; celem>lowerElem; celem--){
1017                 Int_t celemPad=celem+dpad*nTimeBins;
1018                 Double_t qCElem=pdigC[celemPad];
1019                 if ( qCElem<=0) continue;
1020                 //here we substract ion tail    
1021                 Double_t COG=0;
1022                 if (celemPad-nTimeBins>nTimeBins && celemPad+nTimeBins<nElems){   // COG calculation in respect to current pad in pad units
1023                   Double_t sumAmp=pdigC[celemPad-nTimeBins]+pdigC[celemPad]+pdigC[celemPad+nTimeBins];
1024                   COG=(-1.0*pdigC[celemPad-nTimeBins]+pdigC[celemPad+nTimeBins])/sumAmp;
1025                 }
1026                 Int_t indexTRFPRF = (TMath::Nint(TMath::Abs(COG*10.))%20);
1027                 TGraphErrors  *graphTRFPRF = (TGraphErrors*)arrTRF->At(indexTRFPRF);
1028                 if (graphTRFPRF==NULL) continue;
1029                 // here we should get index and point of TRF corresponding to given COG position
1030                 if (graphTRFPRF->GetY()[elem-celem]<0)qIonTail+=qCElem*graphTRFPRF->GetY()[elem-celem];
1031               }
1032             }
1033           }
1034         }
1035       }
1036       //
1037       q -= qXtalk*recoParam->GetCrosstalkCorrection();
1038       q+=qIonTail;
1039       q/=16.;                                              //conversion factor
1040       q+=noise; 
1041       q=TMath::Nint(q);  // round to the nearest integer
1042       
1043       
1044       // fill info for degugging
1045       if (AliTPCReconstructor::StreamLevel()==1 && qOrig > zerosup ) {
1046         TTreeSRedirector &cstream = *fDebugStreamer;
1047         UInt_t uid = AliTPCROC::GetTPCUniqueID(sector, padRow, padNumber);
1048         cstream <<"ionTailXtalk"<<
1049           "uid="<<uid<<                        // globla unique identifier
1050           "sector="<< sector<<   
1051           "globalRowID="<<globalRowID<<
1052           "padRow="<< padRow<<                 //pad row
1053           "wireSegmentID="<< wireSegmentID<<   //wire segment 0-11, 0-3 in IROC 4-10 in OROC 
1054           "localPad="<<localPad<<              // pad number -npads/2 .. npads/2
1055           "padNumber="<<padNumber<<            // pad number 0..npads 
1056           "timeBin="<< timeBin<<               // time bin 
1057           "nPadsPerSegment="<<nPadsPerSegment<<// number of pads per wire segment         
1058           "qTotPerSector="<<qTotPerSector<<    // total charge in sector 
1059           //
1060           "qTotTPC="<<qTotTPC<<                // acumulated charge without crosstalk and ion tail in full TPC
1061           "qOrig="<< qOrig<<                   // charge in given pad-row,pad,time-bin
1062           "q="<<q<<                            // q=qOrig-qXtalk-qIonTail - to check sign of the effects
1063           "qXtalk="<<qXtalk<<                  // crosstal contribtion at given position
1064           "qIonTail="<<qIonTail<<              // ion tail cotribution from past signal
1065           "\n";
1066       }// dump the results to the debug streamer if in debug mode
1067       
1068       if (q > zerosup){ 
1069         if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
1070         //digrow->SetDigitFast((Short_t)q,rows,col);  
1071         *pdig1 =Short_t(q);
1072         for (Int_t tr=0;tr<3;tr++)
1073           {
1074           if (tr<labptr) 
1075             ptr1[tr*nElems] = label[tr];
1076         }
1077       }
1078       pdig1++;
1079       ptr1++;
1080     }
1081      
1082    if (AliTPCReconstructor::StreamLevel()==1 && qOrig > zerosup ) { 
1083       cout << " sector = " << sector  << " row = " << padRow << " localPad = " << localPad << " SegmentID = " << wireSegmentID << " nPadsPerSegment = " <<  nPadsPerSegment << endl;
1084       cout << " qXtalk =   " <<  qXtalk ;
1085       cout << " qOrig = " <<  qOrig ;
1086       cout << " q = "    <<  q ;
1087       cout << " qsec = " <<  qTotPerSector ;
1088       cout << " qTotTPC "<<  qTotTPC << endl;
1089    }
1090     //
1091     //  glitch filter
1092     //
1093     digrow->GlitchFilter();
1094     //
1095     digrow->CompresBuffer(1,zerosup);
1096     digrow->CompresTrackBuffer(1);
1097     tree->Fill();
1098     if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";  
1099   } //for (Int_t n=0; n<param->GetNRowsTotal(); n++) 
1100
1101
1102   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
1103   ogime = orl->GetLoader("TPCLoader");
1104   ogime->WriteDigits("OVERWRITE");
1105
1106   //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
1107
1108   delete digrow;     
1109   for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
1110   delete []masks;
1111   delete []pdig;
1112   delete []ptr;
1113   delete []active;
1114   delete []digarr;  
1115 }