TPC/AliTPCDigitsH.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 /*
  17 $Log$
  18 Revision 1.1.4.2  2000/04/10 11:37:42  kowal2
  19
  20 Digits handling in a new data structure
  21
  22 */
  23
  24 //-----------------------------------------------------------------------------
  25 //
  26 //
  27 //  Author:   Marian Ivanov
  28 //
  29 //  Implementation of class TTPCDigitsH
  30 //
  31 //-----------------------------------------------------------------------------
  32 ////////////////////////////////////////////////
  33 //  Manager class for AliTPCDigitsH             //
  34 ////////////////////////////////////////////////
  35
  36 #include <iostream.h>
  37 #include "TMath.h"
  38 // other include files follow here
  39 #include "TFile.h"
  40 //*KEEP,TFile.
  41 #include "TROOT.h"
  42 #include "TGraph.h"
  43 #include "AliRun.h"
  44 #include "AliDisplay.h"
  45 #include "TPad.h"
  46 #include "TCanvas.h"
  47 #include "TCanvasImp.h"
  48 #include "TPaveText.h"
  49 //*KEEP,TH1.
  50 #include "TH1.h"
  51 //*KEEP,TH2
  52 #include "AliH2F.h"
  53 //*KEEP,TF2.
  54 #include "TF2.h"
  55 //*KEEP,TClonesArray,T=C++.
  56 #include "TClonesArray.h"
  57 #include "TTree.h"
  58 //GALICE includes
  59 #include "GParticle.h"
  60 #include "AliTPC.h"
  61 #include "AliTPCParam.h"
  62 #include "AliTPCD.h"
  63 #include "AliTPCDigitsH.h"
  64
  65
  66 R__EXTERN TSystem *  gSystem;
  67 R__EXTERN AliRun *   gAlice;
  68
  69
  70
  71 ClassImp(AliTPCDigitsH)
  72
  73 AliTPCDigitsH::AliTPCDigitsH()
  74 {
  75 //Begin_Html
  76 /*
  77 <img src="gif/TPCDigitsH.gif">
  78 */
  79 //End_Html
  80
  81   fsec = 1;
  82   frow = 1;
  83   fpad = 1;
  84   fTimeN = 500;
  85   fTimeStart = 0;
  86   fTimeStop = 500;
  87   fOccuN    =25;
  88   fbDelHisto = kFALSE;
  89   fEventN = 0;  //event nuber connected to digit tree
  90   fThreshold  = 5;
  91
  92   fParticles = 0;
  93   fDigits= 0;
  94
  95   fDParam  =0;
  96
  97   fbIOState  = kFALSE;
  98   fbDigState = kFALSE;
  99 }
 100
 101 AliTPCDigitsH::~AliTPCDigitsH()
 102 {
 103
 104 }
 105
 106 AliTPCDigitsH::AliTPCDigitsH(const AliTPCDigitsH &)
 107 {
 108 }
 109
 110 AliTPCDigitsH & AliTPCDigitsH::operator = (const AliTPCDigitsH &)
 111 {
 112    return *this;
 113 }
 114 void  AliTPCDigitsH::SetDParam(AliTPCD * dig)
 115 {
 116   if (dig!=0)
 117     {
 118       fDParam =dig;
 119       fTPCParam = &(fDParam->GetParam());
 120     }
 121   else
 122     {
 123       fTPCParam = 0;
 124       fDParam =0;
 125     }
 126 }
 127
 128 AliTPCParam *&  AliTPCDigitsH::GetParam()
 129 {
 130   return fTPCParam;
 131 }
 132
 133
 134 void AliTPCDigitsH::CloseFiles()
 135 {
 136  if (fin)
 137     {
 138       fin->Close();
 139       //try
 140       // {
 141       delete fin;
 142       //  //    }
 143       //  //catch(...)
 144       //  //    {
 145       //  cout<<"FIN Delete error. Contact autor of root \n";
 146       //        };
 147       fin = 0;
 148     };
 149  if (fout)
 150     {
 151       fout->Close();
 152       //try
 153       //{
 154       //  delete fout;
 155       //}
 156       //      catch(...)
 157       //{
 158       //  cout<<"Out Delete error. Contact autor of root \n";
 159       //};
 160       fout =0;
 161     };
 162  fbIOState  = kFALSE;
 163  fbDigState = kFALSE;
 164 }
 165
 166 Bool_t AliTPCDigitsH::SetIO(const char *  inFile, const char* outFile )
 167 {
 168   ///   Set input and output file
 169   ///   control permisions and so on
 170   ///   if all is OK then set flag fbIOState = kTRUE
 171  fbIOState = kFALSE;
 172  TString s1=inFile;
 173  TString s2=outFile;
 174  //  important ---- it close previious open file if this file exist
 175  if (fin)
 176    {
 177      if (fin == fout) fout = NULL;
 178      fin->Close();
 179      //     try
 180      //{
 181      delete fin;
 182      //}
 183      //catch(...)
 184      //{
 185      //  cout<<"Fin  Delete error. Contact autor of root \n";
 186      //};
 187      fin = 0;
 188    };
 189  //  important ---- it close previous open file if this file exist
 190  if (fout)
 191    {
 192      fout->Close();
 193      //try
 194      //{
 195      delete fout;
 196      //}
 197      //catch(...)
 198      //{
 199      //  cout<<"Fout  Delete error. Contact autor of root \n";
 200      //};
 201      fout = 0;
 202    };
 203
 204  //close the files if exist in root enviroment
 205  TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject((char*)inFile);
 206  if (file) file->Close();
 207  file = (TFile*)gROOT->GetListOfFiles()->FindObject((char*)outFile);
 208  if (file) file->Close();
 209  //if input file is the output file
 210   if (s1 == s2 )
 211     {
 212       if (gSystem->AccessPathName((char*)inFile, kWritePermission ) == kFALSE)
 213          fin  = new TFile((char*)inFile,"UPDATE","trees with digits");
 214       if (!(fin))
 215        {
 216          cout<<"Input file couldn't be open for writing \n";
 217          cout<<"Loook if the file exiest or if you have permision to write \n";
 218          return  kFALSE;
 219        }
 220       else
 221         {
 222           fout = fin;
 223           fbIOState = kTRUE;
 224           return kTRUE;
 225         }
 226     }
 227   //open input file
 228   if (gSystem->AccessPathName((char*)inFile, kReadPermission ) == kFALSE)
 229     fin = new TFile((char*)inFile,"UPDATE","trees with digits");
 230   else
 231     {
 232       cout<<"Input file couldn't be open\n";
 233       cout<<"Maybe the file is not open \n";
 234       return kFALSE ;
 235     }
 236   if (!(fin))
 237     {
 238       cout<<"Input file couldn't be open\n";
 239       cout<<"Probably not root file \n";
 240       return kFALSE ;
 241     }
 242
 243   //open output file
 244   if (gSystem->AccessPathName((char*)outFile, kWritePermission) == kFALSE)
 245     fout = new TFile((char *)outFile,"UPDATE");
 246   else
 247     fout = new TFile((char *)outFile,"NEW");
 248
 249   if (!(fout))
 250     {
 251       cout<<"Output  file couldn't be open\n";
 252       return kFALSE ;
 253     }
 254   //if input and output file is OK set state variable to true
 255   fbIOState = kTRUE;
 256   return kTRUE;
 257 }
 258
 259 Bool_t AliTPCDigitsH::SetEventN(Int_t EventN=0)
 260 {
 261   if (!(fin))
 262     {
 263       cout<<"Warning: Input file not open !!! \n";
 264       return kFALSE;
 265     }
 266   fin->cd();
 267   if (EventN>-1) fEventN = EventN;
 268   fParticles = 0;
 269   if (gAlice)
 270     {
 271       delete gAlice;
 272       gAlice =0;
 273     }
 274   gAlice = (AliRun*) fin->Get("gAlice");
 275   if (gAlice == 0 )
 276     cout<<" Warning : AliRun objects not found in input file \n.";
 277   else
 278     {
 279       gAlice->GetEvent(fEventN);
 280       fParticles = gAlice->Particles();
 281     }
 282   if (fParticles == 0)
 283     return kFALSE;
 284   else
 285     return kTRUE;
 286 }
 287
 288
 289
 290 Bool_t AliTPCDigitsH::SetTree(Int_t eventn  )
 291 {
 292    fbDigState = kFALSE;
 293    ftree = 0;
 294    if (  fbIOState == kFALSE)
 295      {
 296        cout<<"IO files not adjusted \n FIX IT!!!";
 297        return kFALSE;
 298      }
 299    fin->cd();
 300    if (fDParam->SetTree(eventn)==0){
 301      fbDigState = kFALSE;
 302      cout<<"Input tree doesn't exist !!! \n";
 303      return kFALSE;
 304    }
 305    fDigits = fDParam->GetArray();
 306    ftree = fDParam->GetTree();
 307    fbDigState=kTRUE;
 308    return kTRUE;
 309 }
 310
 311 void AliTPCDigitsH::SetSecRowTime(Int_t sec , Int_t row , Int_t TimeN, Float_t TimeStart, Float_t TimeStop )
 312 {
 313   fsec = sec;
 314   frow = row;
 315   fTimeN =TimeN;
 316   fTimeStart = TimeStart;
 317   fTimeStop = TimeStop;
 318 }
 319
 320 void   AliTPCDigitsH::SetParticles(Int_t sec = -1, Int_t row = -1 ,
 321                        Int_t size1 = 30000,Int_t size2=300,
 322                        Bool_t all=kTRUE )
 323 {
 324   if (sec>0) fsec = sec;
 325   if (row>-1) frow =row;
 326   char s[80];
 327   char sh[80];
 328   //  create particles histograms
 329   sprintf(s,"Sector %d   Row %d\n",fsec,frow);
 330   sprintf(sh,"Particles%d_%d",fsec,frow);
 331   fHParticles = new TH1F(sh,s,size1,4,size1);
 332
 333   sprintf(s,"Sector %d   Row %d\n",fsec,frow);
 334   sprintf(sh,"All particles%d_%d",fsec,frow);
 335   fHAllP = new TH1F(sh,s,200,1,25);
 336
 337   sprintf(s,"Sector %d   Row %d\n",fsec,frow);
 338   sprintf(sh,"Secondary Particles%d_%d",fsec,frow);
 339   fHSecondaryP = new TH1F(sh,s,200,1,25);
 340
 341   if (!(fin))
 342   {
 343     cout<<"Input  file not open, open file before \n";
 344   }
 345   else
 346 {
 347   fin->cd();
 348   Int_t sectors_by_rows=(Int_t)ftree->GetEntries();
 349   GParticle * particle;
 350   //  loop over all sectors and rows
 351   for (Int_t n=0; n<sectors_by_rows; n++)
 352     {
 353     if (!ftree->GetEvent(n)) continue;
 354     AliTPCdigit *dig=(AliTPCdigit*)fDigits->UncheckedAt(0);
 355     if (fsec  < dig->fSector) break;
 356     if (fsec != dig->fSector) continue;
 357     if (frow != dig->fPadRow) continue;
 358
 359     Int_t ndigits=fDigits->GetEntriesFast();
 360     //loop over all digits  in sector pad
 361     for (Int_t ndig=0; ndig<ndigits; ndig++)
 362         {
 363       Float_t x,y;
 364         dig=(AliTPCdigit*)fDigits->UncheckedAt(ndig);
 365         fHParticles->Fill(dig->fTracks[0]);
 366         //     get pointer to particle information and fill All and secondary histo
 367         particle = (GParticle*) fParticles->UncheckedAt(dig->fTracks[0]);
 368       fHAllP->Fill(particle->GetKF());
 369       //                Int_t id = particle->GetKF();
 370       x = (Float_t)particle->GetVx();
 371       y = (Float_t)particle->GetVy();
 372       if ( (x*x+y*y ) > 1 )
 373           fHSecondaryP->Fill(particle->GetKF());
 374         if (all==kTRUE)
 375           {
 376             fHParticles->Fill(dig->fTracks[1]);
 377             particle = (GParticle*) fParticles->UncheckedAt(dig->fTracks[1]);
 378             fHAllP->Fill(particle->GetKF());
 379             x = (Float_t)particle->GetVx();
 380             y = (Float_t)particle->GetVy();
 381             if ( (x*x+y*y ) > 1 )
 382               fHSecondaryP->Fill(particle->GetKF());
 383             fHParticles->Fill(dig->fTracks[2]);
 384             particle = (GParticle*) fParticles->UncheckedAt(dig->fTracks[2]);
 385             fHAllP->Fill(particle->GetKF());
 386             x = (Float_t)particle->GetVx();
 387             y = (Float_t)particle->GetVy();
 388             if ( (x*x+y*y ) > 1 )
 389               fHSecondaryP->Fill(particle->GetKF());
 390           }
 391         }
 392       }
 393     //make histogram with multiplicity
 394       char s[80];
 395       char sh[80];
 396       if (all==kTRUE)
 397         sprintf(s,"Number of AliDigits over threshold  per one track in sector %d   Row %d\n (all three most important track recorded)",fsec,frow);
 398       else
 399         sprintf(s,"Number of AliDigits over threshold per sector %d   Row %d\n (only most important track)",fsec,frow);
 400       sprintf(sh,"His_%d_%d",fsec,frow);
 401       fHPartMultiplicity = new TH1F(sh,s,size2,1,size2);
 402       for (Int_t i=1;i<size1;i++)
 403         {
 404           Int_t mul=Int_t(fHParticles->GetBinContent(i));
 405           if (mul>0) fHPartMultiplicity->Fill(mul);
 406         }
 407       if (fout)
 408         {
 409           fout->cd();
 410           fHParticles->Write();
 411           fHPartMultiplicity->Write();
 412           //fHSecondaryP->Write();
 413           //fHAllP->Write();
 414         }
 415   }
 416 }
 417
 418
 419 void AliTPCDigitsH::Anal()
 420 {
 421   if (fbIOState == kFALSE)
 422     {
 423       cout<<"Input output problem. \n Do you initialize IO files ? \n";
 424       return;
 425     }
 426   if (fbDigState == kFALSE)
 427     {
 428       cout<<"Input file doesn't enhalt selected tree \n";
 429       return;
 430     }
 431   fout->cd();
 432
 433 //if we dont want let histogram in memory then we delete old histogram
 434   if ( (fH2Digit) && (fbDelHisto == kTRUE) )
 435     {
 436       //      try
 437       //{
 438           // delete fH2Digit;
 439       //}
 440       //      catch(...)
 441       //{
 442       //  cout<<"Delete error. Contact autor of root \n";
 443       //};
 444       fH2Digit = 0;
 445     }
 446   char s[80];
 447   char sh[80];
 448   sprintf(s,"Sector %d   Row %d\n",fsec,frow);
 449   sprintf(sh,"h%d_%d",fsec,frow);
 450
 451   if  ( (fout) && (fbDelHisto == kFALSE) )
 452     {
 453       fH2Digit = (AliH2F *) fout->Get(sh);
 454       if (fH2Digit) return;
 455     }
 456
 457   Int_t n_of_pads =fTPCParam->GetNPads(fsec,frow);
 458
 459   fH2Digit = new AliH2F(sh, s, fTimeN, fTimeStart, fTimeStop, n_of_pads, 0, n_of_pads-1);
 460
 461   if (!(fout))
 462   {
 463     cout<<"Input  file not open, open file before \n";
 464   }
 465   else
 466   {
 467     //fin->cd();
 468     Int_t sectors_by_rows=(Int_t)ftree->GetEntries();
 469     //loop over all sectors and rows
 470     for (Int_t n=0; n<sectors_by_rows; n++) {
 471       if (!ftree->GetEvent(n)) continue;
 472       AliTPCdigit *dig=(AliTPCdigit*)fDigits->UncheckedAt(0);
 473       if (fsec  < dig->fSector) break;
 474       if (fsec != dig->fSector) continue;
 475       if (frow != dig->fPadRow) continue;
 476
 477       Int_t ndigits=fDigits->GetEntriesFast();
 478       //loop over all digits  in sector pad
 479       for (Int_t ndig=0; ndig<ndigits; ndig++) {
 480         dig=(AliTPCdigit*)fDigits->UncheckedAt(ndig);
 481         fH2Digit->Fill(dig->fTime,dig->fPad,dig->fSignal);
 482       }
 483     }
 484     if (fout) fout->cd();
 485     fH2Digit->Write();
 486   }
 487
 488 }
 489 void  AliTPCDigitsH::Draw(Option_t * opt1 ="cont1"  , Option_t * opt2 = "error",
 490                         Option_t * opt3 = "L" )
 491 {
 492   TString o1 = opt1;
 493   o1.ToLower();
 494
 495   TString o2 = opt2;
 496   o2.ToLower();
 497
 498   TString o3 = opt3;
 499   o3.ToLower();
 500
 501   fcanvas  = new TCanvas("dh","Digits Histograms",700,900);
 502
 503   if (fTitle)
 504     {
 505       //       try
 506       //{
 507           //delete fTitle;
 508       //}
 509       //      catch(...)
 510       //{
 511       //  cout<<"Delete error. Contact autor of root \n";
 512       //};
 513       fTitle = 0;
 514     }
 515   fTitle = new TPaveText(0.2,0.96,0.8,0.995);
 516   fTitle->AddText("Occupancy calculation for TPC");
 517   fTitle->Draw();
 518
 519   fpad1 = new TPad("pad1","",0.05,0.7,0.95,0.95,21);
 520   fpad1->Draw();
 521   fpad2 = new TPad("pad2","",0.05,0.4,0.95,0.65,21);
 522   fpad2->Draw();
 523   fpad3 = new TPad("pad3","",0.05,0.05,0.95,0.35,21);
 524   fpad3->Draw();
 525
 526   fpad1->cd();
 527   //  pad1->TPaveText::title.SetSize(0.1);
 528   if (fH2Digit)
 529     {
 530       fH2DigitBW->Draw(o1);
 531       fH2DigitBW->SetXTitle("time bin");
 532       fH2DigitBW->SetYTitle("pad number");
 533     }
 534   fpad2->cd();
 535   fH1Occu->Fit("pol0");
 536   if (fH1Occu)
 537     {
 538       fH1Occu->Draw(o2);
 539       fH1Occu->SetXTitle("time bin");
 540       fH1Occu->SetYTitle("occupancy");
 541     }
 542   fpad3->cd();
 543   if (fH1Digit)
 544     {
 545       fH1Digit->Draw(o3);
 546       fH1Digit->SetXTitle("time bin");
 547       fH1Digit->SetYTitle("ADC amplitude ");
 548     }
 549
 550 };
 551
 552 void  AliTPCDigitsH::DeleteHisto(const Text_t *namecycle)
 553 {
 554   if (fout) fout->Delete(namecycle);
 555 }
 556
 557 void  AliTPCDigitsH::SetHisto(Int_t pad = 1 )
 558 {
 559
 560   Int_t n_of_pads = fTPCParam->GetNPads(fsec,frow);
 561   if (pad > (n_of_pads-1))
 562     {
 563       cout<<"Pad number is greater then actula number of pads in thi row \n";
 564       cout<<"Noch einmal \n";
 565       return;
 566     }
 567
 568   fpad = pad;
 569   Anal();
 570
 571   //  if ( (fH1Digit) && (fbDelHisto == kTRUE))
 572   //    {
 573   //      try
 574   //        {
 575   //      // delete fH1Digit;
 576   //    }
 577   //      catch(...)
 578   //    {
 579   //      cout<<"Delete error. Contact autor of root \n";
 580   //    };
 581   //      fH1Digit = 0;
 582   //    }
 583
 584   char s[80];
 585   char sh[80];
 586   sprintf(s,"example sector %d   Row %d  Pad %d",fsec,frow,fpad);
 587   sprintf(sh,"h%d_%d_%d",fsec,frow,fpad);
 588  fH2DigitBW = new AliH2F("bw", "", fTimeN, fTimeStart, fTimeStop, n_of_pads, 0, n_of_pads-1);
 589   fH1Digit = new TH1F(sh,s,fTimeN,fTimeStart,fTimeStop);
 590
 591   for (Int_t i = 0;i<fTimeN;i++)
 592     {
 593       Int_t index = fH2Digit->GetBin(i,pad);
 594       Float_t weight = fH2Digit->GetBinContent(index);
 595       fH1Digit->Fill(i,weight);
 596     };
 597
 598
 599   sprintf(s,"Occupancy in sector %d   Row %d  threshold = %d",fsec,frow,fThreshold);
 600   sprintf(sh,"hoccu%d_%d_%d",fsec,frow,fpad);
 601   fH1Occu = new TH1F(sh,s,fOccuN,fTimeStart,fTimeStop);
 602
 603   for (Int_t i = 0;i<fOccuN;i++)
 604     {
 605       Int_t over =0;
 606       Int_t all  =0;
 607       for (int itime = i*(fTimeN/fOccuN); itime<(i+1)*(fTimeN/fOccuN);itime++)
 608         {
 609           for (Int_t ipad = 0; ipad < n_of_pads; ipad++)
 610             {
 611               Int_t index = fH2Digit->GetBin(itime,ipad);
 612               if ( (ipad>3) && ((ipad+3)<n_of_pads)){
 613                 all++;
 614                 if (fH2Digit->GetBinContent(index) >fThreshold) over++ ;
 615               }
 616               if (fH2Digit->GetBinContent(index) >fThreshold)
 617                 fH2DigitBW->Fill(itime,ipad,1);
 618               else
 619                 fH2DigitBW->Fill(itime,ipad,0);
 620             }
 621         }
 622      Float_t occu = ((Float_t)over) /((Float_t) (all));
 623      //     Float_t time =   ((fTimeStop-fTimeStart)/fOccuN)*i+fTimeStart;
 624
 625      fH1Occu->SetBinContent(i,occu);
 626      // Int_t index = fH1Occu->GetBin(i);
 627      Float_t error = sqrt( ((Float_t) ((over)/25+1)) )/((Float_t)(all)/25.);
 628      fH1Occu->SetBinError(i,error);
 629
 630     };
 631
 632
 633 }
 634
 635
 636
 637 void AliTPCDigitsH::Streamer(TBuffer & R__b)
 638 {
 639   if (R__b.IsReading()) {
 640     //      Version_t R__v = R__b.ReadVersion();
 641    } else {
 642       R__b.WriteVersion(AliTPCDigitsH::IsA());
 643    }
 644 }