TPC/Sim/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
  61 using std::cout;
  62 using std::cerr;
  63 using std::endl;
  64 ClassImp(AliTPCDigitizer)
  65
  66 //___________________________________________
  67   AliTPCDigitizer::AliTPCDigitizer() :AliDigitizer(),fDebug(0), fDebugStreamer(0)
  68 {
  69   //
  70 // Default ctor - don't use it
  71 //
  72
  73 }
  74
  75 //___________________________________________
  76 AliTPCDigitizer::AliTPCDigitizer(AliDigitizationInput* digInput)
  77   :AliDigitizer(digInput),fDebug(0), fDebugStreamer(0)
  78 {
  79   //
  80 // ctor which should be used
  81 //
  82   AliDebug(2,"(AliDigitizationInput* digInput) was processed");
  83   if (AliTPCReconstructor::StreamLevel()>0)  fDebugStreamer = new TTreeSRedirector("TPCDigitDebug.root");
  84
  85 }
  86
  87 //------------------------------------------------------------------------
  88 AliTPCDigitizer::~AliTPCDigitizer()
  89 {
  90 // Destructor
  91   if (fDebugStreamer) delete fDebugStreamer;
  92 }
  93
  94
  95
  96 //------------------------------------------------------------------------
  97 Bool_t AliTPCDigitizer::Init()
  98 {
  99 // Initialization
 100
 101  return kTRUE;
 102 }
 103
 104
 105 //------------------------------------------------------------------------
 106 void AliTPCDigitizer::Digitize(Option_t* option)
 107 {
 108   DigitizeFast(option);
 109 }
 110 //------------------------------------------------------------------------
 111 void AliTPCDigitizer::DigitizeFast(Option_t* option)
 112 {
 113
 114   // merge input tree's with summable digits
 115   //output stored in TreeTPCD
 116   char s[100];
 117   char ss[100];
 118   TString optionString = option;
 119   if (!strcmp(optionString.Data(),"deb")) {
 120     cout<<"AliTPCDigitizer:::DigitizeFast called with option deb "<<endl;
 121     fDebug = 3;
 122   }
 123   //get detector and geometry
 124
 125
 126   AliRunLoader *rl, *orl;
 127   AliLoader *gime, *ogime;
 128
 129   if (gAlice == 0x0)
 130    {
 131      Warning("DigitizeFast","gAlice is NULL. Loading from input 0");
 132      rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 133      if (rl == 0x0)
 134       {
 135         Error("DigitizeFast","Can not find Run Loader for input 0. Can not proceed.");
 136         return;
 137       }
 138      rl->LoadgAlice();
 139      rl->GetAliRun();
 140    }
 141   AliTPC *pTPC  = (AliTPC *) gAlice->GetModule("TPC");
 142   AliTPCParam * param = pTPC->GetParam();
 143
 144   //sprintf(s,param->GetTitle());
 145   snprintf(s,100,"%s",param->GetTitle());
 146   //sprintf(ss,"75x40_100x60");
 147   snprintf(ss,100,"75x40_100x60");
 148   if(strcmp(s,ss)==0){
 149     printf("2 pad-length geom hits with 3 pad-lenght geom digits...\n");
 150     delete param;
 151     param=new AliTPCParamSR();
 152   }
 153   else{
 154     //sprintf(ss,"75x40_100x60_150x60");
 155     snprintf(ss,100,"75x40_100x60_150x60");
 156    if(strcmp(s,ss)!=0) {
 157      printf("No TPC parameters found...\n");
 158      exit(2);
 159    }
 160   }
 161
 162   pTPC->GenerNoise(500000); //create table with noise
 163   //
 164   Int_t nInputs = fDigInput->GetNinputs();
 165   Int_t * masks = new Int_t[nInputs];
 166   for (Int_t i=0; i<nInputs;i++)
 167     masks[i]= fDigInput->GetMask(i);
 168   Short_t **pdig= new Short_t*[nInputs];   //pointers to the expanded digits array
 169   Int_t **ptr=  new Int_t*[nInputs];       //pointers to the expanded tracks array
 170   Bool_t *active=  new Bool_t[nInputs];    //flag for active input segments
 171   Char_t phname[100];
 172
 173   //create digits array for given sectors
 174   // make indexes
 175   //
 176   //create branch's in TPC treeD
 177   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
 178   ogime = orl->GetLoader("TPCLoader");
 179   TTree * tree  = ogime->TreeD();
 180   AliSimDigits * digrow = new AliSimDigits;
 181
 182   if (tree == 0x0)
 183    {
 184      ogime->MakeTree("D");
 185      tree  = ogime->TreeD();
 186    }
 187   tree->Branch("Segment","AliSimDigits",&digrow);
 188   //
 189   AliSimDigits ** digarr = new AliSimDigits*[nInputs];
 190   for (Int_t i1=0;i1<nInputs; i1++)
 191     {
 192       digarr[i1]=0;
 193      //    intree[i1]
 194       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
 195       gime = rl->GetLoader("TPCLoader");
 196       gime->LoadSDigits("read");
 197       TTree * treear =  gime->TreeS();
 198
 199       if (!treear)
 200        {
 201         cerr<<"AliTPCDigitizer: Input tree with SDigits not found in"
 202             <<" input "<< i1<<endl;
 203         for (Int_t i2=0;i2<i1+1; i2++){
 204
 205           if(digarr[i2])  delete digarr[i2];
 206         }
 207         delete [] digarr;
 208         delete [] active;
 209         delete []masks;
 210         delete []pdig;
 211         delete []ptr;
 212         return;
 213        }
 214
 215       //sprintf(phname,"lhcphase%d",i1);
 216       snprintf(phname,100,"lhcphase%d",i1);
 217       TParameter<float> *ph = (TParameter<float>*)treear->GetUserInfo()
 218                                ->FindObject("lhcphase0");
 219       if(!ph){
 220         cerr<<"AliTPCDigitizer: LHC phase  not found in"
 221             <<" input "<< i1<<endl;
 222         for (Int_t i2=0;i2<i1+1; i2++){
 223           if(digarr[i2])  delete digarr[i2];
 224         }
 225         delete [] digarr;
 226         delete [] active;
 227         delete []masks;
 228         delete []pdig;
 229         delete []ptr;
 230         return;
 231       }
 232       tree->GetUserInfo()->Add(new TParameter<float>(phname,ph->GetVal()));
 233               //
 234       if (treear->GetIndex()==0)
 235         treear->BuildIndex("fSegmentID","fSegmentID");
 236       treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
 237     }
 238
 239
 240
 241
 242   //
 243
 244   param->SetZeroSup(2);
 245
 246   Int_t zerosup = param->GetZeroSup();
 247   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
 248   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 249   //
 250   //Loop over segments of the TPC
 251
 252   for (Int_t segmentID=0; segmentID<param->GetNRowsTotal(); segmentID++)
 253    {
 254     Int_t sector, padRow;
 255     if (!param->AdjustSectorRow(segmentID,sector,padRow))
 256      {
 257       cerr<<"AliTPC warning: invalid segment ID ! "<<segmentID<<endl;
 258       continue;
 259      }
 260     AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector);  // pad gains per given sector
 261     AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector);  // noise per given sector
 262     digrow->SetID(segmentID);
 263
 264     Int_t nTimeBins = 0;
 265     Int_t nPads = 0;
 266
 267     Bool_t digitize = kFALSE;
 268     for (Int_t i=0;i<nInputs; i++)
 269      {
 270
 271       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
 272       gime = rl->GetLoader("TPCLoader");
 273
 274       if (gime->TreeS()->GetEntryWithIndex(segmentID,segmentID) >= 0) {
 275         digarr[i]->ExpandBuffer();
 276         digarr[i]->ExpandTrackBuffer();
 277         nTimeBins = digarr[i]->GetNRows();
 278         nPads = digarr[i]->GetNCols();
 279         active[i] = kTRUE;
 280         if (!GetRegionOfInterest() || (i == 0)) digitize = kTRUE;
 281       } else {
 282         active[i] = kFALSE;
 283       }
 284       if (GetRegionOfInterest() && !digitize) break;
 285      }
 286     if (!digitize) continue;
 287
 288     digrow->Allocate(nTimeBins,nPads);
 289     digrow->AllocateTrack(3);
 290
 291     Float_t q=0;
 292     Int_t label[1000]; //stack for 300 events
 293     Int_t labptr = 0;
 294
 295     Int_t nElems = nTimeBins*nPads;
 296
 297     for (Int_t i=0;i<nInputs; i++)
 298      if (active[i]) {
 299        pdig[i] = digarr[i]->GetDigits();
 300        ptr[i]  = digarr[i]->GetTracks();
 301       }
 302
 303     Short_t *pdig1= digrow->GetDigits();
 304     Int_t   *ptr1= digrow->GetTracks() ;
 305
 306
 307
 308     for (Int_t elem=0;elem<nElems; elem++)
 309      {
 310
 311        q=0;
 312        labptr=0;
 313        // looop over digits
 314         for (Int_t i=0;i<nInputs; i++) if (active[i])
 315          {
 316           //          q  += digarr[i]->GetDigitFast(rows,col);
 317             q  += *(pdig[i]);
 318
 319            for (Int_t tr=0;tr<3;tr++)
 320             {
 321              //             Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
 322              Int_t lab = ptr[i][tr*nElems];
 323              if ( (lab > 1) && *(pdig[i])>zerosup)
 324               {
 325                 label[labptr]=lab+masks[i];
 326                 labptr++;
 327               }
 328             }
 329            pdig[i]++;
 330            ptr[i]++;
 331          }
 332         q/=16.;  //conversion factor
 333         Float_t gain = gainROC->GetValue(padRow,elem/nTimeBins);  // get gain for given - pad-row pad
 334         //if (gain<0.5){
 335           //printf("problem\n");
 336         //}
 337         q*= gain;
 338         Float_t noisePad = noiseROC->GetValue(padRow,elem/nTimeBins);
 339         //       Float_t noise  = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());
 340         Float_t noise  = pTPC->GetNoise();
 341         q+=noise*noisePad;
 342         q=TMath::Nint(q);
 343         if (q > zerosup)
 344          {
 345           if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
 346           //digrow->SetDigitFast((Short_t)q,rows,col);
 347           *pdig1 =Short_t(q);
 348           for (Int_t tr=0;tr<3;tr++)
 349            {
 350             if (tr<labptr)
 351              ptr1[tr*nElems] = label[tr];
 352            }
 353           }
 354         pdig1++;
 355         ptr1++;
 356      }
 357     //
 358     //  glitch filter
 359     //
 360     digrow->GlitchFilter();
 361     //
 362     digrow->CompresBuffer(1,zerosup);
 363     digrow->CompresTrackBuffer(1);
 364     tree->Fill();
 365     if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";
 366    } //for (Int_t n=0; n<param->GetNRowsTotal(); n++)
 367
 368
 369   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
 370   ogime = orl->GetLoader("TPCLoader");
 371   ogime->WriteDigits("OVERWRITE");
 372
 373   //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
 374
 375   delete digrow;
 376   for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
 377   delete []masks;
 378   delete []pdig;
 379   delete []ptr;
 380   delete []active;
 381   delete []digarr;
 382 }
 383
 384
 385
 386 //------------------------------------------------------------------------
 387 void AliTPCDigitizer::DigitizeSave(Option_t* option)
 388 {
 389   //
 390   // Merge input tree's with summable digits
 391   // Output digits stored in TreeTPCD
 392   //
 393   // Not active for long time.
 394   // Before adding modification (for ion tail calucation and for the crorsstalk) it should be
 395   //  checked one by one with currenlty used AliTPCDigitizer::DigitizeFast
 396   //
 397   TString optionString = option;
 398   if (!strcmp(optionString.Data(),"deb")) {
 399     cout<<"AliTPCDigitizer::Digitize: called with option deb "<<endl;
 400     fDebug = 3;
 401   }
 402   //get detector and geometry
 403   AliRunLoader *rl, *orl;
 404   AliLoader *gime, *ogime;
 405
 406
 407   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
 408   ogime = orl->GetLoader("TPCLoader");
 409
 410   rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 411   //gime = rl->GetLoader("TPCLoader");
 412   rl->GetLoader("TPCLoader");
 413   rl->LoadgAlice();
 414   AliRun* alirun = rl->GetAliRun();
 415
 416   AliTPC *pTPC  = (AliTPC *) alirun->GetModule("TPC");
 417   AliTPCParam * param = pTPC->GetParam();
 418   pTPC->GenerNoise(500000); //create teble with noise
 419   printf("noise %f \n",  param->GetNoise()*param->GetNoiseNormFac());
 420   //
 421   Int_t nInputs = fDigInput->GetNinputs();
 422   // stupid protection...
 423   if (nInputs <= 0) return;
 424   //
 425   Int_t * masks = new Int_t[nInputs];
 426   for (Int_t i=0; i<nInputs;i++)
 427     masks[i]= fDigInput->GetMask(i);
 428
 429   AliSimDigits ** digarr = new AliSimDigits*[nInputs];
 430   for(Int_t ii=0;ii<nInputs;ii++) digarr[ii]=0;
 431
 432   for (Int_t i1=0;i1<nInputs; i1++)
 433    {
 434      //digarr[i1]=0;
 435     //    intree[i1]
 436     rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
 437     gime = rl->GetLoader("TPCLoader");
 438
 439     TTree * treear =  gime->TreeS();
 440     //
 441     if (!treear) {
 442       cerr<<" TPC -  not existing input = \n"<<i1<<" ";
 443       delete [] masks;
 444       for(Int_t i=0; i<nInputs; i++) delete digarr[i];
 445       delete [] digarr;
 446       return;
 447     }
 448     //
 449     TBranch * br = treear->GetBranch("fSegmentID");
 450     if (br) br->GetFile()->cd();
 451     treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
 452    }
 453
 454   rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 455   gime = rl->GetLoader("TPCLoader");
 456   Stat_t nentries = gime->TreeS()->GetEntries();
 457
 458
 459   //create branch's in TPC treeD
 460   AliSimDigits * digrow = new AliSimDigits;
 461   TTree * tree  = ogime->TreeD();
 462
 463   tree->Branch("Segment","AliSimDigits",&digrow);
 464   param->SetZeroSup(2);
 465
 466   Int_t zerosup = param->GetZeroSup();
 467   //Loop over segments of the TPC
 468
 469   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
 470   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 471   for (Int_t n=0; n<nentries; n++) {
 472     rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 473     gime = rl->GetLoader("TPCLoader");
 474     gime->TreeS()->GetEvent(n);
 475
 476     digarr[0]->ExpandBuffer();
 477     digarr[0]->ExpandTrackBuffer();
 478
 479
 480     for (Int_t i=1;i<nInputs; i++){
 481 //      fDigInput->GetInputTreeTPCS(i)->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());
 482       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
 483       gime = rl->GetLoader("TPCLoader");
 484       gime->TreeS()->GetEntryWithIndex(digarr[0]->GetID(),digarr[0]->GetID());
 485       digarr[i]->ExpandBuffer();
 486       digarr[i]->ExpandTrackBuffer();
 487       if ((digarr[0]->GetID()-digarr[i]->GetID())>0)
 488        printf("problem\n");
 489
 490     }
 491
 492     Int_t sector, padRow;
 493     if (!param->AdjustSectorRow(digarr[0]->GetID(),sector,padRow)) {
 494       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr[0]->GetID()<<endl;
 495       continue;
 496     }
 497
 498     AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector);  // pad gains per given sector
 499     AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector);  // noise per given sector
 500     digrow->SetID(digarr[0]->GetID());
 501
 502     Int_t nTimeBins = digarr[0]->GetNRows();
 503     Int_t nPads = digarr[0]->GetNCols();
 504     digrow->Allocate(nTimeBins,nPads);
 505     digrow->AllocateTrack(3);
 506
 507     Float_t q=0;
 508     Int_t label[1000]; //stack for 300 events
 509     Int_t labptr = 0;
 510
 511
 512
 513     for (Int_t iTimeBin=0;iTimeBin<nTimeBins; iTimeBin++){   // iTimeBin
 514       for (Int_t iPad=0;iPad<nPads; iPad++){    // pad
 515
 516        q=0;
 517        labptr=0;
 518        // looop over digits
 519         for (Int_t i=0;i<nInputs; i++){
 520          q  += digarr[i]->GetDigitFast(iTimeBin,iPad);
 521           //q  += *(pdig[i]);
 522
 523           for (Int_t tr=0;tr<3;tr++) {
 524            Int_t lab = digarr[i]->GetTrackIDFast(iTimeBin,iPad,tr);
 525            //Int_t lab = ptr[i][tr*nElems];
 526             if ( (lab > 1) ) {
 527               label[labptr]=lab+masks[i];
 528               labptr++;
 529             }
 530           }
 531          // pdig[i]++;
 532          //ptr[i]++;
 533
 534         }
 535        q/=16.;  //conversion factor
 536        //       Float_t noise  = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());
 537        Float_t gain = gainROC->GetValue(padRow,iPad);
 538        q*= gain;
 539        Float_t noisePad = noiseROC->GetValue(padRow, iPad);
 540
 541        Float_t noise  = pTPC->GetNoise();
 542        q+=noise*noisePad;
 543        //
 544        // here we can get digits from past and add signal
 545        //
 546        //
 547        //for (Int_t jTimeBin=0; jTimeBin<iTimeBin; jTimeBin++)
 548        //  q+=ionTail
 549        //
 550
 551         q=TMath::Nint(q);
 552         if (q > zerosup){
 553
 554          if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
 555          digrow->SetDigitFast((Short_t)q,iTimeBin,iPad);
 556          // *pdig1 =Short_t(q);
 557          for (Int_t tr=0;tr<3;tr++){
 558            if (tr<labptr)
 559              ((AliSimDigits*)digrow)->SetTrackIDFast(label[tr],iTimeBin,iPad,tr);
 560            //ptr1[tr*nElems] = label[tr];
 561            //else
 562              //           ((AliSimDigits*)digrow)->SetTrackIDFast(-1,iTimeBin,iPad,tr);
 563            //  ptr1[tr*nElems] = 1;
 564          }
 565        }
 566        //pdig1++;
 567        //ptr1++;
 568     }
 569     }
 570
 571     digrow->CompresBuffer(1,zerosup);
 572     digrow->CompresTrackBuffer(1);
 573     tree->Fill();
 574     if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";
 575   }
 576 //  printf("end TPC merging - end -Tree %s\t%p\n",fDigInput->GetInputTreeH(0)->GetName(),fDigInput->GetInputTreeH(0)->GetListOfBranches()->At(3));
 577   //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
 578     ogime->WriteDigits("OVERWRITE");
 579
 580     for (Int_t i=1;i<nInputs; i++)
 581      {
 582       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
 583       gime = rl->GetLoader("TPCLoader");
 584       gime->UnloadSDigits();
 585      }
 586     ogime->UnloadDigits();
 587
 588   delete digrow;
 589   for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
 590   delete [] masks;
 591   delete [] digarr;
 592 }
 593
 594
 595
 596 //------------------------------------------------------------------------
 597 void AliTPCDigitizer::DigitizeWithTailAndCrossTalk(Option_t* option)
 598 {
 599   // Modified version of the digitization function
 600   // Modification: adding the ion tail and crosstalk:
 601   //
 602   // pcstream used in order to visually inspect data
 603   //
 604   //
 605   // Crosstalk simulation:
 606   //   1.) Calculate per time bin mean charge (per pad)  within anode wire segment
 607   //   2.) Subsract for the clusters at given time bin fraction of (mean charge) normalized by add hoc constant
 608   //       AliTPCRecoParam::GetCrosstalkCorrection() (0 if not crosstalk, 1 if ideal crosstalk)
 609   //       for simplicity we are assuming that wire segents are related to pad-rows
 610   //       Wire segmentationn is obtatined from the
 611   //       AliTPCParam::GetWireSegment(Int_t sector, Int_t row); // to be implemented
 612   //       AliTPCParam::GetNPadsPerSegment(Int_t segmentID); // to be implemented
 613   //
 614   // Ion tail simulation:
 615   //    1.) Needs signal from pad+-1, taking signal from history
 616   // merge input tree's with summable digits
 617   // output stored in TreeTPCD
 618   //
 619   char s[100];
 620   char ss[100];
 621   TString optionString = option;
 622   if (!strcmp(optionString.Data(),"deb")) {
 623     cout<<"AliTPCDigitizer:::DigitizeFast called with option deb "<<endl;
 624     fDebug = 3;
 625   }
 626   //get detector and geometry
 627
 628
 629   AliRunLoader *rl, *orl;
 630   AliLoader *gime, *ogime;
 631
 632   if (gAlice == 0x0)
 633    {
 634      Warning("DigitizeFast","gAlice is NULL. Loading from input 0");
 635      rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(0));
 636      if (rl == 0x0)
 637       {
 638         Error("DigitizeFast","Can not find Run Loader for input 0. Can not proceed.");
 639         return;
 640       }
 641      rl->LoadgAlice();
 642      rl->GetAliRun();
 643    }
 644   AliTPC *pTPC  = (AliTPC *) gAlice->GetModule("TPC");
 645   AliTPCParam * param = pTPC->GetParam();
 646
 647   //sprintf(s,param->GetTitle());
 648   snprintf(s,100,"%s",param->GetTitle());
 649   //sprintf(ss,"75x40_100x60");
 650   snprintf(ss,100,"75x40_100x60");
 651   if(strcmp(s,ss)==0){
 652     printf("2 pad-length geom hits with 3 pad-lenght geom digits...\n");
 653     delete param;
 654     param=new AliTPCParamSR();
 655   }
 656   else{
 657     //sprintf(ss,"75x40_100x60_150x60");
 658     snprintf(ss,100,"75x40_100x60_150x60");
 659    if(strcmp(s,ss)!=0) {
 660      printf("No TPC parameters found...\n");
 661      exit(2);
 662    }
 663   }
 664
 665   pTPC->GenerNoise(500000); //create table with noise
 666   //
 667   Int_t nInputs = fDigInput->GetNinputs();
 668   Int_t * masks = new Int_t[nInputs];
 669   for (Int_t i=0; i<nInputs;i++)
 670     masks[i]= fDigInput->GetMask(i);
 671   Short_t **pdig= new Short_t*[nInputs];   //pointers to the expanded digits array
 672   Int_t **ptr=  new Int_t*[nInputs];       //pointers to the expanded tracks array
 673   Bool_t *active=  new Bool_t[nInputs];    //flag for active input segments
 674   Char_t phname[100];
 675
 676   //create digits array for given sectors
 677   // make indexes
 678   //
 679   //create branch's in TPC treeD
 680   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
 681   ogime = orl->GetLoader("TPCLoader");
 682   TTree * tree  = ogime->TreeD();
 683   AliSimDigits * digrow = new AliSimDigits;
 684
 685   if (tree == 0x0)
 686    {
 687      ogime->MakeTree("D");
 688      tree  = ogime->TreeD();
 689    }
 690   tree->Branch("Segment","AliSimDigits",&digrow);
 691   //
 692   AliSimDigits ** digarr = new AliSimDigits*[nInputs];
 693   for (Int_t i1=0;i1<nInputs; i1++)
 694     {
 695       digarr[i1]=0;
 696      //    intree[i1]
 697       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i1));
 698       gime = rl->GetLoader("TPCLoader");
 699       gime->LoadSDigits("read");
 700       TTree * treear =  gime->TreeS();
 701
 702       if (!treear)
 703        {
 704         cerr<<"AliTPCDigitizer: Input tree with SDigits not found in"
 705             <<" input "<< i1<<endl;
 706         for (Int_t i2=0;i2<i1+1; i2++){
 707
 708           if(digarr[i2])  delete digarr[i2];
 709         }
 710         delete [] digarr;
 711         delete [] active;
 712         delete []masks;
 713         delete []pdig;
 714         delete []ptr;
 715         return;
 716        }
 717
 718       //sprintf(phname,"lhcphase%d",i1);
 719       snprintf(phname,100,"lhcphase%d",i1);
 720       TParameter<float> *ph = (TParameter<float>*)treear->GetUserInfo()
 721                                ->FindObject("lhcphase0");
 722       if(!ph){
 723         cerr<<"AliTPCDigitizer: LHC phase  not found in"
 724             <<" input "<< i1<<endl;
 725         for (Int_t i2=0;i2<i1+1; i2++){
 726           if(digarr[i2])  delete digarr[i2];
 727         }
 728         delete [] digarr;
 729         delete [] active;
 730         delete []masks;
 731         delete []pdig;
 732         delete []ptr;
 733         return;
 734       }
 735       tree->GetUserInfo()->Add(new TParameter<float>(phname,ph->GetVal()));
 736               //
 737       if (treear->GetIndex()==0)
 738         treear->BuildIndex("fSegmentID","fSegmentID");
 739       treear->GetBranch("Segment")->SetAddress(&digarr[i1]);
 740     }
 741
 742
 743
 744
 745   //
 746
 747   param->SetZeroSup(2);
 748   Int_t zerosup = param->GetZeroSup();
 749   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
 750   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 751   //
 752   // 1.) Make first loop to calculate mean amplitude per pad per segment
 753   //
 754   // TObjArray * crossTalkSignalArray(nsectors);
 755   // TMatrixD  crossTalkSignal(nWireSegments, timeBin);  // structure with mean signal per pad to be filled in first loop
 756   //
 757   //
 758   //2.) Loop over segments (padrows) of the TPC
 759   //
 760   for (Int_t segmentID=0; segmentID<param->GetNRowsTotal(); segmentID++)
 761    {
 762     Int_t sector, padRow;
 763     if (!param->AdjustSectorRow(segmentID,sector,padRow))
 764      {
 765       cerr<<"AliTPC warning: invalid segment ID ! "<<segmentID<<endl;
 766       continue;
 767      }
 768     AliTPCCalROC * gainROC = gainTPC->GetCalROC(sector);  // pad gains per given sector
 769     AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(sector);  // noise per given sector
 770     digrow->SetID(segmentID);
 771
 772     Int_t nTimeBins = 0;
 773     Int_t nPads = 0;
 774
 775     Bool_t digitize = kFALSE;
 776     for (Int_t i=0;i<nInputs; i++)
 777       {
 778
 779       rl = AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(i));
 780       gime = rl->GetLoader("TPCLoader");
 781
 782       if (gime->TreeS()->GetEntryWithIndex(segmentID,segmentID) >= 0) {
 783         digarr[i]->ExpandBuffer();
 784         digarr[i]->ExpandTrackBuffer();
 785         nTimeBins = digarr[i]->GetNRows();
 786         nPads = digarr[i]->GetNCols();
 787         active[i] = kTRUE;
 788         if (!GetRegionOfInterest() || (i == 0)) digitize = kTRUE;
 789       } else {
 790         active[i] = kFALSE;
 791       }
 792       if (GetRegionOfInterest() && !digitize) break;
 793      }
 794     if (!digitize) continue;
 795
 796     digrow->Allocate(nTimeBins,nPads);
 797     digrow->AllocateTrack(3);
 798
 799     Float_t q=0;
 800     Int_t label[1000]; //stack for 300 events
 801     Int_t labptr = 0;
 802
 803     Int_t nElems = nTimeBins*nPads;
 804
 805     for (Int_t i=0;i<nInputs; i++)
 806      if (active[i]) {
 807        pdig[i] = digarr[i]->GetDigits();
 808        ptr[i]  = digarr[i]->GetTracks();
 809       }
 810
 811     Short_t *pdig1= digrow->GetDigits();
 812     Int_t   *ptr1= digrow->GetTracks() ;
 813
 814
 815
 816     for (Int_t elem=0;elem<nElems; elem++)
 817      {
 818        // padNumber=elem/nTimeBins
 819        // timeBin=elem%nTimeBins
 820        q=0;
 821        labptr=0;
 822        // looop over digits
 823         for (Int_t i=0;i<nInputs; i++) if (active[i])
 824          {
 825           //          q  += digarr[i]->GetDigitFast(rows,col);
 826             q  += *(pdig[i]);
 827
 828            for (Int_t tr=0;tr<3;tr++)
 829             {
 830              //             Int_t lab = digarr[i]->GetTrackIDFast(rows,col,tr);
 831              Int_t lab = ptr[i][tr*nElems];
 832              if ( (lab > 1) && *(pdig[i])>zerosup)
 833               {
 834                 label[labptr]=lab+masks[i];
 835                 labptr++;
 836               }
 837             }
 838            pdig[i]++;
 839            ptr[i]++;
 840          }
 841         Int_t padNumber=elem/nTimeBins;
 842
 843         q/=16.;  //conversion factor
 844         Float_t gain = gainROC->GetValue(padRow,padNumber);  // get gain for given - pad-row pad
 845         //if (gain<0.5){
 846           //printf("problem\n");
 847         //}
 848         q*= gain;
 849         // Crosstalk correction:
 850         //      Double_t qCrossTalk=crossTalkSignal(wireSegment,timeBin); // signal matrix from per sector array
 851         //
 852         //
 853         // Ion tail correction:
 854         //   padNumber=elem/nTimeBins
 855         //   timeBin=elem%nTimeBins
 856         //   elem=padNumber*nTimeBins+timeBin;
 857         //    lowerElem=elem-nIonTailBins;
 858         //    if (lowerElem<0) lowerElem=0;
 859         //    if (lowerElem in previospad) lowerElem = padNumber*nTimeBins;
 860         //
 861         // for (Int_t celem=elem-1; celem>lowerElem; celem--){
 862         //  Int_t deltaT=elem-celem
 863         //
 864         // }
 865         //
 866         Float_t noisePad = noiseROC->GetValue(padRow,padNumber);
 867         //       Float_t noise  = gRandom->Gaus(0,param->GetNoise()*param->GetNoiseNormFac());
 868         Float_t noise  = pTPC->GetNoise();
 869         q+=noise*noisePad;
 870         q=TMath::Nint(q);  // round to the nearest integer
 871         /*
 872           fill infor to check consistency of the data
 873          */
 874         if (q > zerosup)
 875          {
 876           if(q >= param->GetADCSat()) q = (Short_t)(param->GetADCSat() - 1);
 877           //digrow->SetDigitFast((Short_t)q,rows,col);
 878           *pdig1 =Short_t(q);
 879           for (Int_t tr=0;tr<3;tr++)
 880            {
 881             if (tr<labptr)
 882              ptr1[tr*nElems] = label[tr];
 883            }
 884           }
 885         pdig1++;
 886         ptr1++;
 887      }
 888     //
 889     //  glitch filter
 890     //
 891     digrow->GlitchFilter();
 892     //
 893     digrow->CompresBuffer(1,zerosup);
 894     digrow->CompresTrackBuffer(1);
 895     tree->Fill();
 896     if (fDebug>0) cerr<<sector<<"\t"<<padRow<<"\n";
 897    } //for (Int_t n=0; n<param->GetNRowsTotal(); n++)
 898
 899
 900   orl = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
 901   ogime = orl->GetLoader("TPCLoader");
 902   ogime->WriteDigits("OVERWRITE");
 903
 904   //fDigInput->GetTreeDTPC()->Write(0,TObject::kOverwrite);
 905
 906   delete digrow;
 907   for (Int_t i1=0;i1<nInputs; i1++) delete digarr[i1];
 908   delete []masks;
 909   delete []pdig;
 910   delete []ptr;
 911   delete []active;
 912   delete []digarr;
 913 }