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