HBTAN/hbtanalysis.C

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