6e3ae1887f856fca3f5023c0099be6ebbfad6f4e
[u/mrichter/AliRoot.git] / HBTAN / hbtanalysis.C
1 void hbtanalysis(Option_t* datatype, Option_t* processopt="TracksAndParticles", 
2                 Int_t first = -1,Int_t last = -1, 
3                 char *outfile = "hbtanalysis.root")
4  {
5 //HBT Anlysis Macro
6 //Anlyzes TPC recontructed tracks and simulated particles that corresponds to them
7
8 //datatype defines type of data to be read
9 //  Kine  - analyzes Kine Tree: simulated particles
10 //  TPC   - analyzes TPC   tracking + particles corresponding to these tracks
11 //  ITSv1 - analyzes ITSv1 ----------------------//--------------------------
12 //  ITSv2 - analyzes ITSv2 ----------------------//--------------------------
13
14 //processopt defines option passed to AliHBTAnlysis::Process method
15 // default: TracksAndParticles - process both recontructed tracks and sim. particles corresponding to them 
16 //          Tracks - process only recontructed tracks
17 //          Particles - process only simulated particles
18
19 //Reads data from diroctories from first to last(including)
20 // For examples if first=3 and last=5 it reads from
21 //  ./3/
22 //  ./4/
23 //  ./5/
24 //if first or last is negative (or both), it reads from current directory
25 //
26 //these names I use when analysis is done directly from CASTOR files via RFIO
27 //  const char* basedir="rfio:/castor/cern.ch/user/s/skowron/";
28 //  const char* serie="standard";
29 //  const char* field = "0.4";
30   const char* basedir=".";
31   const char* serie="";
32   const char* field = "";
33  
34   // Dynamically link some shared libs                    
35   gROOT->LoadMacro("loadlibs.C");
36   loadlibs();
37   
38   gSystem->Load("$(ALICE_ROOT)/lib/tgt_$(ALICE_TARGET)/libHBTAN");
39   
40
41   /***********************************************************/
42   //Create analysis and reader
43   AliHBTAnalysis * analysis = new AliHBTAnalysis();
44   
45   AliHBTReader* reader;
46   Int_t kine = strcmp(datatype,"Kine");
47   Int_t TPC = strcmp(datatype,"TPC");
48   Int_t ITSv1 = strcmp(datatype,"ITSv1");
49   Int_t ITSv2 = strcmp(datatype,"ITSv2");
50   Int_t intern = strcmp(datatype,"Intern");
51   
52   if(!kine)
53    {
54     reader = new AliHBTReaderKineTree();
55     processopt="Particles"; //this reader by definition reads only simulated particles
56    }
57   else if(!TPC)
58    {
59     reader = new AliHBTReaderTPC();
60    }
61   else if(!ITSv1)
62    {
63     reader = new AliHBTReaderITSv1();
64    }
65   else if(!ITSv2)
66    {
67     reader = new AliHBTReaderITSv2();
68    }
69   else if(!intern)
70    {
71     reader = new AliHBTReaderInternal("Kine.sum.root");
72    }
73   else
74    {
75     cerr<<"Option "<<datatype<<"  not recognized. Exiting"<<endl;
76     return;
77    }
78
79   TObjArray* dirs=0;
80   if ( (first >= 0) && (last>=0) && ( (last-first)>=0 ) )
81    {//read from many dirs dirs
82      char buff[50];
83      dirs = new TObjArray(last-first+1);
84      for (Int_t i = first; i<=last; i++)
85       {
86         sprintf(buff,"%s/%s/%s/%d",basedir,field,serie,i);
87         TObjString *odir= new TObjString(buff);
88         dirs->Add(odir);
89       }
90     }
91    reader->SetDirs(dirs);
92    
93   /***********************************************************/    
94   
95   //we want have only low pt pi+ so set a cut to reader
96   AliHBTParticleCut* readerpartcut= new AliHBTParticleCut();
97   readerpartcut->SetPtRange(0.0,1.5);
98   readerpartcut->SetPID(kPiPlus);
99   reader->AddParticleCut(readerpartcut);//read this particle type with this cut
100   
101   analysis->SetReader(reader);
102   /************************************************************/
103   
104   AliHBTPairCut *paircut = new AliHBTPairCut();
105   paircut->SetQInvRange(0.0,0.20);  
106   analysis->SetGlobalPairCut(paircut);
107
108   AliHBTQInvCorrelFctn * qinvcfT= new AliHBTQInvCorrelFctn();
109   AliHBTQInvCorrelFctn * qinvcfP= new AliHBTQInvCorrelFctn();
110   
111   analysis->AddTrackFunction(qinvcfT);
112   analysis->AddParticleFunction(qinvcfP);
113   
114   
115   qinvcfP->Rename("qinvcfP","Particle (simulated) Qinv CF \\pi^{+} \\pi^{+}");
116   qinvcfT->Rename("qinvcfT","Track (recontructed) Qinv CF \\pi^{+} \\pi^{+}");
117   
118   AliHBTQOutCMSLCCorrelFctn* qoutP = new AliHBTQOutCMSLCCorrelFctn();
119   qoutP->Rename("qoutP","Particle (simulated) Q_{out} CF \\pi^{+} \\pi^{+}");
120   AliHBTQOutCMSLCCorrelFctn* qoutT = new AliHBTQOutCMSLCCorrelFctn(); 
121   qoutT->Rename("qoutT","Track (recontructed) Q_{out} CF \\pi^{+} \\pi^{+}");
122
123   AliHBTPairCut *outPairCut = new AliHBTPairCut();
124   outPairCut->SetQOutCMSLRange(0.0,0.15);
125   outPairCut->SetQSideCMSLRange(0.0,0.02);
126   outPairCut->SetQLongCMSLRange(0.0,0.02);
127   qoutP->SetPairCut(outPairCut);
128   qoutT->SetPairCut(outPairCut);
129   
130   AliHBTQSideCMSLCCorrelFctn* qsideP = new AliHBTQSideCMSLCCorrelFctn(); 
131   qsideP->Rename("qsideP","Particle (simulated) Q_{side} CF \\pi^{+} \\pi^{+}");
132   AliHBTQSideCMSLCCorrelFctn* qsideT = new AliHBTQSideCMSLCCorrelFctn(); 
133   qsideT->Rename("qsideT","Track (recontructed) Q_{side} CF \\pi^{+} \\pi^{+}");
134
135   AliHBTPairCut *sidePairCut = new AliHBTPairCut();
136   sidePairCut->SetQOutCMSLRange(0.0,0.02);
137   sidePairCut->SetQSideCMSLRange(0.0,0.15);
138   sidePairCut->SetQLongCMSLRange(0.0,0.02);
139   qsideP->SetPairCut(sidePairCut);
140   qsideT->SetPairCut(sidePairCut);
141
142     
143   AliHBTQLongCMSLCCorrelFctn* qlongP = new AliHBTQLongCMSLCCorrelFctn(); 
144   qlongP->Rename("qlongP","Particle (simulated) Q_{long} CF \\pi^{+} \\pi^{+}");
145   AliHBTQLongCMSLCCorrelFctn* qlongT = new AliHBTQLongCMSLCCorrelFctn(); 
146   qlongT->Rename("qlongT","Track (recontructed) Q_{long} CF \\pi^{+} \\pi^{+}");
147
148   AliHBTPairCut *longPairCut = new AliHBTPairCut();
149   longPairCut->SetQOutCMSLRange(0.0,0.02);
150   longPairCut->SetQSideCMSLRange(0.0,0.02);
151   longPairCut->SetQLongCMSLRange(0.0,0.15);
152   qlongP->SetPairCut(longPairCut);
153   qlongT->SetPairCut(longPairCut);
154
155   analysis->AddTrackFunction(qoutP);
156   analysis->AddParticleFunction(qoutT);
157   
158   analysis->AddTrackFunction(qsideP);
159   analysis->AddParticleFunction(qsideT);
160   
161   analysis->AddTrackFunction(qlongT);
162   analysis->AddParticleFunction(qlongT);
163   
164
165   analysis->Process(processopt);
166   
167   TFile histoOutput(outfile,"recreate"); 
168   analysis->WriteFunctions();
169   histoOutput.Close();
170   
171   delete qinvcfP;
172   delete qinvcfT;
173   delete paircut;
174   delete readerpartcut;
175   if (dirs) 
176    {
177     dirs->SetOwner();
178     delete dirs;
179    }
180   delete reader;
181   delete analysis;
182   
183  }
184