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