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