Bugs in clusterizing and recon fixed
[u/mrichter/AliRoot.git] / HMPID / AliHMPIDSelector.C
1 #include <TF1.h>
2 #include <TH2F.h>
3 #include <TCanvas.h>  //Terminate()
4 #include <TChain.h>
5 #include <TBenchmark.h>
6 #include <TFile.h>    //docosmic()    
7 #include <fstream>    //caf()      
8 #include <TProof.h>   //caf()
9 #include <AliSelector.h>      //base class
10 #include <AliESD.h>           
11 #include <AliBitPacking.h> //HmpidPayload()
12 #include "AliHMPIDDigit.h" 
13 #include "AliHMPIDCluster.h" 
14 #include "AliHMPIDReconstructor.h" //docosmic()
15
16 class 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 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
47 void 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 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
69 void 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 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
85 Bool_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 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
107 void 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 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
123 void 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
155 void loc()
156 {
157   TChain* pChain =new TChain("esdTree");
158   pChain->Add("AliESDs.root");
159
160   pChain->Process("AliHMPIDSelector.C+");       
161 }
162 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
163 void 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 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
190 Int_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 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
227 void 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 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
238 void HmpidPayload(ifstream *pFile,Int_t iDdl,TObjArray *pDigAll)
239 {
240 // payload is 8 sequences with structure: WC36A8 then WC number of w32
241   
242   TClonesArray *pDig1=(TClonesArray*)pDigAll->At(iDdl/2); //get list of digits for requested chamber
243   
244   UInt_t w32=0;
245   Int_t iDigCnt=pDig1->GetEntriesFast();
246   for(Int_t row=1;row<=8;row++){  
247     pFile->read((char*)&w32,4);  Int_t wc=AliBitPacking::UnpackWord(w32,16,31); Int_t ma=AliBitPacking::UnpackWord(w32, 0,15);    
248     if(ma!=0x36a8) Printf("ERROR ERROR ERROR WRONG Marker=0x%x ERROR ERROR ERROR",ma);    
249     for(Int_t i=1;i<=wc;i++){//words loop
250       pFile->read((char*)&w32,4);
251       if(w32&0x08000000) continue; //it's DILOGIC CW
252       AliHMPIDDigit *pDig=new AliHMPIDDigit;
253       pDig->Raw(iDdl,w32);   
254       new ((*pDig1)[iDigCnt++]) AliHMPIDDigit(*pDig);
255     }//words loop 
256   }//rows loop
257 }//HmpidPayload()
258 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
259 void docosmic(const char* name,Int_t iMaxEvt=9999999)
260 {
261   gBenchmark->Start("Cosmic converter");
262   ifstream in(name);
263
264   
265   Bool_t isPrint=kFALSE;
266   
267   TString rooName=name; rooName.Replace(1+rooName.First('.'),4,"root");
268   Int_t ddl=0;    
269   
270   TFile *pOut=new TFile(rooName,"recreate");
271   TTree *pTr=new TTree("cosmic","some for time being");
272   
273   
274   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); 
275   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);
276
277   Int_t iEvt=0;
278   while(1){      
279     Int_t iSize=DateHeader(&in,      isPrint);  if(iSize==68) continue;  //DATE header 
280     if(in.eof()) break;
281     HmpidHeader (&in,      isPrint);    //HMPID header 
282     HmpidPayload(&in,ddl+1,pDigAll);    //HMPID payload
283     HmpidHeader (&in,      isPrint);    //HMPID header 
284     HmpidPayload(&in,ddl  ,pDigAll);    //HMPID payload
285     
286     AliHMPIDReconstructor::Dig2Clu(pDigAll,pCluAll);
287     pTr->Fill();
288     for(Int_t i=0;i<7;i++){
289       ((TClonesArray*)pDigAll->At(i))->Clear(); 
290       ((TClonesArray*)pCluAll->At(i))->Clear(); 
291     }
292     
293     if(!(iEvt%200)) Printf("Event %i processed",iEvt);
294     iEvt++;
295     if(iEvt==iMaxEvt) break;
296   }//events loop
297   
298   pTr->Write();  pOut->Close();
299   pDigAll->Delete(); pCluAll->Delete();
300   in.close();
301   
302   Printf("Total %i events processed",iEvt);
303   gBenchmark->Show("Cosmic converter");
304 }//docosmic()
305 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
306 void cosmic()
307 {
308   TChain *pCh=new TChain("cosmic");  pCh->Add("test.root");
309
310   
311   TH1F *pDigQ=new TH1F("digQ","Digits QDC",500,0,4100);
312   TH1F *pDigO=new TH1F("digO","Digits # per event",500,0,8000);
313   TH2F *pDigM=new TH2F("digM","Digits map",500,0,131,500,0,127);
314     
315   TH1F *pCluQ   =new TH1F("cluQ","Clusters QDC",500,0,12100);
316   TH1F *pCluQmax=new TH1F("cluQmax","Clusters Max QDC",500,0,12100); pCluQmax->SetLineColor(kRed);
317   TH1F *pCluO=new TH1F("cluO","Clusters # per event",500,0,5000);
318   TH2F *pCluM=new TH2F("cluM","Clusters map",500,0,131,500,0,127);
319   
320   TObjArray *pDigAll=0; pCh->SetBranchAddress("Digs",&pDigAll);
321   TObjArray *pCluAll=0; pCh->SetBranchAddress("Clus",&pCluAll);
322
323     
324   for(Int_t iEvt=0;iEvt<pCh->GetEntries();iEvt++){
325     pCh->GetEntry(iEvt);
326     
327     TClonesArray *pDigCh=(TClonesArray*)pDigAll->At(0);
328     TClonesArray *pCluCh=(TClonesArray*)pCluAll->At(0);
329     pDigO->Fill(pDigCh->GetEntriesFast());
330     pCluO->Fill(pCluCh->GetEntriesFast());
331     
332     for(Int_t iDig=0;iDig<pDigCh->GetEntriesFast();iDig++){//digits loop
333       AliHMPIDDigit *pDig=(AliHMPIDDigit*)pDigCh->UncheckedAt(iDig);
334       pDigQ->Fill(pDig->Q());
335       pDigM->Fill(pDig->LorsX(),pDig->LorsY());
336     }//digits loop
337     Int_t qmax=0;    
338     for(Int_t iClu=0;iClu<pCluCh->GetEntriesFast();iClu++){//clusters loop
339       AliHMPIDCluster *pClu=(AliHMPIDCluster*)pCluCh->UncheckedAt(iClu);
340       pCluQ->Fill(pClu->Q());
341       if(pClu->Q()>qmax) qmax=pClu->Q();
342       pCluM->Fill(pClu->X(),pClu->Y());
343     }//digits loop
344     pCluQmax->Fill(qmax);
345   }//entries loop
346   
347   TCanvas *pC=new TCanvas("comic","cosmic"); pC->Divide(2,3);
348   
349   pC->cd(1); pDigM->Draw(); pC->cd(2); pCluM->Draw();
350   pC->cd(3); gPad->SetLogy(); pDigQ->Draw(); pC->cd(4); gPad->SetLogy(); pCluQ->Draw(); pCluQmax->Draw("same");
351   pC->cd(5); gPad->SetLogy(); pDigO->Draw(); pC->cd(6); gPad->SetLogy(); pCluO->Draw();
352 }//cosmic()
353 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
354