HBTAN/hbtanalysis.C

   1   #include "AliHBTAnalysis.h"
   2 //  #include "alles.h"
   3   #include "AliHBTReader.h"
   4   #include "AliHBTReaderKineTree.h"
   5   #include "AliHBTReaderITSv2.h"
   6   #include "AliHBTReaderITSv1.h"
   7   #include "AliHBTReaderTPC.h"
   8   #include "AliHBTReaderInternal.h"
   9   #include "AliHBTParticleCut.h"
  10   #include "AliHBTEvent.h"
  11   #include "AliHBTPairCut.h"
  12   #include "AliHBTQResolutionFctns.h"
  13   #include "AliHBTQDistributionFctns.h"
  14   #include "AliHBTTwoTrackEffFctn.h"
  15   #include "AliHBTCorrelFctn.h"
  16   #include "TSystem.h"
  17   #include "TObjString.h"
  18   #include "TString.h"
  19   #include "TPDGCode.h"
  20   #include "AliHBTMonDistributionFctns.h"
  21   #include "AliHBTMonResolutionFctns.h"
  22
  23 void hbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles",
  24                 Int_t first = -1,Int_t last = -1,
  25                 char *outfile = "hbtanalysis.root")
  26  {
  27 //HBT Anlysis Macro
  28 //Anlyzes TPC recontructed tracks and simulated particles that corresponds to them
  29
  30 //datatype defines type of data to be read
  31 //  Kine  - analyzes Kine Tree: simulated particles
  32 //  TPC   - analyzes TPC   tracking + particles corresponding to these tracks
  33 //  ITSv1 - analyzes ITSv1 ----------------------//--------------------------
  34 //  ITSv2 - analyzes ITSv2 ----------------------//--------------------------
  35
  36 //processopt defines option passed to AliHBTAnlysis::Process method
  37 // default: TracksAndParticles - process both recontructed tracks and sim. particles corresponding to them
  38 //          Tracks - process only recontructed tracks
  39 //          Particles - process only simulated particles
  40
  41 //Reads data from diroctories from first to last(including)
  42 // For examples if first=3 and last=5 it reads from
  43 //  ./3/
  44 //  ./4/
  45 //  ./5/
  46 //if first or last is negative (or both), it reads from current directory
  47 //
  48 //these names I use when analysis is done directly from CASTOR files via RFIO
  49 //  const char* basedir="rfio:/castor/cern.ch/user/s/skowron/";
  50 //  const char* serie="standard";
  51 //  const char* field = "0.4";
  52   const char* basedir=".";
  53   const char* serie="";
  54   const char* field = "";
  55
  56   Bool_t threeDcuts = kFALSE;
  57
  58   // Dynamically link some shared libs
  59 //  gROOT->LoadMacro("loadlibs.C");
  60 //  loadlibs();
  61
  62   gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libHBTAN");
  63
  64
  65   /***********************************************************/
  66   //Create analysis and reader
  67   AliHBTAnalysis * analysis = new AliHBTAnalysis();
  68
  69   AliHBTReader* reader;
  70   Int_t kine = strcmp(datatype,"Kine");
  71   Int_t ESD = strcmp(datatype,"ESD");
  72   Int_t TPC = strcmp(datatype,"TPC");
  73   Int_t ITSv1 = strcmp(datatype,"ITSv1");
  74   Int_t ITSv2 = strcmp(datatype,"ITSv2");
  75   Int_t intern = strcmp(datatype,"Intern");
  76
  77   if(!kine)
  78    {
  79     reader = new AliHBTReaderKineTree();
  80     processopt="Particles"; //this reader by definition reads only simulated particles
  81    }
  82   else if(!ESD)
  83    {
  84     AliHBTReaderESD* esdreader = new AliHBTReaderESD();
  85     esdreader->ReadParticles(kTRUE);
  86     reader = esdreader;
  87    }
  88   else if(!TPC)
  89    {
  90     reader = new AliHBTReaderTPC();
  91    }
  92   else if(!ITSv1)
  93    {
  94     reader = new AliHBTReaderITSv1();
  95    }
  96   else if(!ITSv2)
  97    {
  98     reader = new AliHBTReaderITSv2();
  99    }
 100   else if(!intern)
 101    {
 102     reader = new AliHBTReaderInternal("data.root");
 103    }
 104   else
 105    {
 106     cerr<<"Option "<<datatype<<"  not recognized. Exiting"<<endl;
 107     return;
 108    }
 109
 110   TObjArray* dirs=0;
 111   if ( (first >= 0) && (last>=0) && ( (last-first)>=0 ) )
 112    {//read from many dirs dirs
 113      char buff[50];
 114      dirs = new TObjArray(last-first+1);
 115      for (Int_t i = first; i<=last; i++)
 116       {
 117         sprintf(buff,"%s/%s/%s/%d",basedir,field,serie,i);
 118         TObjString *odir= new TObjString(buff);
 119         dirs->Add(odir);
 120       }
 121     }
 122    reader->SetDirs(dirs);
 123
 124   /***********************************************************/
 125   /*****   R E A D E R ***************************************/
 126   /***********************************************************/
 127
 128   //we want have only low pt pi+ so set a cut to reader
 129   AliHBTParticleCut* readerpartcut= new AliHBTParticleCut();
 130   readerpartcut->SetPtRange(0.0,10.0);
 131   readerpartcut->SetPID(kPiPlus);
 132   reader->AddParticleCut(readerpartcut);//read this particle type with this cut
 133
 134 //  readerpartcut->SetPtRange(0.0,10.0);
 135 //  readerpartcut->SetPID(kPiMinus);
 136 //  reader->AddParticleCut(readerpartcut);//read this particle type with this cut
 137
 138   analysis->SetReader(reader);
 139
 140   /***********************************************************/
 141   /*****   W E I G H T S        ******************************/
 142   /***********************************************************/
 143
 144   AliHBTLLWeights::Instance()->SetParticlesTypes(kPiPlus,kPiPlus);
 145   AliHBTLLWeights::Instance()->SetColoumb(kFALSE);
 146   AliHBTLLWeights::Instance()->SetStrongInterSwitch(kFALSE);
 147
 148   AliHBTLLWeights::Instance()->SetRandomPosition(kFALSE);
 149 //AliHBTLLWeights::Instance()->SetR1dw(8.0);
 150
 151   AliHBTLLWeights::Instance()->Init();
 152   AliHBTLLWeights::Instance()->Set();
 153
 154   //example function
 155   AliHBTWeightTheorQInvFctn  *wqinvcfP = new AliHBTWeightTheorQInvFctn(nbins,qmax);
 156   wqinvcfP->SetNumberOfBinsToScale(binstoscale);
 157   wqinvcfP->Rename("wqinvcfP","Lednicky Q_{inv} Theoretical Correlation Function "+system);
 158   analysis->AddParticleFunction(wqinvcfP);
 159
 160   /************************************************************/
 161   /****   Q INV Correlation Function   ************************/
 162   /************************************************************/
 163
 164   AliHBTPairCut *paircut = new AliHBTPairCut();
 165   Float_t qinvmin = 0.0;
 166   Float_t qinvmax = 0.05;//50MeV
 167   paircut->SetQInvRange(qinvmin,qinvmax);
 168
 169 //  AliHBTParticleCut* partcut= new AliHBTParticleCut();
 170 //  partcut->SetPID(kPiPlus);
 171 //  paircut->SetFirstPartCut(partcut);
 172 //  partcut->SetPID(kPiMinus);
 173 //  paircut->SetSecondPartCut(partcut);
 174
 175 //  analysis->SetGlobalPairCut(paircut);
 176
 177   AliHBTQInvCorrelFctn * qinvcfT = new AliHBTQInvCorrelFctn(100,qinvmax);
 178   AliHBTQInvCorrelFctn * qinvcfP = new AliHBTQInvCorrelFctn(100,qinvmax);
 179
 180   analysis->AddTrackFunction(qinvcfT);
 181   analysis->AddParticleFunction(qinvcfP);
 182
 183
 184   qinvcfP->Rename("qinvcfP","Particle (simulated) Qinv CF \\pi^{+} \\pi^{+}");
 185   qinvcfT->Rename("qinvcfT","Track (recontructed) Qinv CF \\pi^{+} \\pi^{+}");
 186
 187   /************************************************************/
 188   /****   Q OUT Correlation Function   ************************/
 189   /************************************************************/
 190
 191   AliHBTQOutCMSLCCorrelFctn* qoutP = new AliHBTQOutCMSLCCorrelFctn();
 192   qoutP->Rename("qoutP","Particle (simulated) Q_{out} CF \\pi^{+} \\pi^{+}");
 193   AliHBTQOutCMSLCCorrelFctn* qoutT = new AliHBTQOutCMSLCCorrelFctn();
 194   qoutT->Rename("qoutT","Track (recontructed) Q_{out} CF \\pi^{+} \\pi^{+}");
 195
 196   AliHBTPairCut *outPairCut = new AliHBTPairCut();
 197   outPairCut->SetQOutCMSLRange(0.0,0.15);
 198   outPairCut->SetQSideCMSLRange(0.0,0.02);
 199   outPairCut->SetQLongCMSLRange(0.0,0.02);
 200   qoutP->SetPairCut(outPairCut);
 201   qoutT->SetPairCut(outPairCut);
 202
 203   //analysis->AddTrackFunction(qoutP);
 204   //analysis->AddParticleFunction(qoutT);
 205
 206   /************************************************************/
 207   /****   Q SIDE Correlation Function   ***********************/
 208   /************************************************************/
 209
 210   AliHBTQSideCMSLCCorrelFctn* qsideP = new AliHBTQSideCMSLCCorrelFctn(100,qinvmax);
 211   qsideP->Rename("qsideP","Particle (simulated) Q_{side} CF \\pi^{+} \\pi^{+}");
 212   AliHBTQSideCMSLCCorrelFctn* qsideT = new AliHBTQSideCMSLCCorrelFctn(100,qinvmax);
 213   qsideT->Rename("qsideT","Track (recontructed) Q_{side} CF \\pi^{+} \\pi^{+}");
 214
 215   AliHBTPairCut *sidePairCut = new AliHBTPairCut();
 216   sidePairCut->SetQOutCMSLRange(0.0,0.02);
 217   sidePairCut->SetQSideCMSLRange(0.0,0.15);
 218   sidePairCut->SetQLongCMSLRange(0.0,0.02);
 219   qsideP->SetPairCut(sidePairCut);
 220   qsideT->SetPairCut(sidePairCut);
 221
 222   //analysis->AddTrackFunction(qsideP);
 223   //analysis->AddParticleFunction(qsideT);
 224
 225   /************************************************************/
 226   /****   Q LONG Correlation Function   ***********************/
 227   /************************************************************/
 228
 229   AliHBTQLongCMSLCCorrelFctn* qlongP = new AliHBTQLongCMSLCCorrelFctn(100,qinvmax);
 230   qlongP->Rename("qlongP","Particle (simulated) Q_{long} CF \\pi^{+} \\pi^{+}");
 231   AliHBTQLongCMSLCCorrelFctn* qlongT = new AliHBTQLongCMSLCCorrelFctn(100,qinvmax);
 232   qlongT->Rename("qlongT","Track (recontructed) Q_{long} CF \\pi^{+} \\pi^{+}");
 233
 234   AliHBTPairCut *longPairCut = new AliHBTPairCut();
 235   longPairCut->SetQOutCMSLRange(0.0,0.02);
 236   longPairCut->SetQSideCMSLRange(0.0,0.02);
 237   longPairCut->SetQLongCMSLRange(0.0,0.15);
 238   qlongP->SetPairCut(longPairCut);
 239   qlongT->SetPairCut(longPairCut);
 240
 241   //analysis->AddTrackFunction(qlongT);
 242   //analysis->AddParticleFunction(qlongT);
 243
 244
 245   /************************************************************/
 246   /************************************************************/
 247   /************************************************************/
 248   /****   R E S O L U T I O N S         ***********************/
 249   /************************************************************/
 250   /************************************************************/
 251   /************************************************************/
 252
 253
 254
 255   AliHBTQInvResolVsKtFctn qinvVsktF(200,1.0,0.0,300,0.015,-0.015);    //QInvCMSLC  Res   Vs   Kt
 256   if (threeDcuts)  qinvVsktF.SetPairCut(paircut);
 257   //analysis->AddResolutionFunction(&qinvVsktF);
 258
 259   AliHBTQOutResolVsKtFctn qoutVsktF(200,1.0,0.0,300,0.05,-0.05);    //QOutCMSLC  Res   Vs   Kt
 260   if (threeDcuts)  qoutVsktF.SetPairCut(outPairCut);
 261   //analysis->AddResolutionFunction(&qoutVsktF);
 262
 263   AliHBTQSideResolVsKtFctn qsideVsktF(200,1.0,0.0,300,0.015,-0.015);   //QSideCMSLC Res   Vs   Kt
 264   if (threeDcuts)  qsideVsktF.SetPairCut(sidePairCut);
 265   //analysis->AddResolutionFunction(&qsideVsktF);
 266
 267   AliHBTQLongResolVsKtFctn qlongVsktF(200,1.0,0.0,300,0.015,-0.015);   //QLongCMSLC Res   Vs   Kt
 268   if (threeDcuts)  qlongVsktF.SetPairCut(longPairCut);
 269   //analysis->AddResolutionFunction(&qlongVsktF);
 270
 271   AliHBTQInvResolVsQInvFctn qInvVsqinvF(200,qinvmax,qinvmin,300,0.015,-0.015); //QInvCMSLC Res   Vs   QInvCMSLC
 272   //analysis->AddResolutionFunction(&qInvVsqinvF);
 273
 274   AliHBTQOutResolVsQInvFctn qoutVsqinvF(200,qinvmax,qinvmin,300,0.05,-0.05);  //QOutCMSLC  Res   Vs   QInvCMSLC
 275   if (threeDcuts)  qoutVsqinvF.SetPairCut(outPairCut);
 276   //analysis->AddResolutionFunction(&qoutVsqinvF);
 277
 278   AliHBTQSideResolVsQInvFctn qsideVsqinvF(200,qinvmax,qinvmin,300,0.015,-0.015); //QSideCMSLC Res   Vs   QInvCMSLC
 279   if (threeDcuts)  qsideVsqinvF.SetPairCut(sidePairCut);
 280   //analysis->AddResolutionFunction(&qsideVsqinvF);
 281
 282   AliHBTQLongResolVsQInvFctn qlongVsqinvF(200,qinvmax,qinvmin,300,0.015,-0.015); //QLongCMSLC Res   Vs   QInvCMSLC
 283   if (threeDcuts)  qlongVsqinvF.SetPairCut(longPairCut);
 284   //analysis->AddResolutionFunction(&qlongVsqinvF);
 285
 286
 287   AliHBTPairThetaResolVsQInvFctn pairThetaVsqinv(200,qinvmax,qinvmin,300,0.015,-0.015);
 288   if (threeDcuts) pairThetaVsqinv.SetPairCut(longPairCut);
 289   //analysis->AddResolutionFunction(&pairThetaVsqinv);
 290
 291   AliHBTPairThetaResolVsKtFctn pairThetaVsKt(200,1.0,0.0,300,0.015,-0.015);
 292   if (threeDcuts) pairThetaVsKt.SetPairCut(longPairCut);
 293   //analysis->AddResolutionFunction(&pairThetaVsKt);
 294
 295   AliHBTPairCut phipc;
 296   phipc.SetQLongCMSLRange(0.0,0.02);
 297
 298   AliHBTPairPhiResolVsQInvFctn pairPhiVsqinv(200,qinvmax,qinvmin,300,0.015,-0.015);
 299   if (threeDcuts) pairPhiVsqinv.SetPairCut(&phipc);
 300   //analysis->AddResolutionFunction(&pairPhiVsqinv);
 301
 302   AliHBTPairPhiResolVsKtFctn pairPhiVsKt(200,1.0,0.0,300,0.015,-0.015);
 303   if (threeDcuts) pairPhiVsKt.SetPairCut(&phipc);
 304   //analysis->AddResolutionFunction(&pairPhiVsKt);
 305
 306
 307   AliHBTQOutResolVsQOutFctn qoutVsqoutF;  //QOutCMSLC  Res   Vs   QOut
 308   if (threeDcuts) qoutVsqoutF.SetPairCut(outPairCut);
 309   //analysis->AddResolutionFunction(&qoutVsqoutF);
 310
 311   AliHBTQSideResolVsQSideFctn qsideVsqsideF;//QSideCMSLC Res   Vs   QSide
 312   if (threeDcuts) qsideVsqsideF.SetPairCut(sidePairCut);
 313   //analysis->AddResolutionFunction(&qsideVsqsideF);
 314
 315
 316   AliHBTQLongResolVsQLongFctn qlongVsqlongF;//QLongCMSLC Res   Vs   QLong
 317   if (threeDcuts) qlongVsqlongF.SetPairCut(longPairCut);
 318   //analysis->AddResolutionFunction(&qlongVsqlongF);
 319
 320
 321   /************************************************************/
 322   /************************************************************/
 323   /************************************************************/
 324   /****   D I S T R B U T I O N S         ***********************/
 325   /************************************************************/
 326   /************************************************************/
 327   /************************************************************/
 328
 329   /************************************************************/
 330   /****   OUT VS KT                 ***********************/
 331   /************************************************************/
 332   AliHBTQOutDistributionVsKtFctn qoutVsKtDistP;
 333   qoutVsKtDistP.Rename("qoutVsKtDistP",qoutVsKtDistP.GetTitle());
 334   AliHBTQOutDistributionVsKtFctn qoutVsKtDistT;
 335   qoutVsKtDistT.Rename("qoutVsKtDistT",qoutVsKtDistT.GetTitle());
 336   if (threeDcuts)  qoutVsKtDistP.SetPairCut(outPairCut);
 337   if (threeDcuts)  qoutVsKtDistT.SetPairCut(outPairCut);
 338   //analysis->AddParticleFunction(&qoutVsKtDistP);
 339   //analysis->AddTrackFunction(&qoutVsKtDistT);
 340   /************************************************************/
 341   /****   Side VS KT                 ***********************/
 342   /************************************************************/
 343   AliHBTQSideDistributionVsKtFctn qSideVsKtDistP;
 344   qSideVsKtDistP.Rename("qSideVsKtDistP",qSideVsKtDistP.GetTitle());
 345   AliHBTQSideDistributionVsKtFctn qSideVsKtDistT;
 346   qSideVsKtDistT.Rename("qSideVsKtDistT",qSideVsKtDistT.GetTitle());
 347   if (threeDcuts)  qSideVsKtDistP.SetPairCut(sidePairCut);
 348   if (threeDcuts)  qSideVsKtDistT.SetPairCut(sidePairCut);
 349   //analysis->AddParticleFunction(&qSideVsKtDistP);
 350   //analysis->AddTrackFunction(&qSideVsKtDistT);
 351   /************************************************************/
 352   /****   Long VS KT                 ***********************/
 353   /************************************************************/
 354   AliHBTQLongDistributionVsKtFctn qLongVsKtDistP;
 355   qLongVsKtDistP.Rename("qLongVsKtDistP",qLongVsKtDistP.GetTitle());
 356   AliHBTQLongDistributionVsKtFctn qLongVsKtDistT;
 357   qLongVsKtDistT.Rename("qLongVsKtDistT",qLongVsKtDistT.GetTitle());
 358   if (threeDcuts)  qLongVsKtDistP.SetPairCut(longPairCut);
 359   if (threeDcuts)  qLongVsKtDistT.SetPairCut(longPairCut);
 360   //analysis->AddParticleFunction(&qLongVsKtDistP);
 361   //analysis->AddTrackFunction(&qLongVsKtDistT);
 362
 363   /************************************************************/
 364   /************************************************************/
 365   /************************************************************/
 366   /****    M O N I T O R  R E S O L U T I O N S   *************/
 367   /************************************************************/
 368   /************************************************************/
 369   /************************************************************/
 370
 371   AliHBTMonPxResolutionVsPtFctn PxVsPtResF(200,1.0,0.0,300,0.05,-0.05);
 372   //analysis->AddParticleAndTrackMonitorFunction(&PxVsPtResF);
 373
 374   AliHBTMonPyResolutionVsPtFctn PyVsPtResF(200,1.0,0.0,300,0.05,-0.05);
 375   //analysis->AddParticleAndTrackMonitorFunction(&PyVsPtResF);
 376
 377   AliHBTMonPzResolutionVsPtFctn PzVsPtResF(200,1.0,0.0,300,0.05,-0.05);
 378   //analysis->AddParticleAndTrackMonitorFunction(&PzVsPtResF);
 379
 380   AliHBTMonPResolutionVsPtFctn PVsPtResF(200,1.0,0.0,300,0.05,-0.05);
 381   //analysis->AddParticleAndTrackMonitorFunction(&PVsPtResF);
 382
 383   AliHBTMonPtResolutionVsPtFctn PtVsPtResF(400,4.0,0.0,300,0.07,-0.07);
 384   analysis->AddParticleAndTrackMonitorFunction(&PtVsPtResF);
 385
 386   AliHBTMonPhiResolutionVsPtFctn PhiVsPtResF(200,1.0,0.0,300,0.02,-0.02);
 387   analysis->AddParticleAndTrackMonitorFunction(&PhiVsPtResF);
 388
 389   AliHBTMonThetaResolutionVsPtFctn ThetaVsPtResF(200,1.0,0.0,300,0.02,-0.02);
 390   analysis->AddParticleAndTrackMonitorFunction(&ThetaVsPtResF);
 391
 392
 393
 394
 395   /************************************************************/
 396   /****   P R O C E S S                 ***********************/
 397   /************************************************************/
 398   analysis->SetBufferSize(2);
 399   analysis->Process(processopt);
 400
 401   TFile histoOutput(outfile,"recreate");
 402   analysis->WriteFunctions();
 403   histoOutput.Close();
 404
 405   delete qinvcfP;
 406   delete qinvcfT;
 407   delete paircut;
 408   delete readerpartcut;
 409   if (dirs)
 410    {
 411     dirs->SetOwner();
 412     delete dirs;
 413    }
 414   delete reader;
 415   delete analysis;
 416
 417  }
 418