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