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