HBTAN/hbtanalysis.C

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