TPC/AliTPCDigitizer.cxx

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