]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HBTAN/hbtanalysis.C
const X& X::operator=(const X&) --> X& X::operator=(const X&)
[u/mrichter/AliRoot.git] / HBTAN / hbtanalysis.C
index 15e1eb59e8aee5f30608211d8a54b15a7c3b6b86..7fcc2e439172a5972cee6cdc4d3cab8e8b2ff955 100644 (file)
@@ -1,7 +1,33 @@
+  #include "AliHBTAnalysis.h"
+//  #include "alles.h"
+  #include "AliHBTReader.h"
+  #include "AliHBTReaderKineTree.h"
+  #include "AliHBTReaderITSv2.h"
+  #include "AliHBTReaderITSv1.h"
+  #include "AliHBTReaderTPC.h"
+  #include "AliHBTReaderInternal.h"
+  #include "AliHBTParticleCut.h"
+  #include "AliHBTEvent.h"
+  #include "AliHBTPairCut.h"
+  #include "AliHBTQResolutionFctns.h"
+  #include "AliHBTQDistributionFctns.h"
+  #include "AliHBTTwoTrackEffFctn.h"
+  #include "AliHBTCorrelFctn.h"
+  #include "TSystem.h"
+  #include "TObjString.h"
+  #include "TString.h"
+  #include "TPDGCode.h"
+  #include "AliHBTMonDistributionFctns.h"
+  #include "AliHBTMonResolutionFctns.h"
+  
 void hbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles", 
                 Int_t first = -1,Int_t last = -1, 
                 char *outfile = "hbtanalysis.root")
  {
+   
+//  AliHBTTrackPoints::SetDebug(1);
+  AliHBTParticle::SetDebug(0);
+  
 //HBT Anlysis Macro
 //Anlyzes TPC recontructed tracks and simulated particles that corresponds to them
 
@@ -27,13 +53,16 @@ void hbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles",
 //  const char* basedir="rfio:/castor/cern.ch/user/s/skowron/";
 //  const char* serie="standard";
 //  const char* field = "0.4";
+
   const char* basedir=".";
   const char* serie="";
   const char* field = "";
  
+  Bool_t threeDcuts = kFALSE;
+
   // Dynamically link some shared libs                    
-  gROOT->LoadMacro("loadlibs.C");
-  loadlibs();
+//  gROOT->LoadMacro("loadlibs.C");
+//  loadlibs();
   
   gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libHBTAN");
   
@@ -41,21 +70,32 @@ void hbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles",
   /***********************************************************/
   //Create analysis and reader
   AliHBTAnalysis * analysis = new AliHBTAnalysis();
+  analysis->SetCutsOnTracks();
   
   AliHBTReader* reader;
   Int_t kine = strcmp(datatype,"Kine");
+  Int_t ESD = strcmp(datatype,"ESD");
   Int_t TPC = strcmp(datatype,"TPC");
   Int_t ITSv1 = strcmp(datatype,"ITSv1");
   Int_t ITSv2 = strcmp(datatype,"ITSv2");
-
+  Int_t intern = strcmp(datatype,"Intern");
+  
   if(!kine)
    {
     reader = new AliHBTReaderKineTree();
     processopt="Particles"; //this reader by definition reads only simulated particles
    }
+  else if(!ESD)
+   {
+    AliHBTReaderESD* esdreader = new AliHBTReaderESD();
+    esdreader->ReadParticles(kTRUE);
+    reader = esdreader;
+   }
   else if(!TPC)
    {
     reader = new AliHBTReaderTPC();
+    //((AliHBTReaderTPC*)reader)->SetNumberOfTrackPoints(5,30.);
+    ((AliHBTReaderTPC*)reader)->SetClusterMap();
    }
   else if(!ITSv1)
    {
@@ -65,6 +105,10 @@ void hbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles",
    {
     reader = new AliHBTReaderITSv2();
    }
+  else if(!intern)
+   {
+    reader = new AliHBTReaderInternal("data.root");
+   }
   else
    {
     cerr<<"Option "<<datatype<<"  not recognized. Exiting"<<endl;
@@ -78,7 +122,6 @@ void hbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles",
      dirs = new TObjArray(last-first+1);
      for (Int_t i = first; i<=last; i++)
       {
-       if(i == 24) continue;
         sprintf(buff,"%s/%s/%s/%d",basedir,field,serie,i);
         TObjString *odir= new TObjString(buff);
         dirs->Add(odir);
@@ -87,31 +130,286 @@ void hbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles",
    reader->SetDirs(dirs);
    
   /***********************************************************/    
+  /*****   R E A D E R ***************************************/    
+  /***********************************************************/    
   
   //we want have only low pt pi+ so set a cut to reader
   AliHBTParticleCut* readerpartcut= new AliHBTParticleCut();
-  readerpartcut->SetPtRange(0.0,1.5);
+  readerpartcut->SetPtRange(0.0,10.0);
   readerpartcut->SetPID(kPiPlus);
   reader->AddParticleCut(readerpartcut);//read this particle type with this cut
   
-  analysis->SetReader(reader);
-  /************************************************************/
+//  readerpartcut->SetPtRange(0.0,10.0);
+//  readerpartcut->SetPID(kPiMinus);
+//  reader->AddParticleCut(readerpartcut);//read this particle type with this cut
   
+  analysis->SetReader(reader);
+
   AliHBTPairCut *paircut = new AliHBTPairCut();
-  paircut->SetQInvRange(0.0,0.20);  
+  Float_t qinvmin = 0.0;
+  Float_t qinvmax = 0.05;//50MeV
+  paircut->SetQInvRange(qinvmin,qinvmax);  
+  
+ // paircut->SetAvSeparationRange(10.); //AntiMerging
+  paircut->SetClusterOverlapRange(-.5,1.0);//AntiSplitting
+  
+//  AliHBTParticleCut* partcut= new AliHBTParticleCut();
+//  partcut->SetPID(kPiPlus);
+//  paircut->SetFirstPartCut(partcut);
+//  partcut->SetPID(kPiMinus);
+//  paircut->SetSecondPartCut(partcut);
+  
   analysis->SetGlobalPairCut(paircut);
+  
+  /***********************************************************/    
+  /*****   W E I G H T S        ******************************/    
+  /***********************************************************/    
 
-  AliHBTQInvCorrelFctn * qinvcfT= new AliHBTQInvCorrelFctn();
-  AliHBTQInvCorrelFctn * qinvcfP= new AliHBTQInvCorrelFctn();
+  AliHBTLLWeights::Instance()->SetParticlesTypes(kPiPlus,kPiPlus);
+  AliHBTLLWeights::Instance()->SetColoumb(kFALSE);
+  AliHBTLLWeights::Instance()->SetStrongInterSwitch(kFALSE);
+
+  AliHBTLLWeights::Instance()->SetRandomPosition(kFALSE);
+//AliHBTLLWeights::Instance()->SetR1dw(8.0);
+
+  AliHBTLLWeights::Instance()->Init();
+  AliHBTLLWeights::Instance()->Set();
+
+  //example function
+  AliHBTWeightTheorQInvFctn  *wqinvcfP = new AliHBTWeightTheorQInvFctn(100,qinvmax);
+  wqinvcfP->SetNumberOfBinsToScale(30);
+  wqinvcfP->Rename("wqinvcfP","Lednicky Q_{inv} Theoretical Correlation Function");
+  analysis->AddParticleFunction(wqinvcfP);
+
+  /************************************************************/
+  /****   Q INV Correlation Function   ************************/
+  /************************************************************/
+   
+
+  AliHBTQInvCorrelFctn * qinvcfT = new AliHBTQInvCorrelFctn(100,qinvmax);
+  AliHBTQInvCorrelFctn * qinvcfP = new AliHBTQInvCorrelFctn(100,qinvmax);
   
   analysis->AddTrackFunction(qinvcfT);
   analysis->AddParticleFunction(qinvcfP);
   
-  analysis->Process(processopt);
   
   qinvcfP->Rename("qinvcfP","Particle (simulated) Qinv CF \\pi^{+} \\pi^{+}");
   qinvcfT->Rename("qinvcfT","Track (recontructed) Qinv CF \\pi^{+} \\pi^{+}");
+  
+  /************************************************************/
+  /****   Q OUT Correlation Function   ************************/
+  /************************************************************/
+  
+  AliHBTQOutCMSLCCorrelFctn* qoutP = new AliHBTQOutCMSLCCorrelFctn();
+  qoutP->Rename("qoutP","Particle (simulated) Q_{out} CF \\pi^{+} \\pi^{+}");
+  AliHBTQOutCMSLCCorrelFctn* qoutT = new AliHBTQOutCMSLCCorrelFctn(); 
+  qoutT->Rename("qoutT","Track (recontructed) Q_{out} CF \\pi^{+} \\pi^{+}");
+
+  AliHBTPairCut *outPairCut = new AliHBTPairCut();
+  outPairCut->SetQOutCMSLRange(0.0,0.15);
+  outPairCut->SetQSideCMSLRange(0.0,0.02);
+  outPairCut->SetQLongCMSLRange(0.0,0.02);
+  qoutP->SetPairCut(outPairCut);
+  qoutT->SetPairCut(outPairCut);
+
+  //analysis->AddTrackFunction(qoutP);
+  //analysis->AddParticleFunction(qoutT);
+
+  /************************************************************/
+  /****   Q SIDE Correlation Function   ***********************/
+  /************************************************************/
+  
+  AliHBTQSideCMSLCCorrelFctn* qsideP = new AliHBTQSideCMSLCCorrelFctn(100,qinvmax); 
+  qsideP->Rename("qsideP","Particle (simulated) Q_{side} CF \\pi^{+} \\pi^{+}");
+  AliHBTQSideCMSLCCorrelFctn* qsideT = new AliHBTQSideCMSLCCorrelFctn(100,qinvmax); 
+  qsideT->Rename("qsideT","Track (recontructed) Q_{side} CF \\pi^{+} \\pi^{+}");
+
+  AliHBTPairCut *sidePairCut = new AliHBTPairCut();
+  sidePairCut->SetQOutCMSLRange(0.0,0.02);
+  sidePairCut->SetQSideCMSLRange(0.0,0.15);
+  sidePairCut->SetQLongCMSLRange(0.0,0.02);
+  qsideP->SetPairCut(sidePairCut);
+  qsideT->SetPairCut(sidePairCut);
+
+  //analysis->AddTrackFunction(qsideP);
+  //analysis->AddParticleFunction(qsideT);
+
+  /************************************************************/
+  /****   Q LONG Correlation Function   ***********************/
+  /************************************************************/
+    
+  AliHBTQLongCMSLCCorrelFctn* qlongP = new AliHBTQLongCMSLCCorrelFctn(100,qinvmax);
+  qlongP->Rename("qlongP","Particle (simulated) Q_{long} CF \\pi^{+} \\pi^{+}");
+  AliHBTQLongCMSLCCorrelFctn* qlongT = new AliHBTQLongCMSLCCorrelFctn(100,qinvmax); 
+  qlongT->Rename("qlongT","Track (recontructed) Q_{long} CF \\pi^{+} \\pi^{+}");
+
+  AliHBTPairCut *longPairCut = new AliHBTPairCut();
+  longPairCut->SetQOutCMSLRange(0.0,0.02);
+  longPairCut->SetQSideCMSLRange(0.0,0.02);
+  longPairCut->SetQLongCMSLRange(0.0,0.15);
+  qlongP->SetPairCut(longPairCut);
+  qlongT->SetPairCut(longPairCut);
+
+  //analysis->AddTrackFunction(qlongT);
+  //analysis->AddParticleFunction(qlongT);
+
+
+  /************************************************************/
+  /************************************************************/
+  /************************************************************/
+  /****   R E S O L U T I O N S         ***********************/
+  /************************************************************/
+  /************************************************************/
+  /************************************************************/
+
 
+
+  AliHBTQInvResolVsKtFctn qinvVsktF(200,1.0,0.0,300,0.015,-0.015);    //QInvCMSLC  Res   Vs   Kt
+  if (threeDcuts)  qinvVsktF.SetPairCut(paircut);
+  //analysis->AddResolutionFunction(&qinvVsktF);
+
+  AliHBTQOutResolVsKtFctn qoutVsktF(200,1.0,0.0,300,0.05,-0.05);    //QOutCMSLC  Res   Vs   Kt
+  if (threeDcuts)  qoutVsktF.SetPairCut(outPairCut);
+  //analysis->AddResolutionFunction(&qoutVsktF);
+
+  AliHBTQSideResolVsKtFctn qsideVsktF(200,1.0,0.0,300,0.015,-0.015);   //QSideCMSLC Res   Vs   Kt
+  if (threeDcuts)  qsideVsktF.SetPairCut(sidePairCut);
+  //analysis->AddResolutionFunction(&qsideVsktF);
+
+  AliHBTQLongResolVsKtFctn qlongVsktF(200,1.0,0.0,300,0.015,-0.015);   //QLongCMSLC Res   Vs   Kt
+  if (threeDcuts)  qlongVsktF.SetPairCut(longPairCut);
+  //analysis->AddResolutionFunction(&qlongVsktF);
+
+  AliHBTQInvResolVsQInvFctn qInvVsqinvF(200,qinvmax,qinvmin,300,0.015,-0.015); //QInvCMSLC Res   Vs   QInvCMSLC
+  //analysis->AddResolutionFunction(&qInvVsqinvF);
+
+  AliHBTQOutResolVsQInvFctn qoutVsqinvF(200,qinvmax,qinvmin,300,0.05,-0.05);  //QOutCMSLC  Res   Vs   QInvCMSLC
+  if (threeDcuts)  qoutVsqinvF.SetPairCut(outPairCut);
+  //analysis->AddResolutionFunction(&qoutVsqinvF);
+
+  AliHBTQSideResolVsQInvFctn qsideVsqinvF(200,qinvmax,qinvmin,300,0.015,-0.015); //QSideCMSLC Res   Vs   QInvCMSLC
+  if (threeDcuts)  qsideVsqinvF.SetPairCut(sidePairCut);
+  //analysis->AddResolutionFunction(&qsideVsqinvF);
+
+  AliHBTQLongResolVsQInvFctn qlongVsqinvF(200,qinvmax,qinvmin,300,0.015,-0.015); //QLongCMSLC Res   Vs   QInvCMSLC
+  if (threeDcuts)  qlongVsqinvF.SetPairCut(longPairCut);
+  //analysis->AddResolutionFunction(&qlongVsqinvF);
+
+
+  AliHBTPairThetaResolVsQInvFctn pairThetaVsqinv(200,qinvmax,qinvmin,300,0.015,-0.015);
+  if (threeDcuts) pairThetaVsqinv.SetPairCut(longPairCut);
+  //analysis->AddResolutionFunction(&pairThetaVsqinv);
+
+  AliHBTPairThetaResolVsKtFctn pairThetaVsKt(200,1.0,0.0,300,0.015,-0.015);
+  if (threeDcuts) pairThetaVsKt.SetPairCut(longPairCut);
+  //analysis->AddResolutionFunction(&pairThetaVsKt);
+
+  AliHBTPairCut phipc;
+  phipc.SetQLongCMSLRange(0.0,0.02);
+
+  AliHBTPairPhiResolVsQInvFctn pairPhiVsqinv(200,qinvmax,qinvmin,300,0.015,-0.015);
+  if (threeDcuts) pairPhiVsqinv.SetPairCut(&phipc);
+  //analysis->AddResolutionFunction(&pairPhiVsqinv);
+
+  AliHBTPairPhiResolVsKtFctn pairPhiVsKt(200,1.0,0.0,300,0.015,-0.015);
+  if (threeDcuts) pairPhiVsKt.SetPairCut(&phipc);
+  //analysis->AddResolutionFunction(&pairPhiVsKt);
+
+
+  AliHBTQOutResolVsQOutFctn qoutVsqoutF;  //QOutCMSLC  Res   Vs   QOut
+  if (threeDcuts) qoutVsqoutF.SetPairCut(outPairCut);
+  //analysis->AddResolutionFunction(&qoutVsqoutF);
+
+  AliHBTQSideResolVsQSideFctn qsideVsqsideF;//QSideCMSLC Res   Vs   QSide
+  if (threeDcuts) qsideVsqsideF.SetPairCut(sidePairCut);
+  //analysis->AddResolutionFunction(&qsideVsqsideF);
+
+  
+  AliHBTQLongResolVsQLongFctn qlongVsqlongF;//QLongCMSLC Res   Vs   QLong
+  if (threeDcuts) qlongVsqlongF.SetPairCut(longPairCut);
+  //analysis->AddResolutionFunction(&qlongVsqlongF);
+
+
+  /************************************************************/
+  /************************************************************/
+  /************************************************************/
+  /****   D I S T R B U T I O N S         ***********************/
+  /************************************************************/
+  /************************************************************/
+  /************************************************************/
+
+  /************************************************************/
+  /****   OUT VS KT                 ***********************/
+  /************************************************************/
+  AliHBTQOutDistributionVsKtFctn qoutVsKtDistP;
+  qoutVsKtDistP.Rename("qoutVsKtDistP",qoutVsKtDistP.GetTitle());
+  AliHBTQOutDistributionVsKtFctn qoutVsKtDistT;
+  qoutVsKtDistT.Rename("qoutVsKtDistT",qoutVsKtDistT.GetTitle());
+  if (threeDcuts)  qoutVsKtDistP.SetPairCut(outPairCut);
+  if (threeDcuts)  qoutVsKtDistT.SetPairCut(outPairCut);
+  //analysis->AddParticleFunction(&qoutVsKtDistP);
+  //analysis->AddTrackFunction(&qoutVsKtDistT);
+  /************************************************************/
+  /****   Side VS KT                 ***********************/
+  /************************************************************/
+  AliHBTQSideDistributionVsKtFctn qSideVsKtDistP;
+  qSideVsKtDistP.Rename("qSideVsKtDistP",qSideVsKtDistP.GetTitle());
+  AliHBTQSideDistributionVsKtFctn qSideVsKtDistT;
+  qSideVsKtDistT.Rename("qSideVsKtDistT",qSideVsKtDistT.GetTitle());
+  if (threeDcuts)  qSideVsKtDistP.SetPairCut(sidePairCut);
+  if (threeDcuts)  qSideVsKtDistT.SetPairCut(sidePairCut);
+  //analysis->AddParticleFunction(&qSideVsKtDistP);
+  //analysis->AddTrackFunction(&qSideVsKtDistT);
+  /************************************************************/
+  /****   Long VS KT                 ***********************/
+  /************************************************************/
+  AliHBTQLongDistributionVsKtFctn qLongVsKtDistP;
+  qLongVsKtDistP.Rename("qLongVsKtDistP",qLongVsKtDistP.GetTitle());
+  AliHBTQLongDistributionVsKtFctn qLongVsKtDistT;
+  qLongVsKtDistT.Rename("qLongVsKtDistT",qLongVsKtDistT.GetTitle());
+  if (threeDcuts)  qLongVsKtDistP.SetPairCut(longPairCut);
+  if (threeDcuts)  qLongVsKtDistT.SetPairCut(longPairCut);
+  //analysis->AddParticleFunction(&qLongVsKtDistP);
+  //analysis->AddTrackFunction(&qLongVsKtDistT);
+
+  /************************************************************/
+  /************************************************************/
+  /************************************************************/
+  /****    M O N I T O R  R E S O L U T I O N S   *************/
+  /************************************************************/
+  /************************************************************/
+  /************************************************************/
+  
+  AliHBTMonPxResolutionVsPtFctn PxVsPtResF(200,1.0,0.0,300,0.05,-0.05);
+  //analysis->AddParticleAndTrackMonitorFunction(&PxVsPtResF);
+  
+  AliHBTMonPyResolutionVsPtFctn PyVsPtResF(200,1.0,0.0,300,0.05,-0.05);
+  //analysis->AddParticleAndTrackMonitorFunction(&PyVsPtResF);
+  
+  AliHBTMonPzResolutionVsPtFctn PzVsPtResF(200,1.0,0.0,300,0.05,-0.05);
+  //analysis->AddParticleAndTrackMonitorFunction(&PzVsPtResF);
+
+  AliHBTMonPResolutionVsPtFctn PVsPtResF(200,1.0,0.0,300,0.05,-0.05);
+  //analysis->AddParticleAndTrackMonitorFunction(&PVsPtResF);
+
+  AliHBTMonPtResolutionVsPtFctn PtVsPtResF(400,4.0,0.0,300,0.07,-0.07);
+  analysis->AddParticleAndTrackMonitorFunction(&PtVsPtResF);
+
+  AliHBTMonPhiResolutionVsPtFctn PhiVsPtResF(200,1.0,0.0,300,0.02,-0.02);
+  analysis->AddParticleAndTrackMonitorFunction(&PhiVsPtResF);
+
+  AliHBTMonThetaResolutionVsPtFctn ThetaVsPtResF(200,1.0,0.0,300,0.02,-0.02);
+  analysis->AddParticleAndTrackMonitorFunction(&ThetaVsPtResF);
+  
+  
+  
+  
+  /************************************************************/
+  /****   P R O C E S S                 ***********************/
+  /************************************************************/
+  analysis->SetBufferSize(2);
+  analysis->Process(processopt);
+  
   TFile histoOutput(outfile,"recreate"); 
   analysis->WriteFunctions();
   histoOutput.Close();