PHOS/AliPHOSDApi0mip.cxx

   1 // --
   2 // --
   3 // Implementation for TTree output in PHOS DA
   4 // for calibrating energy by pi0 and MIP.
   5 // --
   6 // -- Author: Hisayuki Torii (Hiroshima Univ.)
   7 // --
   8
   9 #include <stdio.h>
  10 #include <iostream>
  11 #include <math.h>
  12 #include "AliPHOSDApi0mip.h"
  13 ClassImp(AliPHOSDApi0mip)
  14 //----------------------------------------------------------------
  15 AliPHOSDApi0mip::AliPHOSDApi0mip(int module,int iterid,char* fopt) :
  16   TNamed(), fCreateTree(false), fCreateHist(false), fMod(0), fIterId(0),
  17   fTFile(0), fTTree(0), fEvent(0), fEventClustered(0), fTime(0),
  18   fH1Time(0), fH1DigitNum(0), fH1ClusterNum(0),
  19   fH2EneDigitId(0), fH2MipDigitId(0), fH2Pi0DigitId(0), fH3Pi0AsymPt(0) {
  20   // Constructor
  21
  22   char hname[1024], htitle[1024];
  23   sprintf(hname,"AliPHOSDApi0mip_mod%d_iter%d",module,iterid);
  24   SetName(hname);
  25   sprintf(htitle,"PHOS MIP/pi0 Calibration DA for Module:%d Iteration:%d",module,iterid);
  26   SetTitle(htitle);
  27
  28   fMod = module;
  29   fIterId = iterid;
  30   fEvent = 0;
  31   fEventClustered = false;
  32   fTime = 0;
  33   fCreateTree = false;
  34   fCreateHist = false;
  35
  36   char fname[1024];
  37   sprintf(fname,"AliPHOSDApi0mip_mod%d.root",module);
  38   fTFile = TFile::Open(fname,fopt);
  39 }
  40 //----------------------------------------------------------------
  41 AliPHOSDApi0mip::AliPHOSDApi0mip(const AliPHOSDApi0mip& da):
  42   TNamed(da), fCreateTree(false), fCreateHist(false), fMod(da.fMod), fIterId(da.fIterId),
  43   fTFile(0), fTTree(0), fEvent(0), fEventClustered(0), fTime(0),
  44   fH1Time(0), fH1DigitNum(0), fH1ClusterNum(0),
  45   fH2EneDigitId(0), fH2MipDigitId(0), fH2Pi0DigitId(0), fH3Pi0AsymPt(0) {
  46   // Copy Constructor
  47
  48   char fname[1024], hname[1024], htitle[1024];;
  49   sprintf(fname,"%s.root",GetName());
  50   fTFile = TFile::Open(fname,"RECREATE");
  51   fEvent = 0;
  52   fEventClustered = false;
  53   fTime = 0;
  54
  55   if( da.fCreateHist ){
  56     fCreateHist = true;
  57     fH1Time = (TH1I*) da.fH1Time->Clone();
  58     fH1DigitNum = (TH1F*) da.fH1DigitNum->Clone();
  59     fH1ClusterNum = (TH1F*) da.fH1ClusterNum->Clone();
  60     fH2EneDigitId = (TH2F*) da.fH2EneDigitId->Clone();
  61     fH2MipDigitId = (TH2F*) da.fH2MipDigitId->Clone();
  62     fH2Pi0DigitId = (TH2F*) da.fH2Pi0DigitId->Clone();
  63     fH3Pi0AsymPt = (TH3F*) da.fH3Pi0AsymPt->Clone();
  64   } else {
  65     fCreateHist = false;
  66     fH1Time = 0;
  67     fH1DigitNum = 0;
  68     fH1ClusterNum = 0;
  69     fH2EneDigitId = 0;
  70     fH2MipDigitId = 0;
  71     fH2Pi0DigitId = 0;
  72     fH3Pi0AsymPt = 0;
  73   }
  74   if( da.fCreateTree ){
  75     // Create new ttree instead of copy.
  76     fCreateTree = true;
  77     sprintf(hname,"tevt_mod%d_iter%d",fMod,fIterId);
  78     sprintf(htitle,"Calibration for Module:%d Iteration:%d",fMod,fIterId);
  79     fTTree = new TTree(hname,htitle);
  80     fTTree->Branch("AliPHOSDATreeEvent","AliPHOSDATreeEvent",&fEvent);
  81   } else {
  82     fCreateTree = false;
  83     fTTree = 0;
  84   }
  85 }
  86 //-------------------------------------------------------------------
  87 AliPHOSDApi0mip::~AliPHOSDApi0mip(){
  88   // Destructor
  89
  90   //std::cout<<" DEBUGDEBUG :: AliPHOSDApi0mip:: deconstructor "<<std::endl;
  91   if( fCreateHist ){
  92     fH1Time->Write(0,TObject::kOverwrite);
  93     fH1DigitNum->Write(0,TObject::kOverwrite);
  94     fH1ClusterNum->Write(0,TObject::kOverwrite);
  95     fH2EneDigitId->Write(0,TObject::kOverwrite);
  96     fH2MipDigitId->Write(0,TObject::kOverwrite);
  97     fH2Pi0DigitId->Write(0,TObject::kOverwrite);
  98     fH3Pi0AsymPt->Write(0,TObject::kOverwrite);
  99     delete fH1Time;
 100     delete fH1DigitNum;
 101     delete fH1ClusterNum;
 102     delete fH2EneDigitId;
 103     delete fH2MipDigitId;
 104     delete fH2Pi0DigitId;
 105     delete fH3Pi0AsymPt;
 106   }
 107   if( fCreateTree ){
 108     fTTree->Write();
 109   }
 110   if( fEvent ) delete fEvent;
 111   fTFile->Close();
 112   delete fTFile;
 113 }
 114 //-------------------------------------------------------------------
 115 void AliPHOSDApi0mip::FillDigit(float adc,int row,int col){
 116   // Put new digit information in event
 117   //
 118
 119   if( fEvent==0 ) fEvent = new AliPHOSDATreeEvent;
 120   float energy = adc * 0.005; // 5MeV/ch is fixed.
 121   if( energy > 0.020 ){  // Threshold at 20MeV.
 122     //std::cout<<" AliPHOSDApi0mip::AppendDigit() DEBUG:: (row,col)=("<<row<<","<<col<<") adc="<<adc<<" energy="<<energy<<std::endl;
 123     fEvent->Fill(energy,row,col);
 124   }
 125 }
 126 //-------------------------------------------------------------------
 127 void AliPHOSDApi0mip::NewEvent(){
 128   // Initiator must be called in the beginning of event
 129   //
 130
 131   if(fEvent) fEvent->Reset();
 132   fEventClustered = false;
 133 }
 134 //----------------------------------------------------------------
 135 bool AliPHOSDApi0mip::CreateTree(){
 136   if(! fTFile ) {
 137     return false;
 138   }
 139   fTFile->cd();
 140   char hname[1024], htitle[1024];
 141   sprintf(hname,"trevt_mod%d_iter%d",fMod,fIterId);
 142   sprintf(htitle,"Calibration for Module:%d Iteration:%d",fMod,fIterId);
 143   fTTree = new TTree(hname,htitle);
 144   fTTree->Branch("AliPHOSDATreeEvent","AliPHOSDATreeEvent",&fEvent);
 145   fCreateTree = true;
 146   return true;
 147 }
 148 //----------------------------------------------------------------
 149 bool AliPHOSDApi0mip::CreateHist(){
 150   // Create histogram routine called by the constructor only.
 151
 152   char hname[1024], htitle[1024];
 153
 154   sprintf(hname,"h1_time_mod%d_iter%d",fMod,fIterId);
 155   sprintf(htitle,"Time : Mod:%d Iter:%d",fMod,fIterId);
 156   fH1Time = (TH1I*) gDirectory->Get(hname);
 157   if( fH1Time>0 ){
 158     std::cout<<" AliPHOSDApi0mip:Warning!! Output object already exist : "<<fH1Time->GetName()<<std::endl;
 159     return false;
 160   }
 161   fH1Time = new TH1I(hname,htitle,2,0,2);
 162
 163   sprintf(hname,"h1_digitnum_mod%d_iter%d",fMod,fIterId);
 164   sprintf(htitle,"Number of Digits : Mod:%d Iter:%d",fMod,fIterId);
 165   fH1DigitNum = new TH1F(hname,htitle,100,0,100);
 166
 167   sprintf(hname,"h1_clusternum_mod%d_iter%d",fMod,fIterId);
 168   sprintf(htitle,"Number of Clusters : Mod:%d Iter:%d",fMod,fIterId);
 169   fH1ClusterNum = new TH1F(hname,htitle,100,0,100);
 170
 171   sprintf(hname,"h2_pi0digitid_mod%d_iter%d",fMod,fIterId);
 172   sprintf(htitle,"PHOS pi0 mass vs Digit Id : Mod:%d Iter:%d",fMod,fIterId);
 173   fH2Pi0DigitId = new TH2F(hname,htitle,17920,0,17920,120,0,0.3);
 174
 175   sprintf(hname,"h2_mipdigitid_mod%d_iter%d",fMod,fIterId);
 176   sprintf(htitle,"PHOS MIP vs Digit Id : Mod:%d Iter:%d",fMod,fIterId);
 177   fH2MipDigitId = new TH2F(hname,htitle,17920,0,17920,100,0-0.0025,0.5-0.0025);
 178   //fH2MipDigitId = new TH2F(hname,htitle,17920,0,17920,50,0.5,1.5);
 179
 180   sprintf(hname,"h2_enedigitid_mod%d_iter%d",fMod,fIterId);
 181   sprintf(htitle,"PHOS MIP vs Digit Id : Mod:%d Iter:%d",fMod,fIterId);
 182   //fH2EneDigitId = new TH2F(hname,htitle,17920,0,17920,100,0-0.0025,0.5-0.0025);
 183   fH2EneDigitId = new TH2F(hname,htitle,17920,0,17920,50,0,5.0);
 184
 185   sprintf(hname,"h3_pi0asympt_mod%d_iter%d",fMod,fIterId);
 186   sprintf(htitle,"PHOS pi0 mass vs pT vs Asym : Mod:%d Iter:%d",fMod,fIterId);
 187   fH3Pi0AsymPt = new TH3F(hname,htitle,20,0,1,20,0,10,200,0,1);
 188
 189   fCreateHist = true;
 190   return true;
 191 }
 192 //-------------------------------------------------------------------
 193 void AliPHOSDApi0mip::FillTree(AliPHOSDATreeEvent* event){
 194   // Fill all information (AliPHOSDATreeEvent) into TTree
 195   // if event==NULL, the obtained information (=fEvent) is used.
 196
 197   if(! fCreateTree ) CreateTree();
 198   if( event == 0 ){
 199     if( fEvent == 0 ){
 200       fEvent = new AliPHOSDATreeEvent();
 201     }
 202     // No need to make clustering for tree saving.
 203     //if(! fEventClustered ){
 204     //fEvent->ExecuteClustering();
 205     //fEventClustered = true;
 206     //}
 207     fEvent->SetTime(fTime);
 208     fTTree->Fill();
 209   } else {
 210     if( fEvent ) delete fEvent;
 211     fEvent = event;
 212     fTTree->Fill();
 213     fEvent = 0;
 214   }
 215 }
 216 //-------------------------------------------------------------------
 217 void AliPHOSDApi0mip::FillHist(AliPHOSDATreeEvent* event){
 218   // Fill all information (AliPHOSDATreeEvent) into histograms
 219   // if event==NULL, the obtained information (=fEvent) is used.
 220
 221   if(! fCreateHist ) CreateHist();
 222   if( event == 0 ){
 223     if( fEvent==0 ) fEvent = new AliPHOSDATreeEvent;
 224     if(! fEventClustered ){
 225       fEvent->ExecuteClustering();
 226       fEventClustered = true;
 227     }
 228     fEvent->SetTime(fTime);
 229     event = fEvent;
 230   }
 231
 232   // -- Constant
 233   int nclt1, nclt2;
 234   int ndigits;
 235   int ndigitsall = 0;
 236   AliPHOSDATreeDigit digit1, digit2;
 237   const float kDistanceFromVertex = 460.;
 238   float dist1, dist2, px, py;
 239   float mass, cosine, pt, asym;
 240
 241   // -- Filling into Time Histogram
 242   double start = fH1Time->GetBinContent(1);
 243   double stop = fH1Time->GetBinContent(2);
 244   double currenttime = (double) (event->GetTime());
 245   if( start == 0 && stop == 0 ) start = stop = currenttime;
 246   if( currenttime < start ) start = currenttime;
 247   if( currenttime > stop ) stop = currenttime;
 248   fH1Time->SetBinContent(1,start);
 249   fH1Time->SetBinContent(2,stop);
 250
 251   // -- Start looping clusters
 252   nclt1 = event->GetNClusters();
 253   fH1ClusterNum->Fill(nclt1);
 254   //std::cout<<" AliPHOSDApi0mip::FillHisto() DEBUG:: Number of clusters: "<<nclt1<<std::endl;
 255   while( nclt1-- ){
 256     AliPHOSDATreeCluster& clt1 = event->GetCluster(nclt1);
 257     ndigits = clt1.GetNDigits();
 258     while( ndigits-- ){
 259       AliPHOSDATreeDigit& digit = clt1.GetDigit(ndigits);
 260       fH2EneDigitId->Fill(digit.GetDigitId(),digit.GetEnergy());
 261       ndigitsall++;
 262     }
 263     digit1 = clt1.GetMaxDigit();
 264     if( clt1.GetEnergy() > 0 && clt1.GetNDigits() >= 2 ){
 265       fH2MipDigitId->Fill(digit1.GetDigitId(),clt1.GetEnergy());
 266     }
 267     if( clt1.GetEnergy() > 0.2 ){
 268       nclt2 = nclt1;
 269       while( nclt2-- ){
 270         AliPHOSDATreeCluster& clt2 = event->GetCluster(nclt2);
 271         digit2 = clt2.GetMaxDigit();
 272         if( clt2.GetEnergy() > 0.2 ){
 273           //
 274           // Following approximation is applied in order to avoid many calling "sqrt" function.
 275           // sqrt( 1 + x^2 ) =~ 1 + 0.5x^2 (for x<<1)
 276           //
 277           //cosine = ((clt1.GetRow()-31.5)*(clt2.GetRow()-31.5)*4.84 + (clt1.GetCol()-27.5)*(clt2.GetCol()-27.5)*4.84 + kDistanceFromVertex*kDistanceFromVertex)
 278           //        / ( kDistanceFromVertex + (clt1.GetRow()-31.5)*(clt1.GetRow()-31.5)*2.42/kDistanceFromVertex + (clt1.GetCol()-27.5)*(clt1.GetCol()-27.5)*2.42/kDistanceFromVertex )
 279           //        / ( kDistanceFromVertex + (clt2.GetRow()-31.5)*(clt2.GetRow()-31.5)*2.42/kDistanceFromVertex + (clt2.GetCol()-27.5)*(clt2.GetCol()-27.5)*2.42/kDistanceFromVertex );
 280           //
 281
 282
 283           dist1 = sqrt( kDistanceFromVertex*kDistanceFromVertex + (clt1.GetRow()-31.5)*(clt1.GetRow()-31.5)*4.84 + (clt1.GetCol()-27.5)*(clt1.GetCol()-27.5)*4.84 );
 284           dist2 = sqrt( kDistanceFromVertex*kDistanceFromVertex + (clt2.GetRow()-31.5)*(clt2.GetRow()-31.5)*4.84 + (clt2.GetCol()-27.5)*(clt2.GetCol()-27.5)*4.84 );
 285           cosine = ((clt1.GetRow()-31.5)*(clt2.GetRow()-31.5)*4.84 + (clt1.GetCol()-27.5)*(clt2.GetCol()-27.5)*4.84 + kDistanceFromVertex*kDistanceFromVertex)
 286             / dist1 / dist2;
 287           mass = sqrt( 2.0 * clt1.GetEnergy() * clt2.GetEnergy() * (1 - cosine) );
 288           //std::cout<<" DEBUG: dist1:"<<dist1<<" dist2:"<<dist2<<" cosine:"<<cosine<<" c1:"<<clt1.GetEnergy()<<" c2:"<<clt2.GetEnergy()<<" mass:"<<mass<<std::endl;
 289           //std::cout<<"      : clt1(row,col)=("<<clt1.GetRow()<<","<<clt1.GetCol()<<") clt2(row,col)=("<<clt2.GetRow()<<","<<clt2.GetCol()<<") "<<std::endl;
 290           asym = fabs(clt1.GetEnergy() - clt2.GetEnergy())/(clt1.GetEnergy() + clt2.GetEnergy());
 291           py = clt1.GetEnergy() * kDistanceFromVertex / dist1 + clt2.GetEnergy() * kDistanceFromVertex / dist2;
 292           px = clt1.GetEnergy() * (clt1.GetCol()-27.5)*2.42/dist1 + clt2.GetEnergy() * (clt2.GetCol()-27.5)*2.42/ dist2;
 293           pt = sqrt(px*px+py*py);
 294           fH3Pi0AsymPt->Fill(asym,pt,mass);
 295           if( digit1.IsValid() && digit1.IsValid() &&
 296               digit1.GetDigitId() != digit2.GetDigitId() &&
 297               //pt > 1.5 && clt1.GetEnergy() > 0.3 && clt2.GetEnergy() > 0.3 ){
 298               //pt > 1.3 && clt1.GetEnergy() > 0.25 && clt2.GetEnergy() > 0.25 ){
 299               pt > 1.0 && clt1.GetEnergy() > 0.2 && clt2.GetEnergy() > 0.2 ){
 300             fH2Pi0DigitId->Fill(digit1.GetDigitId(),mass);
 301             fH2Pi0DigitId->Fill(digit2.GetDigitId(),mass);
 302           }
 303         } // End of energy cut
 304       } // End of second cluster loop (clt2)
 305     } // End of energy cut
 306   } // End of first cluster loop (clt1)
 307   fH1DigitNum->Fill(ndigitsall);
 308   //
 309 }
 310 //-------------------------------------------------------------------
 311 void AliPHOSDApi0mip::Print(char* opt){
 312   // Print Out
 313
 314   //fTFile->ls();
 315   //fTTree->Print();
 316   if( fEvent ){
 317     //fEvent->ExecuteClustering();
 318     fEvent->Print(opt);
 319   }
 320 }
 321 //-------------------------------------------------------------------