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