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