PWGCF/FEMTOSCOPY/macros/Local/AliFemtoSplittingMergingQA.C

   1 // AliFemtoSplittingMergingQA.C - macro to create the splitting/merging
   2 // test with the pion correlation function.
   3 // As a default the anti-splitting and anti-merging cuts are open
   4 // and the two correlation functions:
   5 // AliFemtoShareQualityCorrFctn and AliFemtoTPCInnerCorrFctn
   6 // can be used to study the splitting (former) and merging (latter) effect
   7 // If ones needs to produce a "clean" sample with both effects removed,
   8 // one needs to change the cut values to the "reasonable" ones, or perform
   9 // the full systematic analysis with the above-mentioned functions
  10
  11 // Author: Adam Kisiel. Adam.Kisiel@cern.ch
  12
  13 // parameters:
  14 // listFileName - a text file containing a list of ESDs (with full paths)
  15 // to analyze
  16 void AliFemtoSplittingMergingQA(const char *listFileName)
  17 {
  18   double PionMass = 0.13956995;
  19   double chargePi = 1.0;
  20
  21   // Load the neccessary libraries
  22   gSystem->Load("libTree");
  23   gSystem->Load("libGeom");
  24   gSystem->Load("libEG");
  25   gSystem->Load("libCint");
  26   gSystem->Load("libESD");
  27   gSystem->Load("libPWG2femtoscopy");
  28   gSystem->Load("libPWG2femtoscopyUser");
  29
  30   // Set-up the reader for ALICE ESD
  31   AliFemtoEventReaderESD* Reader=new AliFemtoEventReaderESD();
  32   // Read only constrained momenta - primordial particles
  33   Reader->SetConstrained(true);
  34   Reader->SetReadTPCInner(true);
  35   // Read the file list from the filename supplied by the user
  36   Reader->SetInputFile(listFileName);
  37
  38   // Setup the manager
  39   AliFemtoManager* Manager=new AliFemtoManager();
  40   // Point to the data source - the reader
  41   Manager->SetEventReader(Reader);
  42
  43   // Setup the analysis
  44   AliFemtoSimpleAnalysis* an =new AliFemtoSimpleAnalysis();
  45   // Number of events to construct the background
  46   an->SetNumEventsToMix(3);
  47
  48   // The event selector
  49   AliFemtoBasicEventCut* mec = new AliFemtoBasicEventCut();
  50   // Accept events with the given multiplicity
  51   mec->SetEventMult(0,100000);
  52   // and z-vertex distance to the center of the TPC
  53   mec->SetVertZPos(-1000,1000);
  54
  55   // The track selector
  56   AliFemtoESDTrackCut* dtc = new AliFemtoESDTrackCut();
  57   // We want positive pions
  58   dtc->SetPidProbPion(0.2,1.001);
  59   dtc->SetPidProbMuon(0.0,0.8);
  60   dtc->SetPidProbKaon(0.0,0.1);
  61   dtc->SetPidProbProton(0.0,0.1);
  62   dtc->SetMostProbablePion();
  63   dtc->SetCharge(1.0);
  64   // so we set the correct mass
  65   dtc->SetMass(PionMass);
  66   // we select low pt
  67   dtc->SetPt(0.1,0.7);
  68   dtc->SetStatus(AliESDtrack::kTPCrefit|AliESDtrack::kITSrefit);
  69   dtc->SetminTPCncls(95);
  70   dtc->SetRemoveKinks(kTRUE);
  71   dtc->SetLabel(kFALSE);
  72   dtc->SetMaxITSChiNdof(3.0);
  73   dtc->SetMaxTPCChiNdof(2.0);
  74   dtc->SetMaxSigmaToVertex(3.0);
  75
  76   AliFemtoCutMonitorParticleYPt *cutPass = new AliFemtoCutMonitorParticleYPt("cutPass", 0.13957);
  77   AliFemtoCutMonitorParticleYPt *cutFail = new AliFemtoCutMonitorParticleYPt("cutFail", 0.13957);
  78   dtc->AddCutMonitor(cutPass, cutFail);
  79
  80   // Pair selector
  81   AliFemtoShareQualityTPCEntranceSepPairCut *sqpc = new AliFemtoShareQualityTPCEntranceSepPairCut();
  82   // remove split track pairs and pairs that share hits
  83
  84   // Set maximim allowed "quality" for the pair
  85   //  1.0 - accept all pairs
  86   // -0.5 - reject all pairs
  87   // a reasonable value should lie between 0.0 and 0.5
  88   sqpc->SetShareQualityMax(1.0);
  89
  90   // Set maximum allowed shared hits fraction per pair
  91   //  1.0 - accept all pairs
  92   //  0.0 - reject all pairs
  93   // a reasonable value is small but nno-zero (0.05)
  94   sqpc->SetShareFractionMax(1.0);
  95
  96   // Set minimum allowed separation between nominal TPC entrance points
  97   // of the two tracks in the pair
  98   // 0.0 - accept all pairs
  99   // a reasonable value is 3.0 [cm]
 100   sqpc->SetTPCEntranceSepMinimum(0.0);
 101   sqpc->SetRemoveSameLabel(kFALSE);
 102
 103   // Add the cuts to the analysis
 104   an->SetEventCut(mec);
 105   an->SetFirstParticleCut(dtc);
 106   an->SetSecondParticleCut(dtc);
 107   an->SetPairCut(sqpc);
 108
 109   // Setup correlation functions
 110   // A simple qinv correlation function
 111   AliFemtoQinvCorrFctn *cqinv= new AliFemtoQinvCorrFctn("qinvcf",75,0.0,0.4);
 112
 113   // A correlation function to monitor the splitting and cluster sharing in TPC
 114   AliFemtoShareQualityCorrFctn *csqqinv= new AliFemtoShareQualityCorrFctn("sqqinvcf",75,0.0,0.4);
 115
 116   // A correlation function to monitor the distance at the entrance to the TPC
 117   AliFemtoTPCInnerCorrFctn *tpcin = new AliFemtoTPCInnerCorrFctn("tpcin",80, 0.0, 0.4);
 118
 119   // add the correlation functions to the analysis
 120   an->AddCorrFctn(cqinv);
 121   an->AddCorrFctn(csqqinv);
 122   an->AddCorrFctn(tpcin);
 123
 124   // Add the analysis to the manager
 125   Manager->AddAnalysis(an);
 126
 127   // Run the event loop
 128   long nE= 100000;
 129   if (Manager->Init())
 130     cout<<" Problem"<<endl;
 131
 132   int Status=0;
 133   long int nEP=0;
 134   while ((!Status)&&(nEP<nE))
 135     {
 136       nEP++;
 137       cout<<" next event "<<nEP<<endl;
 138       Status=Manager->ProcessEvent();
 139     }
 140
 141   // Save the results
 142   char ofname[200];
 143   sprintf(ofname, "QinvCF.In.root");
 144
 145   TFile f1 (ofname,"RECREATE","Data");
 146   cqinv->Numerator()->Write();
 147   cqinv->Denominator()->Write();
 148   csqqinv->WriteHistos();
 149   tpcin->WriteHistos();
 150   cutPass->Write();
 151   cutFail->Write();
 152
 153   // Save the cut settings in the output file
 154   f1.mkdir("Settings");
 155   f1.cd("Settings");
 156
 157   TList *tListSettings = an->ListSettings();
 158   tListSettings->Write();
 159
 160   f1.Close();
 161
 162   TFile f2("Listout.root","RECREATE");
 163   TList *tOutList = an->GetOutputList();
 164   tOutList->Write();
 165 }