Small modifications by Alexandru
[u/mrichter/AliRoot.git] / HMPID / AliHMPIDSelector.C
CommitLineData
d3da6dc4 1#include <TF1.h>
2#include <TH2F.h>
3#include <TCanvas.h> //Terminate()
4#include <TChain.h>
5#include <TBenchmark.h>
da08475b 6#include <TFile.h> //docosmic()
d3da6dc4 7#include <fstream> //caf()
8#include <TProof.h> //caf()
9#include <AliSelector.h> //base class
10#include <AliESD.h>
da08475b 11#include <AliBitPacking.h> //HmpidPayload()
12#include "AliHMPIDDigit.h"
13#include "AliHMPIDCluster.h"
14#include "AliHMPIDReconstructor.h" //docosmic()
d3da6dc4 15
16class AliHMPIDSelector : public AliSelector {
17 public :
18 AliHMPIDSelector():AliSelector(),fChain(0),fEsd(0),fCkovP(0),fMipXY(0),fDifXY(0),fSigP(0) {for(Int_t i=0;i<5;i++) fProb[i]=0;}
19 virtual ~AliHMPIDSelector() {delete fEsd;}
20
21
22 virtual Int_t Version () const {return 2;}
23 virtual void Begin (TTree *) {}
24 virtual void SlaveBegin (TTree *tree);
25 virtual void Init (TTree *tree);
26 virtual Bool_t Notify () {return kTRUE;}
27 virtual Bool_t Process (Long64_t entry);
28 virtual void SetOption (const char *option) { fOption = option; }
29 virtual void SetObject (TObject *obj) { fObject = obj; }
30 virtual void SetInputList (TList *input) {fInput = input;}
31 virtual TList *GetOutputList () const { return fOutput; }
32 virtual void SlaveTerminate ();
33 virtual void Terminate ();
34
35 private:
36 TTree *fChain ; //!pointer to the analyzed TTree or TChain
37 AliESD *fEsd ; //!
38
39 TH2F *fCkovP,*fMipXY,*fDifXY,*fSigP; //!
40 TH1F *fProb[5]; //!
41
42 ClassDef(AliHMPIDSelector,0);
43};
44
45
46//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
47void AliHMPIDSelector::SlaveBegin(TTree *tree)
48{
49// The SlaveBegin() function is called after the Begin() function. When running with PROOF SlaveBegin() is called on each slave server.
50// The tree argument is deprecated (on PROOF 0 is passed).
51
52 Init(tree);
53
54 TString option = GetOption();
55
56 // create histograms on each slave server
57 fCkovP = new TH2F("CkovP" , "#theta_{c}, [rad];P, [GeV]", 150, 0, 7 ,100, -3, 1);
58 fSigP = new TH2F("SigP" ,"#sigma_{#theta_c}" , 150, 0, 7 ,100, 0, 1e20);
59 fMipXY = new TH2F("MipXY" ,"mip position" , 260, 0,130 ,252,0,126);
60 fDifXY = new TH2F("DifXY" ,"diff" , 260, -10, 10 ,252,-10,10);
61
62 fProb[0] = new TH1F("PidE" ,"PID: e yellow #mu magenta" ,100,0,1); fProb[0]->SetLineColor(kYellow);
63 fProb[1] = new TH1F("PidMu","pid of #mu" ,100,0,1); fProb[1]->SetLineColor(kMagenta);
64 fProb[2] = new TH1F("PidPi","PID: #pi red K green p blue",100,0,1); fProb[2]->SetLineColor(kRed);
65 fProb[3] = new TH1F("PidK" ,"pid of K" ,100,0,1); fProb[3]->SetLineColor(kGreen);
66 fProb[4] = new TH1F("PidP" ,"pid of p" ,100,0,1); fProb[4]->SetLineColor(kBlue);
67}
68//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
69void AliHMPIDSelector::Init(TTree *pTr)
70{
71 // The Init() function is called when the selector needs to initialize
72 // a new tree or chain. Typically here the branch addresses of the tree
73 // will be set. It is normaly not necessary to make changes to the
74 // generated code, but the routine can be extended by the user if needed.
75 // Init() will be called many times when running with PROOF.
76
77 // Set branch addresses
78 if ( !pTr ) return ;
79 fChain = pTr ;
80 fChain->SetBranchAddress("ESD", &fEsd) ;
81 fChain->SetBranchStatus("*", 0);
82 fChain->SetBranchStatus("fTracks.*", 1);
83}
84//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
85Bool_t AliHMPIDSelector::Process(Long64_t entry)
86{
87
88 fChain->GetTree()->GetEntry(entry);
89
90 for(Int_t iTrk=0;iTrk<fEsd->GetNumberOfTracks();iTrk++){
91 AliESDtrack *pTrk=fEsd->GetTrack(iTrk);
92
93// if(pTrk->GetHMPIDsignal()<0) continue;
94
95 fCkovP->Fill(pTrk->GetP(),pTrk->GetHMPIDsignal()) ;
96 fSigP ->Fill(pTrk->GetP(),TMath::Sqrt(pTrk->GetHMPIDchi2()));
97
98// Float_t xm,ym; Int_t q,np; pTrk->GetHMPIDmip(xm,ym,q,np); fMipXY->Fill(xm,ym); //mip info
99// Float_t xd,yd,th,ph; pTrk->GetHMPIDtrk(xd,yd,th,ph); fDifXY->Fill(xd,yd); //track info
100
101 Double_t pid[5]; pTrk->GetHMPIDpid(pid); for(Int_t i =0;i<5;i++) fProb[i]->Fill(pid[i]);
102 }//tracks loop
103
104 return kTRUE;
105}//Process()
106//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
107void AliHMPIDSelector::SlaveTerminate()
108{
109 // The SlaveTerminate() function is called after all entries or objects
110 // have been processed. When running with PROOF SlaveTerminate() is called
111 // on each slave server.
112
113 // Add the histograms to the output on each slave server
114
115 fOutput->Add(fCkovP);
116 fOutput->Add(fSigP);
117 fOutput->Add(fMipXY);
118 fOutput->Add(fDifXY);
119
120 for(Int_t i=0;i<5;i++) fOutput->Add(fProb[i]);
121}//SlaveTerminate()
122//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
123void AliHMPIDSelector::Terminate()
124{
125 // The Terminate() function is the last function to be called during
126 // a query. It always runs on the client, it can be used to present
127 // the results graphically or save the results to file.
128
129 fCkovP = dynamic_cast<TH2F*>(fOutput->FindObject("CkovP")) ;
130 fSigP = dynamic_cast<TH2F*>(fOutput->FindObject("SigP")) ;
131 fMipXY = dynamic_cast<TH2F*>(fOutput->FindObject("MipXY")) ;
132 fDifXY = dynamic_cast<TH2F*>(fOutput->FindObject("DifXY")) ;
133
134 fProb[0] = dynamic_cast<TH1F*>(fOutput->FindObject("PidE")) ;
135 fProb[1] = dynamic_cast<TH1F*>(fOutput->FindObject("PidMu")) ;
136 fProb[2] = dynamic_cast<TH1F*>(fOutput->FindObject("PidPi")) ;
137 fProb[3] = dynamic_cast<TH1F*>(fOutput->FindObject("PidK")) ;
138 fProb[4] = dynamic_cast<TH1F*>(fOutput->FindObject("PidP")) ;
139
140 Float_t n=1.292; //mean freon ref idx
141 TF1 *pPi=new TF1("RiPiTheo","acos(sqrt(x*x+[0]*[0])/(x*[1]))",1.2,7); pPi->SetLineWidth(1); pPi->SetParameter(1,n);
142 AliPID ppp; pPi->SetLineColor(kRed); pPi->SetParameter(0,AliPID::ParticleMass(AliPID::kPion)); //mass
143 TF1 *pK=(TF1*)pPi->Clone(); pK ->SetLineColor(kGreen); pK ->SetParameter(0,AliPID::ParticleMass(AliPID::kKaon));
144 TF1 *pP=(TF1*)pPi->Clone(); pP ->SetLineColor(kBlue); pP ->SetParameter(0,AliPID::ParticleMass(AliPID::kProton));
145
146 TCanvas *pC=new TCanvas("c1","ESD QA");pC->SetFillColor(10); pC->SetHighLightColor(10); pC->Divide(3,2);
147 pC->cd(1); fCkovP->Draw(); pPi->Draw("same"); pK->Draw("same"); pP->Draw("same"); pC->cd(2); fMipXY->Draw(); pC->cd(3); fProb[0]->Draw(); fProb[1]->Draw("same");
148 pC->cd(4); fSigP ->Draw(); pC->cd(5); fDifXY->Draw(); pC->cd(6); fProb[2]->Draw(); fProb[3]->Draw("same"); fProb[4]->Draw("same");
149
150}
151//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
152
153
154
155void loc()
156{
157 TChain* pChain =new TChain("esdTree");
158 pChain->Add("AliESDs.root");
159
160 pChain->Process("AliHMPIDSelector.C+");
161}
162//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
163void caf()
164{
165 gBenchmark->Start("PRooF exec");
166 TChain* pChain =new TChain("esdTree");
167
168 ifstream list; list.open("list.txt");
169
170 TString file;
171 while(list.good()) {
172 list>>file;
173 if (!file.Contains("root")) continue; //it's wrong file name
174 pChain->Add(file.Data());
175 }
176 list.close();
177
178 pChain->GetListOfFiles()->Print();
179
180 TVirtualProof *pProof=TProof::Open("kir@lxb6046.cern.ch");
181 pProof->UploadPackage("ESD.par");
182 pProof->EnablePackage("ESD");
183
184 pChain->SetProof(pProof);
185 pChain->Process("AliHMPIDSelector.C+");
186
187 gBenchmark->Show("PRooF exec");
188}
189//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
da08475b 190Int_t DateHeader(ifstream *pFile,Bool_t isPrint=0)
191{
192 Int_t iSize=-1;
193 pFile->read((char*)&iSize,4);
194 if(!isPrint)
195 pFile->seekg(16*4,ios::cur);
196 else{
197 Int_t w32=-1;
198 Printf("");
199 Printf("Event size %i bytes",iSize); //1 common DATE header 17 words
200 pFile->read((char*)&w32,4); Printf("Event magic 0x%x" ,w32); //2
201 pFile->read((char*)&w32,4); Printf("Event head size %i bytes",w32); //3
202 pFile->read((char*)&w32,4); Printf("Event version 0x%x" ,w32); //4
203 pFile->read((char*)&w32,4); Printf("Event type %i (%s)" ,w32,(w32==7)? "physics":"SOR"); //5
204 pFile->read((char*)&w32,4); Printf("Run number %i" ,w32); //6
205
206 pFile->read((char*)&w32,4); Printf("Event ID 1 %i" ,w32); //7
207 pFile->read((char*)&w32,4); Printf("Event ID 2 %i" ,w32); //8
208
209 pFile->read((char*)&w32,4); Printf("Trigger pattern 1 %i" ,w32); //9
210 pFile->read((char*)&w32,4); Printf("Trigger pattern 2 %i" ,w32); //10
211
212 pFile->read((char*)&w32,4); Printf("Detector pattern %i" ,w32); //11
213
214 pFile->read((char*)&w32,4); Printf("Type attribute 1 %i" ,w32); //12
215 pFile->read((char*)&w32,4); Printf("Type attribute 2 %i" ,w32); //13
216 pFile->read((char*)&w32,4); Printf("Type attribute 3 %i" ,w32); //14
217
218 pFile->read((char*)&w32,4); Printf("LDC ID %i" ,w32); //15
219 pFile->read((char*)&w32,4); Printf("GDC ID %i" ,w32); //16
220 pFile->read((char*)&w32,4); TDatime time(w32); time.Print(); //17
221
222 Printf("");
223 }
224 return iSize;
225}
226//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
227void HmpidHeader(ifstream *pFile,Bool_t isPrint=kFALSE)
228{
229// Prints hmpid trailer and returns number of words for this trailer
230 if(!isPrint) {pFile->seekg(15*4,ios::cur);return;}
231 Int_t w32=-1;
232 Printf("\nHMPID Header:");//private HEADER is 15 words
233 for(Int_t i=1;i<=11;i++) { pFile->read((char*)&w32,4); Printf("Word #%2i=%12i meaningless",i,w32);}
234 pFile->read((char*)&w32,4); Printf("Word #12=%12i event counter",w32);
235 for(Int_t i=13;i<=15;i++){ pFile->read((char*)&w32,4); Printf("Word #%2i=%12i empty",i,w32);}
236}
237//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
238void HmpidPayload(ifstream *pFile,Int_t iDdl,TClonesArray *pDigLst)
239{
240// payload is 8 sequences with structure: WC36A8 then WC number of w32
241 UInt_t w32=0;
242 Int_t iDigCnt=pDigLst->GetEntriesFast();
243 for(Int_t row=1;row<=8;row++){
244 pFile->read((char*)&w32,4); Int_t wc=AliBitPacking::UnpackWord(w32,16,31); Int_t ma=AliBitPacking::UnpackWord(w32, 0,15);
245 if(ma!=0x36a8) Printf("ERROR ERROR ERROR WRONG Marker=0x%x ERROR ERROR ERROR",ma);
246 for(Int_t i=1;i<=wc;i++){//words loop
247 pFile->read((char*)&w32,4);
248 if(w32&0x08000000) continue; //it's DILOGIC CW
249 AliHMPIDDigit *pDig=new AliHMPIDDigit;
250 pDig->Raw(iDdl,w32);
251 new ((*pDigLst)[iDigCnt++]) AliHMPIDDigit(*pDig);
252 }//words loop
253 }//rows loop
254}//HmpidPayload()
255//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
256void docosmic(const char* name)
257{
258 gBenchmark->Start("Cosmic converter");
259 ifstream in(name);
260
261
262 Bool_t isPrint=kFALSE;
263
264 TString rooName=name; rooName.Replace(1+rooName.First('.'),4,"root");
265 Int_t ddl=0;
266
267 TFile *pOut=new TFile(rooName,"recreate");
268 TTree *pTr=new TTree("cosmic","some for time being");
269
270 TClonesArray *pDigLst=new TClonesArray("AliHMPIDDigit"); pTr->Branch("Digs",&pDigLst);
271 TClonesArray *pCluLst=new TClonesArray("AliHMPIDCluster"); pTr->Branch("Clus",&pCluLst);
272
273 while(1){
274 Int_t iSize=DateHeader(&in, isPrint); if(iSize==68) continue; //DATE header
275 if(in.eof()) break;
276 HmpidHeader (&in, isPrint); //HMPID header
277 HmpidPayload(&in,ddl+1,pDigLst); //HMPID payload
278 HmpidHeader (&in, isPrint); //HMPID header
279 HmpidPayload(&in,ddl ,pDigLst); //HMPID payload
280
281 AliHMPIDReconstructor::Dig2Clu(pDigLst,pCluLst);
282 pTr->Fill();
283 pDigLst->Clear(); pCluLst->Clear();
284 }
285
286 pTr->Write();
287 pOut->Close();
288 delete pDigLst;
289 in.close();
290 gBenchmark->Show("Cosmic converter");
291}//docosmic()
292//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
293void cosmic()
294{
295 TChain *pCh=new TChain("cosmic");
296 pCh->Add("test.root");
297
298
299 TH1F *pDigQ=new TH1F("digQ","Digits QDC",500,0,4100);
300 TH1F *pDigO=new TH1F("digO","Digits # per event",500,0,8000);
301 TH2F *pDigM=new TH2F("digM","Digits map",500,0,131,500,0,127);
302
303 TH1F *pCluQ=new TH1F("cluQ","Clusters QDC",500,0,12100);
304 TH1F *pCluO=new TH1F("cluO","Clusters # per event",500,0,5000);
305 TH2F *pCluM=new TH2F("cluM","Clusters map",500,0,131,500,0,127);
306
307 TClonesArray *pDigLst=new TClonesArray("AliHMPIDDigit"); pCh->SetBranchAddress("Digs",&pDigLst);
308 TClonesArray *pCluLst=new TClonesArray("AliHMPIDCluster"); pCh->SetBranchAddress("Clus",&pCluLst);
309
310 for(Int_t iEvt=0;iEvt<pCh->GetEntries();iEvt++){
311 pCh->GetEntry(iEvt);
312
313 pDigO->Fill(pDigLst->GetEntriesFast());
314 pCluO->Fill(pCluLst->GetEntriesFast());
315
316 for(Int_t iDig=0;iDig<pDigLst->GetEntriesFast();iDig++){//digits loop
317 AliHMPIDDigit *pDig=(AliHMPIDDigit*)pDigLst->UncheckedAt(iDig);
318 pDigQ->Fill(pDig->Q());
319 pDigM->Fill(pDig->LorsX(),pDig->LorsY());
320 }//digits loop
321
322 for(Int_t iClu=0;iClu<pCluLst->GetEntriesFast();iClu++){//clusters loop
323 AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluLst->UncheckedAt(iClu);
324 pCluQ->Fill(pClu->Q());
325 pCluM->Fill(pClu->X(),pClu->Y());
326 }//digits loop
327 }//entries loop
328
329 TCanvas *pC=new TCanvas("comic","cosmic"); pC->Divide(2,3);
330
331 pC->cd(1); pDigM->Draw(); pC->cd(2); pCluM->Draw();
332 pC->cd(3); pDigQ->Draw(); pC->cd(4); pCluQ->Draw();
333 pC->cd(5); pDigO->Draw(); pC->cd(6); pCluO->Draw();
334}//cosmic()
335//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
d3da6dc4 336