PWGCF/FEMTOSCOPY/AliFemtoUser/AliFemtoChi2CorrFctn.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 ///                                                                          ///
   3 /// AliFemtoChi2CorrFctn - A correlation function that saves the correlation ///
   4 /// function as a function of single track quality (chi2/ndof) for its and   ///
   5 /// tpc                                                                      ///
   6 /// Authors: Adam Kisiel kisiel@mps.ohio-state.edu                           ///
   7 ///                                                                          ///
   8 ////////////////////////////////////////////////////////////////////////////////
   9
  10 #include "AliFemtoChi2CorrFctn.h"
  11 //#include "AliFemtoHisto.hh"
  12 #include <cstdio>
  13
  14 #ifdef __ROOT__
  15 ClassImp(AliFemtoChi2CorrFctn)
  16 #endif
  17
  18 //____________________________
  19 AliFemtoChi2CorrFctn::AliFemtoChi2CorrFctn(char* title, const int& nbins, const float& QinvLo, const float& QinvHi):
  20   AliFemtoCorrFctn(),
  21   fChi2ITSSUMNumerator(0),
  22   fChi2ITSSUMDenominator(0),
  23   fChi2TPCSUMNumerator(0),
  24   fChi2TPCSUMDenominator(0),
  25   fChi2ITSONENumerator(0),
  26   fChi2ITSONEDenominator(0),
  27   fChi2TPCONENumerator(0),
  28   fChi2TPCONEDenominator(0),
  29   fSigmaToVertexNumerator(0),
  30   fSigmaToVertexDenominator(0)
  31 {
  32   // set up numerator
  33   char tTitNum[101] = "NumChi2ITSSUM";
  34   strncat(tTitNum,title, 100);
  35   fChi2ITSSUMNumerator = new TH2D(tTitNum,title,nbins,QinvLo,QinvHi,200,0.0,6.0);
  36   // set up denominator
  37   char tTitDen[101] = "DenChi2ITSSUM";
  38   strncat(tTitDen,title, 100);
  39   fChi2ITSSUMDenominator = new TH2D(tTitDen,title,nbins,QinvLo,QinvHi,200,0.0,6.0);
  40
  41   // set up numerator
  42   char tTit2Num[101] = "NumChi2TPCSUM";
  43   strncat(tTit2Num,title, 100);
  44   fChi2TPCSUMNumerator = new TH2D(tTit2Num,title,nbins,QinvLo,QinvHi,200,0.0,6.0);
  45   // set up denominator
  46   char tTit2Den[101] = "DenChi2TPCSUM";
  47   strncat(tTit2Den,title, 100);
  48   fChi2TPCSUMDenominator = new TH2D(tTit2Den,title,nbins,QinvLo,QinvHi,200,0.0,6.0);
  49
  50   // to enable error bar calculation...
  51   fChi2ITSSUMNumerator->Sumw2();
  52   fChi2ITSSUMDenominator->Sumw2();
  53
  54   fChi2TPCSUMNumerator->Sumw2();
  55   fChi2TPCSUMDenominator->Sumw2();
  56   // set up numerator
  57   snprintf(tTitNum , 100, "%s%s","NumChi2ITSONE",title);
  58   fChi2ITSONENumerator = new TH2D(tTitNum,title,nbins,QinvLo,QinvHi,200,0.0,6.0);
  59   // set up denominator
  60   snprintf(tTitDen , 100, "%s%s", "DenChi2ITSONE", title);
  61   fChi2ITSONEDenominator = new TH2D(tTitDen,title,nbins,QinvLo,QinvHi,200,0.0,6.0);
  62
  63   // set up numerator
  64   snprintf(tTit2Num , 100, "%s%s","NumChi2TPCONE",title);
  65   fChi2TPCONENumerator = new TH2D(tTit2Num,title,nbins,QinvLo,QinvHi,200,0.0,6.0);
  66   // set up denominator
  67   snprintf(tTit2Den , 100, "%s%s", "DenChi2TPCONE", title);
  68   fChi2TPCONEDenominator = new TH2D(tTit2Den,title,nbins,QinvLo,QinvHi,200,0.0,6.0);
  69
  70   // set up numerator
  71   snprintf(tTit2Num , 100, "%s%s","NumSigmaToVertex",title);
  72   fSigmaToVertexNumerator = new TH2D(tTit2Num,title,nbins,QinvLo,QinvHi,200,0.0,6.0);
  73   // set up denominator
  74   snprintf(tTit2Den , 100, "%s%s", "DenSigmaToVertex", title);
  75   fSigmaToVertexDenominator = new TH2D(tTit2Den,title,nbins,QinvLo,QinvHi,200,0.0,6.0);
  76
  77   // to enable error bar calculation...
  78   fChi2ITSONENumerator->Sumw2();
  79   fChi2ITSONEDenominator->Sumw2();
  80
  81   fChi2TPCONENumerator->Sumw2();
  82   fChi2TPCONEDenominator->Sumw2();
  83 }
  84
  85 //____________________________
  86 AliFemtoChi2CorrFctn::AliFemtoChi2CorrFctn(const AliFemtoChi2CorrFctn& aCorrFctn) :
  87   AliFemtoCorrFctn(),
  88   fChi2ITSSUMNumerator(0),
  89   fChi2ITSSUMDenominator(0),
  90   fChi2TPCSUMNumerator(0),
  91   fChi2TPCSUMDenominator(0),
  92   fChi2ITSONENumerator(0),
  93   fChi2ITSONEDenominator(0),
  94   fChi2TPCONENumerator(0),
  95   fChi2TPCONEDenominator(0),
  96   fSigmaToVertexNumerator(0),
  97   fSigmaToVertexDenominator(0)
  98 {
  99   // copy constructor
 100   if (aCorrFctn.fChi2ITSSUMNumerator)
 101     fChi2ITSSUMNumerator = new TH2D(*aCorrFctn.fChi2ITSSUMNumerator);
 102   if (aCorrFctn.fChi2ITSSUMDenominator)
 103     fChi2ITSSUMDenominator = new TH2D(*aCorrFctn.fChi2ITSSUMDenominator);
 104   if (aCorrFctn.fChi2TPCSUMNumerator)
 105     fChi2TPCSUMNumerator = new TH2D(*aCorrFctn.fChi2TPCSUMNumerator);
 106   if (aCorrFctn.fChi2TPCSUMDenominator)
 107     fChi2TPCSUMDenominator = new TH2D(*aCorrFctn.fChi2TPCSUMDenominator);
 108   if (aCorrFctn.fChi2ITSONENumerator)
 109     fChi2ITSONENumerator = new TH2D(*aCorrFctn.fChi2ITSONENumerator);
 110   if (aCorrFctn.fChi2ITSONEDenominator)
 111     fChi2ITSONEDenominator = new TH2D(*aCorrFctn.fChi2ITSONEDenominator);
 112   if (aCorrFctn.fChi2TPCONENumerator)
 113     fChi2TPCONENumerator = new TH2D(*aCorrFctn.fChi2TPCONENumerator);
 114   if (aCorrFctn.fChi2TPCONEDenominator)
 115     fChi2TPCONEDenominator = new TH2D(*aCorrFctn.fChi2TPCONEDenominator);
 116   if (aCorrFctn.fSigmaToVertexNumerator)
 117     fSigmaToVertexNumerator = new TH2D(*aCorrFctn.fSigmaToVertexNumerator);
 118   if (aCorrFctn.fSigmaToVertexDenominator)
 119     fSigmaToVertexDenominator = new TH2D(*aCorrFctn.fSigmaToVertexDenominator);
 120 }
 121 //____________________________
 122 AliFemtoChi2CorrFctn::~AliFemtoChi2CorrFctn(){
 123   // destructor
 124   delete fChi2ITSSUMNumerator;
 125   delete fChi2ITSSUMDenominator;
 126   delete fChi2TPCSUMNumerator;
 127   delete fChi2TPCSUMDenominator;
 128   delete fChi2ITSONENumerator;
 129   delete fChi2ITSONEDenominator;
 130   delete fChi2TPCONENumerator;
 131   delete fChi2TPCONEDenominator;
 132   delete fSigmaToVertexNumerator;
 133   delete fSigmaToVertexDenominator;
 134 }
 135 //_________________________
 136 AliFemtoChi2CorrFctn& AliFemtoChi2CorrFctn::operator=(const AliFemtoChi2CorrFctn& aCorrFctn)
 137 {
 138   // assignment operator
 139   if (this == &aCorrFctn)
 140     return *this;
 141
 142   if (aCorrFctn.fChi2ITSSUMNumerator)
 143     fChi2ITSSUMNumerator = new TH2D(*aCorrFctn.fChi2ITSSUMNumerator);
 144   else
 145     fChi2ITSSUMNumerator = 0;
 146   if (aCorrFctn.fChi2ITSSUMDenominator)
 147     fChi2ITSSUMDenominator = new TH2D(*aCorrFctn.fChi2ITSSUMDenominator);
 148   else
 149     fChi2ITSSUMDenominator = 0;
 150   if (aCorrFctn.fChi2TPCSUMNumerator)
 151     fChi2TPCSUMNumerator = new TH2D(*aCorrFctn.fChi2TPCSUMNumerator);
 152   else
 153     fChi2TPCSUMNumerator = 0;
 154   if (aCorrFctn.fChi2TPCSUMDenominator)
 155     fChi2TPCSUMDenominator = new TH2D(*aCorrFctn.fChi2TPCSUMDenominator);
 156   else
 157     fChi2TPCSUMDenominator = 0;
 158   if (aCorrFctn.fChi2ITSONENumerator)
 159     fChi2ITSONENumerator = new TH2D(*aCorrFctn.fChi2ITSONENumerator);
 160   else
 161     fChi2ITSONENumerator = 0;
 162   if (aCorrFctn.fChi2ITSONEDenominator)
 163     fChi2ITSONEDenominator = new TH2D(*aCorrFctn.fChi2ITSONEDenominator);
 164   else
 165     fChi2ITSONEDenominator = 0;
 166   if (aCorrFctn.fChi2TPCONENumerator)
 167     fChi2TPCONENumerator = new TH2D(*aCorrFctn.fChi2TPCONENumerator);
 168   else
 169     fChi2TPCONENumerator = 0;
 170   if (aCorrFctn.fChi2TPCONEDenominator)
 171     fChi2TPCONEDenominator = new TH2D(*aCorrFctn.fChi2TPCONEDenominator);
 172   else
 173     fChi2TPCONEDenominator = 0;
 174   if (aCorrFctn.fSigmaToVertexNumerator)
 175     fSigmaToVertexNumerator = new TH2D(*aCorrFctn.fSigmaToVertexNumerator);
 176   else
 177     fSigmaToVertexNumerator = 0;
 178   if (aCorrFctn.fSigmaToVertexDenominator)
 179     fSigmaToVertexDenominator = new TH2D(*aCorrFctn.fSigmaToVertexDenominator);
 180   else
 181     fSigmaToVertexDenominator = 0;
 182
 183   return *this;
 184 }
 185 //_________________________
 186 void AliFemtoChi2CorrFctn::Finish(){
 187   // here is where we should normalize, fit, etc...
 188   // we should NOT Draw() the histos (as I had done it below),
 189   // since we want to insulate ourselves from root at this level
 190   // of the code.  Do it instead at root command line with browser.
 191   //  mShareNumerator->Draw();
 192   //mShareDenominator->Draw();
 193   //mRatio->Draw();
 194
 195 }
 196
 197 //____________________________
 198 AliFemtoString AliFemtoChi2CorrFctn::Report(){
 199   // create report
 200   string stemp = "ITS and TPC quality Correlation Function Report:\n";
 201   char ctemp[100];
 202   snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fChi2ITSSUMNumerator->GetEntries());
 203   stemp += ctemp;
 204   snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fChi2ITSSUMDenominator->GetEntries());
 205   stemp += ctemp;
 206   //  stemp += mCoulombWeight->Report();
 207   AliFemtoString returnThis = stemp;
 208   return returnThis;
 209 }
 210 //____________________________
 211 void AliFemtoChi2CorrFctn::AddRealPair( AliFemtoPair* pair){
 212   // add real (effect) pair
 213   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
 214
 215   if ((pair->Track1()->Track()->ITSncls() == 0) && (pair->Track2()->Track()->ITSncls() == 0))
 216     fChi2ITSSUMNumerator->Fill(tQinv, 1000.0);
 217   else
 218     fChi2ITSSUMNumerator->Fill(tQinv,
 219                                (pair->Track1()->Track()->ITSchi2() +
 220                                 pair->Track2()->Track()->ITSchi2())/
 221                                (pair->Track1()->Track()->ITSncls() +
 222                                 pair->Track2()->Track()->ITSncls()));
 223   if ((pair->Track1()->Track()->TPCncls() == 0) && (pair->Track2()->Track()->TPCncls() == 0))
 224     fChi2TPCSUMNumerator->Fill(tQinv, 1000.0);
 225   else
 226     fChi2TPCSUMNumerator->Fill(tQinv,
 227                                (pair->Track1()->Track()->TPCchi2() +
 228                                 pair->Track2()->Track()->TPCchi2())/
 229                                (pair->Track1()->Track()->TPCncls() +
 230                                 pair->Track2()->Track()->TPCncls()));
 231   double chi2perpointITS1, chi2perpointITS2;
 232   if (pair->Track1()->Track()->ITSncls() == 0)
 233     chi2perpointITS1 = 1000.0;
 234   else
 235     chi2perpointITS1 = pair->Track1()->Track()->ITSchi2()/pair->Track1()->Track()->ITSncls();
 236
 237   if (pair->Track2()->Track()->ITSncls() == 0)
 238     chi2perpointITS2 = 1000.0;
 239   else
 240     chi2perpointITS2 = pair->Track2()->Track()->ITSchi2()/pair->Track2()->Track()->ITSncls();
 241
 242
 243   if (chi2perpointITS1 > chi2perpointITS2) {
 244     fChi2ITSONENumerator->Fill(tQinv, chi2perpointITS1);
 245   }
 246   else {
 247     fChi2ITSONENumerator->Fill(tQinv, chi2perpointITS2);
 248   }
 249
 250   double chi2perpointTPC1, chi2perpointTPC2;
 251   if (pair->Track1()->Track()->TPCncls() == 0)
 252     chi2perpointTPC1 = 1000.0;
 253   else
 254     chi2perpointTPC1 = pair->Track1()->Track()->TPCchi2()/pair->Track1()->Track()->TPCncls();
 255
 256   if (pair->Track2()->Track()->TPCncls() == 0)
 257     chi2perpointTPC2 = 1000.0;
 258   else
 259     chi2perpointTPC2 = pair->Track2()->Track()->TPCchi2()/pair->Track2()->Track()->TPCncls();
 260
 261
 262   if (chi2perpointTPC1 > chi2perpointTPC2) {
 263     fChi2TPCONENumerator->Fill(tQinv, chi2perpointTPC1);
 264   }
 265   else {
 266     fChi2TPCONENumerator->Fill(tQinv, chi2perpointTPC2);
 267   }
 268
 269   if (pair->Track1()->Track()->SigmaToVertex() > pair->Track2()->Track()->SigmaToVertex()) {
 270     fSigmaToVertexNumerator->Fill(tQinv,
 271                                   pair->Track1()->Track()->SigmaToVertex());
 272   }
 273   else {
 274     fSigmaToVertexNumerator->Fill(tQinv,
 275                                   pair->Track2()->Track()->SigmaToVertex());
 276   }
 277 }
 278 //____________________________
 279 void AliFemtoChi2CorrFctn::AddMixedPair( AliFemtoPair* pair){
 280   // add mixed (background) pair
 281   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
 282
 283   if ((pair->Track1()->Track()->ITSncls() == 0) && (pair->Track2()->Track()->ITSncls() == 0))
 284     fChi2ITSSUMDenominator->Fill(tQinv, 1000.0);
 285   else
 286     fChi2ITSSUMDenominator->Fill(tQinv,
 287                                  (pair->Track1()->Track()->ITSchi2() +
 288                                   pair->Track2()->Track()->ITSchi2())/
 289                                  (pair->Track1()->Track()->ITSncls() +
 290                                   pair->Track2()->Track()->ITSncls()));
 291   if ((pair->Track1()->Track()->TPCncls() == 0) && (pair->Track2()->Track()->TPCncls() == 0))
 292     fChi2TPCSUMDenominator->Fill(tQinv, 1000.0);
 293   else
 294     fChi2TPCSUMDenominator->Fill(tQinv,
 295                                  (pair->Track1()->Track()->TPCchi2() +
 296                                   pair->Track2()->Track()->TPCchi2())/
 297                                  (pair->Track1()->Track()->TPCncls() +
 298                                   pair->Track2()->Track()->TPCncls()));
 299   double chi2perpointITS1, chi2perpointITS2;
 300   if (pair->Track1()->Track()->ITSncls() == 0)
 301     chi2perpointITS1 = 1000.0;
 302   else
 303     chi2perpointITS1 = pair->Track1()->Track()->ITSchi2()/pair->Track1()->Track()->ITSncls();
 304
 305   if (pair->Track2()->Track()->ITSncls() == 0)
 306     chi2perpointITS2 = 1000.0;
 307   else
 308     chi2perpointITS2 = pair->Track2()->Track()->ITSchi2()/pair->Track2()->Track()->ITSncls();
 309
 310
 311   if (chi2perpointITS1 > chi2perpointITS2) {
 312     fChi2ITSONEDenominator->Fill(tQinv, chi2perpointITS1);
 313   }
 314   else {
 315     fChi2ITSONEDenominator->Fill(tQinv, chi2perpointITS2);
 316   }
 317
 318   double chi2perpointTPC1, chi2perpointTPC2;
 319   if (pair->Track1()->Track()->TPCncls() == 0)
 320     chi2perpointTPC1 = 1000.0;
 321   else
 322     chi2perpointTPC1 = pair->Track1()->Track()->TPCchi2()/pair->Track1()->Track()->TPCncls();
 323
 324   if (pair->Track2()->Track()->TPCncls() == 0)
 325     chi2perpointTPC2 = 1000.0;
 326   else
 327     chi2perpointTPC2 = pair->Track2()->Track()->TPCchi2()/pair->Track2()->Track()->TPCncls();
 328
 329
 330   if (chi2perpointTPC1 > chi2perpointTPC2) {
 331     fChi2TPCONEDenominator->Fill(tQinv, chi2perpointTPC1);
 332   }
 333   else {
 334     fChi2TPCONEDenominator->Fill(tQinv, chi2perpointTPC2);
 335   }
 336   if (pair->Track1()->Track()->SigmaToVertex() > pair->Track2()->Track()->SigmaToVertex()) {
 337     fSigmaToVertexDenominator->Fill(tQinv,
 338                                   pair->Track1()->Track()->SigmaToVertex());
 339   }
 340   else {
 341     fSigmaToVertexDenominator->Fill(tQinv,
 342                                   pair->Track2()->Track()->SigmaToVertex());
 343   }
 344 }
 345
 346
 347 void AliFemtoChi2CorrFctn::WriteHistos()
 348 {
 349   // Write out result histograms
 350   fChi2ITSSUMNumerator->Write();
 351   fChi2ITSSUMDenominator->Write();
 352   fChi2TPCSUMNumerator->Write();
 353   fChi2TPCSUMDenominator->Write();
 354   fChi2ITSONENumerator->Write();
 355   fChi2ITSONEDenominator->Write();
 356   fChi2TPCONENumerator->Write();
 357   fChi2TPCONEDenominator->Write();
 358   fSigmaToVertexNumerator->Write();
 359   fSigmaToVertexDenominator->Write();
 360
 361 }
 362
 363 TList* AliFemtoChi2CorrFctn::GetOutputList()
 364 {
 365   // Prepare the list of objects to be written to the output
 366   TList *tOutputList = new TList();
 367
 368   tOutputList->Add(fChi2ITSSUMNumerator);
 369   tOutputList->Add(fChi2ITSSUMDenominator);
 370   tOutputList->Add(fChi2TPCSUMNumerator);
 371   tOutputList->Add(fChi2TPCSUMDenominator);
 372   tOutputList->Add(fChi2ITSONENumerator);
 373   tOutputList->Add(fChi2ITSONEDenominator);
 374   tOutputList->Add(fChi2TPCONENumerator);
 375   tOutputList->Add(fChi2TPCONEDenominator);
 376   tOutputList->Add(fSigmaToVertexNumerator);
 377   tOutputList->Add(fSigmaToVertexDenominator);
 378
 379   return tOutputList;
 380
 381 }