HMPID/Hmenu.C

   1 AliRun *a;     AliRunLoader *al;   TGeoManager *g; //globals for easy manual manipulations
   2 AliHMPID   *r; AliLoader    *rl; AliHMPIDParam *rp;
   3 Bool_t isGeomType=kFALSE;
   4
   5 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   6 void GetParam()
   7 {
   8   isGeomType=!isGeomType;
   9   if(g) delete g;  if(rp) delete rp; //delete current TGeoManager and AliHMPIDParam
  10   if(isGeomType) g=TGeoManager::Import("geometry.root");
  11   else           g=TGeoManager::Import("misaligned_geometry.root");
  12   rp=AliHMPIDParam::Instance();
  13 }
  14 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  15 void Hmenu()
  16 {
  17   TString status="Status: ";
  18   if(gSystem->IsFileInIncludePath("galice.root")){
  19     status+="galice.root found";
  20     al=AliRunLoader::Open();                                                //try to open galice.root from current dir
  21     if(gAlice) delete gAlice;                                               //in case we execute this in aliroot delete default AliRun object
  22     al->LoadgAlice(); a=al->GetAliRun();                                    //take new AliRun object from galice.root
  23     rl=al->GetDetectorLoader("HMPID");  r=(AliHMPID*)a->GetDetector("HMPID");  //get HMPID object from galice.root
  24
  25     status+=Form(" with %i event(s)",al->GetNumberOfEvents()); status+=(r)? " with HMPID": " without HMPID";
  26   }else
  27     status+="No galice.root";
  28
  29   GetParam();
  30
  31   TControlBar *pMenu = new TControlBar("horizontal",status.Data(),0,0);
  32     pMenu->AddButton("                     ","","");
  33     pMenu->AddButton(" General  "           ,"General()"  ,"general items which do not depend on any files");
  34     pMenu->AddButton("                     ",""           ,"");
  35     pMenu->AddButton(" Sim data "           ,"SimData()"  ,"items which expect to have simulated files"    );
  36     pMenu->AddButton("                     ",""           ,"");
  37     pMenu->AddButton(" Raw data "           ,"RawData()"  ,"items which expect to have raw files"          );
  38     pMenu->AddButton("                     ","print()"    ,"");
  39     pMenu->AddButton("Test"                 ,"Test()"     ,"all test utilities");
  40     pMenu->AddButton("                     ","GetParam()" ,"");
  41     pMenu->AddButton("Quit"                 ,".q"         ,"close session"                                 );
  42   pMenu->Show();
  43 }
  44 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  45 void General()
  46 {
  47   TControlBar *pMenu = new TControlBar("vertical","Sim data",60,50);
  48     pMenu->AddButton("Debug ON"           ,"don();"                    ,"Switch debug on-off"                        );
  49     pMenu->AddButton("Debug OFF"          ,"doff();"                   ,"Switch debug on-off"                        );
  50     pMenu->AddButton("Geo GUI"            ,"geo();"                    ,"Shows geometry"                             );
  51     pMenu->AddButton("Browser"            ,"new TBrowser;"             ,"Start ROOT TBrowser"                        );
  52   pMenu->Show();
  53 }//menu()
  54 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  55 void SimData()
  56 {
  57   TControlBar *pMenu = new TControlBar("vertical","Sim",190,50);
  58     pMenu->AddButton("Display ALL chambers"            ,"ed();"     , "Display Fast");
  59     pMenu->AddButton("HITS QA"                         ,"hqa()"     ,"QA plots for hits: hqa()");
  60
  61     pMenu->AddButton("Print hits"                      ,"hp();"      ,"To print hits:     hp(evt)");
  62     pMenu->AddButton("Print sdigits"                   ,"sp();"      ,"To print sdigits:  sp(evt)");
  63     pMenu->AddButton("Print digits"                    ,"dp();"      ,"To print digits:   dp(evt)");
  64     pMenu->AddButton("Print clusters"                  ,"cp();"      ,"To print clusters: cp(evt)");
  65
  66   pMenu->Show();
  67 }
  68 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  69 void RawData()
  70 {
  71   TControlBar *pMenu = new TControlBar("vertical","Raw",350,50);
  72     pMenu->AddButton("ESD print"                       ,"ep();"                  ,"To print ESD info: ep()"         );
  73     pMenu->AddButton("ESD QA"                          ,"eq();"                  ,"To draw ESD hists: eq()"         );
  74     pMenu->AddButton("Clusters print"                  ,"cp();"                  ,"To print clusters: cp()"         );
  75     pMenu->AddButton("Clusters QA"                     ,"cq();"                  ,"To draw clusters hists: cq()"    );
  76     pMenu->AddButton("Print Matrix"                    ,"mp();"                  ,"To print prob matrix: mp()"      );
  77     pMenu->AddButton("Print occupancy"                 ,"r->OccupancyPrint(-1);" ,"To print occupancy"              );
  78     pMenu->AddButton("Print event summary  "           ,"r->SummaryOfEvent();"   ,"To print a summary of the event" );
  79   pMenu->Show();
  80 }//RawData
  81 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  82 void Test()
  83 {
  84   TControlBar *pMenu = new TControlBar("vertical","Test",400,50);
  85     pMenu->AddButton("Hits->Digits"       ,"thd();"                    ,"test hits->sdigits->digits"                 );
  86     pMenu->AddButton("Segmentation"       ,"ts()"                      ,"test segmentation methods"                  );
  87     pMenu->AddButton("Test response"      ,"AliHMPIDParam::TestResp();","Test AliHMPIDParam response methods"         );
  88     pMenu->AddButton("Print map"          ,"PrintMap();"               ,"Test AliHMPIDParam transformation methods"   );
  89     pMenu->AddButton("Test Recon"         ,"rec();"                    ,"Test AliHMPIDRecon"                          );
  90   pMenu->Show();
  91 }//Test()
  92 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  93
  94
  95 void doff(){  Printf("DebugOFF");  AliLog::SetGlobalDebugLevel(0);}
  96 void don() {  Printf("DebugON");   AliLog::SetGlobalDebugLevel(AliLog::kDebug);}
  97
  98 void geo (                       ) {gGeoManager->GetTopVolume()->Draw("ogl");}
  99
 100 void du  (                       ) {r->Dump         (   );}                //utility display
 101
 102 void hp  (Int_t evt=0            ) {r->HitPrint  (evt);}   //print hits for requested event
 103 void hq  (                       ) {r->HitQA     (   );}   //hits QA plots for all events
 104 void sp  (Int_t evt=0            ) {r->SdiPrint  (evt);}   //print sdigits for requested event
 105 void sq  (Int_t evt=0            ) {r->SdiPrint  (evt);}   //print sdigits for requested event
 106 void dp  (Int_t evt=0            ) {r->DigPrint  (evt);}   //print digits for requested event
 107 void dq  (                       ) {AliHMPIDReconstructor::DigQA     (al );}   //digits QA plots for all events
 108 void cp  (Int_t evt=0            ) {r->CluPrint  (evt);                   }   //print clusters for requested event
 109 void cq  (                       ) {AliHMPIDReconstructor::CluQA     (al );}   //clusters QA plots for all events
 110
 111 void ep  (                       ) {AliHMPIDTracker::EsdQA(1);            }
 112 void eq  (                       ) {AliHMPIDTracker::EsdQA();             }
 113 void mp  (Double_t probCut=0.7   ) {AliHMPIDTracker::MatrixPrint(probCut);}
 114
 115
 116 void t   (Int_t evt=0          )   {AliHMPIDParam::Stack(evt);}
 117 void tid (Int_t tid,Int_t evt=0)   {AliHMPIDParam::Stack(evt,tid);}
 118 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 119 void tst()
 120 {
 121 //create all lists
 122   TClonesArray hit("AliHMPIDHit"),sdi("AliHMPIDDigit");
 123   TObjArray    dig,clu,cog;      for(Int_t i=0;i<7;i++) {dig.AddAt(new TClonesArray("AliHMPIDDigit"),i);
 124                                                          clu.AddAt(new TClonesArray("AliHMPIDCluster"),i);
 125                                                          cog.AddAt(new TClonesArray("AliHMPIDCluster"),i);}
 126 //simulate track
 127   TLorentzVector mom; mom.SetPtEtaPhiM(3,0,30*TMath::DegToRad(),0.938);
 128   TLorentzVector vtx; vtx.SetXYZT(0,0,0,0);
 129
 130   AliESD *pESD=new AliESD;  pESD->SetMagneticField(0.2);
 131   pESD->AddTrack(new AliESDtrack(new TParticle(kProton,0,-1,-1,-1,-1,mom,vtx)));
 132   pESD->Print();
 133
 134   AliTracker::SetFieldMap(new AliMagF,1);
 135   AliHMPIDTracker *pTracker=new AliHMPIDTracker;
 136
 137
 138   return;
 139   Double_t th=8*TMath::DegToRad();                              //gRandom->Rndm()*TMath::PiOver4();
 140   Double_t ph=gRandom->Rndm()*TMath::TwoPi();
 141   Double_t radx=gRandom->Rndm()*AliHMPIDDigit::SizeAllX();
 142   Double_t rady=gRandom->Rndm()*AliHMPIDDigit::SizeAllY();
 143
 144   Int_t iHitCnt=0;
 145
 146
 147
 148
 149
 150
 151              AliHMPIDv1::Hit2Sdi(&hitLst,&sdiLst);
 152       AliHMPIDDigitizer::Sdi2Dig(&sdiLst,&digLst);
 153   AliHMPIDReconstructor::Dig2Clu(&digLst,&cluLst);   AliHMPIDReconstructor::Dig2Clu(&digLst,&cogLst,0);
 154
 155
 156   Int_t iDigN=0,iCluN=0,iCogN=0; for(Int_t i=0;i<7;i++){ iDigN+=((TClonesArray*)digLst.At(i))->GetEntries();
 157                                                          iCluN+=((TClonesArray*)cluLst.At(i))->GetEntries();
 158                                                          iCogN+=((TClonesArray*)cogLst.At(i))->GetEntries(); }
 159
 160   Printf("SDI------SDI---------SDI--------SDI------SDI------SDI #%i",sdiLst.GetEntries());sdiLst.Print();Printf("");
 161   Printf("DIG------DIG---------DIG--------DIG------DIG------DIG #%i",iDigN              );digLst.Print();Printf("");
 162   Printf("HIT------HIT---------HIT--------HIT------HIT------HIT #%i",hitLst.GetEntries());hitLst.Print();Printf("");
 163   Printf("CLU------CLU---------CLU--------CLU------CLU------CLU #%i",iCluN              );cluLst.Print();Printf("");
 164   Printf("COG------COG---------COG--------COG------COG------COG #%i",iCogN              );cogLst.Print();Printf("");
 165 }
 166
 167
 168
 169 void PrintMap()
 170 {
 171
 172   Double_t r2d=TMath::RadToDeg();
 173
 174   Double_t x=AliHMPIDDigit::SizeAllX(),y=AliHMPIDDigit::SizeAllY();
 175
 176   Printf("\n\n\n");
 177
 178   for(int ch=6;ch>=0;ch--){
 179     AliHMPIDDigit dL(ch,0,1,0,0),dR(ch,0,1,67,0);
 180     TVector3 lt=rp->Lors2Mars(ch,0,y);                                              TVector3 rt=rp->Lors2Mars(ch,x,y);
 181                                        TVector3 ce=rp->Lors2Mars(ch,x/2,y/2);
 182     TVector3 lb=rp->Lors2Mars(ch,0,0);                                              TVector3 rb=rp->Lors2Mars(ch,x,0);
 183
 184     Printf(" ____________________________");
 185     Printf("|%6.2fcm            %6.2fcm|"         ,lt.Mag()                             , rt.Mag()       );
 186     Printf("|%6.2fde            %6.2fde|"         ,lt.Theta()*r2d                       , rt.Theta()*r2d );
 187     Printf("|%6.2fde            %6.2fde|"         ,lt.Phi()*r2d                         , rt.Phi()*r2d   );
 188     Printf("|                            |"                                                       );
 189     Printf("|DDL %2i    %7.2fcm   DDL %2i|"       ,dL.DdlIdx()    ,  ce.Mag()           , dR.DdlIdx()    );
 190     Printf("| 0x%x    %7.2fdeg   0x%x|"           ,dL.DdlId()     ,  ce.Theta()*r2d     , dR.DdlId()     );
 191     Printf("|          %7.2fdeg        |"                         ,  ce.Phi()*r2d                        );
 192     Printf("|                            |");
 193     Printf("|%6.2fcm            %6.2fcm|"         ,lb.Mag()                             , rb.Mag()       );
 194     Printf("|%6.2fde            %6.2fde|"         ,lb.Theta()*r2d                       , rb.Theta()*r2d );
 195     Printf("|%6.2fde     Ch%i    %6.2fde|"        ,lb.Phi()*r2d   ,  ch                 , rb.Phi()*r2d   );
 196     Printf(" ----------------------------");
 197   }
 198
 199   Double_t m[3];
 200   for(int i=0;i<1000;i++){
 201     Float_t xout=0,xin=gRandom->Rndm()*130.60;
 202     Float_t yout=0,yin=gRandom->Rndm()*126.16;
 203     Int_t   c=gRandom->Rndm()*6;
 204     rp->Lors2Mars(c,xin,yin,m);
 205     rp->Mars2Lors(c,m,xout,yout);
 206     if( (xin-xout) != 0) Printf("Problem in X");
 207     if( (yin-yout) != 0) Printf("Problem in Y");
 208   }
 209
 210   Int_t ddl,r,d,a,ch,raw,pc,px,py; AliHMPIDDigit dd;
 211
 212   ddl=0;raw=0x2214000;r= 8;d=8;a=20;
 213   ddl=1;raw=0x2214000;r= 8;d=8;a=20;
 214
 215
 216   ddl=2;raw=0x08d6000;r= 2;d=3;a=22;
 217   ddl=3;raw=0x08d6000;r= 2;d=3;a=22;
 218
 219
 220   ddl=6;raw=0x592e000;r=22;d=4;a=46;ch=3;pc=4;px=55;py=5;dd.Raw(ddl,raw);
 221   Printf("(ch=%i,pc=%i,x=%2i,y=%2i) ddl=%i raw=0x%h (r=%2i,d=%2i,a=%2i)",
 222            ch,   pc,  px,   py,     ddl,   raw,      r,    d,    a); dd.Print();
 223   ddl=7;raw=0x592e000;r=22;d=4;a=46;ch=3;pc=1;
 224
 225
 226 }//PrintMap()
 227
 228
 229 void rec()
 230 {
 231
 232
 233   Double_t ckovMax=0.75,ckovSim;
 234   Int_t nSim=0;
 235   while(nSim<3){
 236     ckovSim=gRandom->Rndm()*ckovMax;//0.6468;
 237     nSim=20*TMath::Power(TMath::Sin(ckovSim)/TMath::Sin(ckovMax),2); //scale number of photons
 238   }
 239
 240
 241   TClonesArray *pCluLst=new TClonesArray("AliHMPIDCluster");
 242   TPolyMarker  *pMipMap=new TPolyMarker();   pMipMap->SetMarkerStyle(8);  pMipMap->SetMarkerColor(kRed);
 243   TPolyMarker  *pPhoMap=new TPolyMarker();   pPhoMap->SetMarkerStyle(4);  pPhoMap->SetMarkerColor(kRed);
 244   TPolyMarker  *pBkgMap=new TPolyMarker();   pBkgMap->SetMarkerStyle(25); pBkgMap->SetMarkerColor(kRed);
 245   TPolyLine    *pRing  =new TPolyLine;                                    pRing->SetLineColor(kGreen);
 246   TPolyLine    *pOld   =new TPolyLine;                                    pOld->SetLineColor(kBlue);
 247
 248   Int_t iCluCnt=0; pMipMap->SetPoint(iCluCnt,x,y); new((*pCluLst)[iCluCnt++])   AliHMPIDCluster(1,x,y,200); //add mip cluster
 249
 250 //  for(int j=0;j<30;j++){                                                                                   //add bkg photons
 251 //    Float_t bkgX=gRandom->Rndm()*AliHMPIDDigit::SizeAllX();
 252 //    Float_t bkgY=gRandom->Rndm()*AliHMPIDDigit::SizeAllY();
 253 //    pBkgMap->SetPoint(iCluCnt,bkgX,bkgY); new((*pCluLst)[iCluCnt++]) AliHMPIDCluster(1,bkgX,bkgY,35);
 254 //  }
 255
 256
 257
 258
 259   AliHMPIDRecon    rec; rec.SetTrack(th,ph,x,y);
 260
 261   TVector2 pos;
 262   for(int i=0;i<nSim;i++){
 263     rec.TracePhot(ckovSim,gRandom->Rndm()*2*TMath::Pi(),pos);                                   //add photons
 264     if(AliHMPIDDigit::IsInDead(pos.X(),pos.Y())) continue;
 265     pPhoMap->SetPoint(iCluCnt,pos.X(),pos.Y());    new((*pCluLst)[iCluCnt++]) AliHMPIDCluster(1,pos.X(),pos.Y(),35);
 266   }
 267
 268
 269   Int_t nRec=0,nOld=0;
 270   Double_t ckovRec=rec.CkovAngle(pCluLst,nRec); Double_t err=TMath::Sqrt(rec.CkovSigma2());
 271   Double_t ckovOld=old.CkovAngle(pCluLst,nOld);
 272
 273   Printf("---------------- Now reconstructed --------------------");
 274
 275
 276   for(int j=0;j<100;j++){rec.TracePhot(ckovRec,j*0.0628,pos); pRing->SetPoint(j,pos.X(),pos.Y());}
 277   for(int j=0;j<100;j++){rec.TracePhot(ckovOld,j*0.0628,pos); pOld->SetPoint(j,pos.X(),pos.Y());}
 278
 279   new TCanvas;  AliHMPIDDigit::DrawPc();  pMipMap->Draw(); pPhoMap->Draw(); pBkgMap->Draw(); pRing->Draw();  pOld->Draw();
 280
 281   TLatex txt; txt.SetTextSize(0.03);
 282   txt.DrawLatex(65,127,Form("#theta=%.4f#pm%.5f with %2i #check{C}"                             ,ckovSim, 0.,nSim             ));
 283   txt.DrawLatex(65,122,Form("#theta=%.4f#pm%.5f with %2i #check{C} Old=%.4f with %i #check{C}"  ,ckovRec,err,nRec,ckovOld,nOld));
 284   txt.DrawLatex(0 ,127,Form("#theta=%.2f#circ   #phi=%.2f#circ @(%.2f,%.2f) ",th*TMath::RadToDeg(),ph*TMath::RadToDeg(),x,y));
 285
 286 //  for(int i=0;i<35;i++){
 287 //    Double_t ckov=0.1+i*0.02;
 288 //    Printf("Ckov=%.2f Old=%.3f New=%.3f",ckov,old.FindRingArea(ckov),rec.FindRingArea(ckov));
 289 //  }
 290
 291
 292 }//rec()
 293
 294
 295 void HitQA(Double_t cut,Double_t cutele,Double_t cutR)
 296 {
 297 // Provides a set of control plots intended primarily for charged particle flux analisys
 298 // Arguments: cut (GeV)    - cut on momentum of any charged particles but electrons,
 299 //            cetele (GeV) - the same for electrons-positrons
 300 //            cutR (cm)    - cut on production vertex radius (cylindrical system)
 301   gBenchmark->Start("HitsAna");
 302
 303   Double_t cutPantiproton    =cut;
 304   Double_t cutPkaonminus     =cut;
 305   Double_t cutPpionminus     =cut;
 306   Double_t cutPmuonminus     =cut;
 307   Double_t cutPpositron      =cutele;
 308
 309   Double_t cutPelectron      =cutele;
 310   Double_t cutPmuonplus      =cut;
 311   Double_t cutPpionplus      =cut;
 312   Double_t cutPkaonplus      =cut;
 313   Double_t cutPproton        =cut;
 314
 315
 316   TH2F *pEleHitRZ    =new TH2F("EleHitRZ"    ,Form("e^{+} e^{-} hit %s;z[cm];R[cm]" ,GetName())     , 400,-300,300 ,400,-500,500);   //R-z plot 0cm<R<550cm -300cm<z<300cm
 317   TH2F *pEleHitRP    =new TH2F("EleHitRP"    ,Form("e^{+} e^{-} hit %s;p[GeV];R[cm]",GetName())     ,1000,-1  ,1   ,400,   0,550);   //R-p plot 0cm<R<550cm -1GeV<p<1GeV
 318   TH1F *pEleAllP     =new TH1F("EleAllP"     ,     "e^{+} e^{-} all;p[GeV]"                         ,1000,-1  ,1                );
 319   TH1F *pEleHitP     =new TH1F("EleHitP"     ,Form("e^{+} e^{-} hit %s;p[GeV]"      ,GetName())     ,1000,-1  ,1                );
 320   TH1F *pMuoHitP     =new TH1F("MuoHitP"     ,Form("#mu^{-} #mu^{+} hit %s;p[GeV]"  ,GetName())     ,1000,-4  ,4                );
 321   TH1F *pPioHitP     =new TH1F("PioHitP"     ,Form("#pi^{-} #pi^{+} hit %s;p[GeV]"  ,GetName())     ,1000,-4  ,4                );
 322   TH1F *pKaoHitP     =new TH1F("KaoHitP"     ,Form("K^{-} K^{+} hit %s;p[GeV]"      ,GetName())     ,1000,-4  ,4                );
 323   TH1F *pProHitP     =new TH1F("ProHitP"     ,Form("p^{-} p^{+} hit %s;p[GeV]"      ,GetName())     ,1000,-4  ,4                );
 324   TH2F *pFlux        =new TH2F("flux"        ,Form("%s flux with Rvertex<%.1fcm"    ,GetName(),cutR),10  ,-5  ,5   , 10,0    ,10); //special text hist
 325   TH2F *pVertex      =new TH2F("vertex"      ,Form("%s 2D vertex of HMPID hit;x;y"   ,GetName())     ,120 ,0   ,600 ,120,0    ,600); //special text hist
 326   TH1F *pRho         =new TH1F("rho"         ,Form("%s r of HMPID hit"               ,GetName())     ,600 ,0   ,600); //special text hist
 327   pFlux->SetStats(0);
 328   pFlux->GetXaxis()->SetBinLabel(1 ,Form("p^{-}>%.3fGeV/c"   ,cutPantiproton));
 329   pFlux->GetXaxis()->SetBinLabel(2 ,Form("K^{-}>%.3fGeV/c"   ,cutPkaonminus ));
 330   pFlux->GetXaxis()->SetBinLabel(3 ,Form("#pi^{-}>%.3fGeV/c" ,cutPpionminus ));
 331   pFlux->GetXaxis()->SetBinLabel(4 ,Form("#mu^{-}>%.3fGeV/c" ,cutPmuonminus ));
 332   pFlux->GetXaxis()->SetBinLabel(5 ,Form("e^{+}>%.3fGeV/c"   ,cutPpositron  ));
 333
 334   pFlux->GetXaxis()->SetBinLabel(6 ,Form("e^{-}>%.3fGeV/c"   ,cutPelectron  ));
 335   pFlux->GetXaxis()->SetBinLabel(7 ,Form("#mu^{+}>%.3fGeV/c" ,cutPmuonplus  ));
 336   pFlux->GetXaxis()->SetBinLabel(8 ,Form("#pi^{+}>%.3fGeV/c" ,cutPpionplus  ));
 337   pFlux->GetXaxis()->SetBinLabel(9 ,Form("K^{+}>%.3fGeV/c"   ,cutPkaonplus  ));
 338   pFlux->GetXaxis()->SetBinLabel(10,Form("p^{+}>%.3fGeV/c"   ,cutPproton    ));
 339
 340   pFlux->GetYaxis()->SetBinLabel(1,"sum");
 341   pFlux->GetYaxis()->SetBinLabel(2,"ch1");
 342   pFlux->GetYaxis()->SetBinLabel(3,"ch2");
 343   pFlux->GetYaxis()->SetBinLabel(4,"ch3");
 344   pFlux->GetYaxis()->SetBinLabel(5,"ch4");
 345   pFlux->GetYaxis()->SetBinLabel(6,"ch5");
 346   pFlux->GetYaxis()->SetBinLabel(7,"ch6");
 347   pFlux->GetYaxis()->SetBinLabel(8,"ch7");
 348   pFlux->GetYaxis()->SetBinLabel(9,"prim");
 349   pFlux->GetYaxis()->SetBinLabel(10,"tot");
 350
 351 //end of hists definition
 352
 353   Int_t iNevents=fLoader->GetRunLoader()->GetAliRun()->GetEventsPerRun(),iCntPrimParts=0,iCntTotParts=0;
 354 //load all needed trees
 355   fLoader->LoadHits();
 356   fLoader->GetRunLoader()->LoadHeader();
 357   fLoader->GetRunLoader()->LoadKinematics();
 358
 359   for(Int_t iEvtN=0;iEvtN < iNevents;iEvtN++){//events loop
 360     fLoader->GetRunLoader()->GetEvent(iEvtN);
 361     AliInfo(Form(" %i event processes",fLoader->GetRunLoader()->GetEventNumber()));
 362     AliStack *pStack= fLoader->GetRunLoader()->Stack();
 363
 364     for(Int_t iParticleN=0;iParticleN<pStack->GetNtrack();iParticleN++){//stack loop
 365       TParticle *pPart=pStack->Particle(iParticleN);
 366
 367       if(iParticleN%10000==0) AliInfo(Form(" %i particles read",iParticleN));
 368
 369       switch(pPart->GetPdgCode()){
 370         case kProtonBar: pFlux->Fill(-4.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-4.5,8); break;
 371         case kKMinus:    pFlux->Fill(-3.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-3.5,8); break;
 372         case kPiMinus:   pFlux->Fill(-2.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-2.5,8); break;
 373         case kMuonMinus: pFlux->Fill(-1.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-1.5,8); break;
 374         case kPositron:  pFlux->Fill(-0.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-0.5,8); pEleAllP->Fill(-pPart->P()); break;
 375
 376         case kElectron:  pFlux->Fill( 0.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 0.5,8); pEleAllP->Fill( pPart->P()); break;
 377         case kMuonPlus:  pFlux->Fill( 1.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 1.5,8); break;
 378         case kPiPlus:    pFlux->Fill( 2.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 2.5,8); break;
 379         case kKPlus:     pFlux->Fill( 3.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 3.5,8); break;
 380         case kProton:    pFlux->Fill( 4.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 4.5,8); break;
 381       }//switch
 382     }//stack loop
 383 //now hits analiser
 384     for(Int_t iEntryN=0;iEntryN < fLoader->TreeH()->GetEntries();iEntryN++){//TreeH loop
 385       fLoader->TreeH()->GetEntry(iEntryN);                                  //get current entry (prim)
 386       for(Int_t iHitN=0;iHitN < Hits()->GetEntries();iHitN++){//hits loop
 387         AliHMPIDHit *pHit = (AliHMPIDHit*)Hits()->At(iHitN);            //get current hit
 388         TParticle  *pPart=pStack->Particle(pHit->GetTrack());      //get stack particle which produced the current hit
 389
 390         if(pPart->GetPDG()->Charge()!=0&&pPart->Rho()>0.1) pVertex->Fill(pPart->Vx(),pPart->Vy()); //safe margin for sec.
 391         if(pPart->GetPDG()->Charge()!=0) pRho->Fill(pPart->Rho()); //safe margin for sec.
 392         if(pPart->R()>cutR) continue;                                   //cut on production radius (cylindrical system)
 393
 394         switch(pPart->GetPdgCode()){
 395           case kProtonBar: if(pPart->P()>cutPantiproton) {pProHitP->Fill(-pPart->P()); pFlux->Fill(-4.5,pHit->Ch());}break;
 396           case kKMinus   : if(pPart->P()>cutPkaonminus)  {pKaoHitP->Fill(-pPart->P()); pFlux->Fill(-3.5,pHit->Ch());}break;
 397           case kPiMinus  : if(pPart->P()>cutPpionminus)  {pPioHitP->Fill(-pPart->P()); pFlux->Fill(-2.5,pHit->Ch());}break;
 398           case kMuonMinus: if(pPart->P()>cutPmuonminus)  {pMuoHitP->Fill(-pPart->P()); pFlux->Fill(-1.5,pHit->Ch());}break;
 399           case kPositron : if(pPart->P()>cutPpositron)   {pEleHitP->Fill(-pPart->P()); pFlux->Fill(-0.5,pHit->Ch());
 400                pEleHitRP->Fill(-pPart->P(),pPart->R());  pEleHitRZ->Fill(pPart->Vz(),pPart->R()); }break;
 401
 402           case kElectron : if(pPart->P()>cutPelectron)   {pEleHitP->Fill( pPart->P()); pFlux->Fill( 0.5,pHit->Ch());
 403                pEleHitRP->Fill( pPart->P(),pPart->R());  pEleHitRZ->Fill(pPart->Vz(),pPart->R()); }break;
 404           case kMuonPlus : if(pPart->P()>cutPmuonplus)   {pMuoHitP->Fill( pPart->P()); pFlux->Fill( 1.5,pHit->Ch());}break;
 405           case kPiPlus   : if(pPart->P()>cutPpionplus)   {pPioHitP->Fill( pPart->P()); pFlux->Fill( 2.5,pHit->Ch());}break;
 406           case kKPlus    : if(pPart->P()>cutPkaonplus)   {pKaoHitP->Fill( pPart->P()); pFlux->Fill( 3.5,pHit->Ch());}break;
 407           case kProton   : if(pPart->P()>cutPproton)     {pProHitP->Fill( pPart->P()); pFlux->Fill( 4.5,pHit->Ch());}break;
 408         }
 409       }//hits loop
 410     }//TreeH loop
 411     iCntPrimParts +=pStack->GetNprimary();
 412     iCntTotParts  +=pStack->GetNtrack();
 413   }//events loop
 414 //unload all loaded staff
 415   fLoader->UnloadHits();
 416   fLoader->GetRunLoader()->UnloadHeader();
 417   fLoader->GetRunLoader()->UnloadKinematics();
 418 //Calculater some sums
 419   Stat_t sum=0;
 420 //sum row, sum over rows
 421   for(Int_t i=1;i<=pFlux->GetNbinsX();i++){
 422     sum=0; for(Int_t j=2;j<=8;j++)    sum+=pFlux->GetBinContent(i,j);
 423     pFlux->SetBinContent(i,1,sum);
 424   }
 425
 426 //display everything
 427   new TCanvas("canvas1",Form("Events %i Nprims=%i Nparticles=%i",iNevents,iCntPrimParts,iCntTotParts),1000,900); pFlux->Draw("text");  gPad->SetGrid();
 428 //total prims and particles
 429   TLatex latex; latex.SetTextSize(0.02);
 430   sum=0; for(Int_t i=1;i<=pFlux->GetNbinsX();i++) sum+=pFlux->GetBinContent(i,10);    latex.DrawLatex(5.1,9.5,Form("%.0f",sum));
 431   sum=0; for(Int_t i=1;i<=pFlux->GetNbinsX();i++) sum+=pFlux->GetBinContent(i,9);     latex.DrawLatex(5.1,8.5,Form("%.0f",sum));
 432   for(Int_t iCh=0;iCh<7;iCh++) {
 433     sum=0; for(Int_t i=1;i<=pFlux->GetNbinsX();i++) sum+=pFlux->GetBinContent(i,iCh+2);latex.DrawLatex(5.1,iCh+1.5,Form("%.0f",sum));
 434   }
 435   sum=0; for(Int_t i=1;i<=pFlux->GetNbinsX();i++) sum+=pFlux->GetBinContent(i,1);    latex.DrawLatex(5.1,0.5,Form("%.0f",sum));
 436
 437   new TCanvas("cEleAllP"   ,"e" ,200,100); pEleAllP->Draw();
 438   new TCanvas("cEleHitRP"  ,"e" ,200,100); pEleHitRP->Draw();
 439   new TCanvas("cEleHitRZ"  ,"e" ,200,100); pEleHitRZ->Draw();
 440   new TCanvas("cEleHitP"   ,"e" ,200,100); pEleHitP->Draw();
 441   new TCanvas("cMuoHitP"   ,"mu",200,100); pMuoHitP->Draw();
 442   new TCanvas("cPioHitP"   ,"pi",200,100); pPioHitP->Draw();
 443   new TCanvas("cKaoHitP"   ,"K" ,200,100); pKaoHitP->Draw();
 444   new TCanvas("cProHitP"   ,"p" ,200,100); pProHitP->Draw();
 445   new TCanvas("cVertex"    ,"2d vertex" ,200,100); pVertex->Draw();
 446   new TCanvas("cRho"    ,"Rho of sec" ,200,100); pRho->Draw();
 447
 448   gBenchmark->Show("HitsPlots");
 449 }//HitQA()
 450
 451
 452 void RecWithStack()
 453 {
 454   al->LoadHeader();al->LoadKinematics();
 455   AliStack *pStk=al->Stack();
 456
 457   AliESD *pEsd=new AliESD;
 458   for(Int_t iTrk=0;iTrk<pStk->GetNtrack();iTrk++){//stack loop
 459     TParticle *pPart=pStk->Particle(iTrk);
 460     if(pPart->GetPDG()->Charge()==0) continue; //neutral particles are not reconstructed
 461     pEsd->AddTrack(new AliESDtrack(pPart));
 462   }//stack loop
 463
 464   pEsd->Print();
 465   AliTracker::SetFieldMap(new AliMagF,1);
 466   AliHMPIDTracker t;
 467   rl->LoadRecPoints();
 468   t.LoadClusters(rl->TreeR());
 469   t.PropagateBack(pEsd);
 470   rl->UnloadRecPoints();
 471 }
 472
 473
 474 void AliHMPIDReconstructor::CluQA(AliRunLoader *pAL)
 475 {
 476 // Quality assesment plots for clusters.
 477 // This methode takes list of digits and form list of clusters again in order to
 478 // calculate cluster shape and cluster particle mixture
 479   AliLoader *pRL=pAL->GetDetectorLoader("HMPID");  AliHMPID *pRich=(AliHMPID*)pAL->GetAliRun()->GetDetector("HMPID");//get pointers for HMPID and HMPID loader
 480   Int_t iNevt=pAL->GetNumberOfEvents();  if(iNevt==0)             {AliInfoClass("No events");return;}
 481                                          if(pRL->LoadDigits())    {AliInfoClass("No digits file");return;}
 482                                             pAL->LoadHeader();
 483                                             pAL->LoadKinematics();
 484 //                                            AliStack *pStack=pAL->Stack();
 485   TH1::AddDirectory(kFALSE);
 486
 487
 488   TH1F*    pQ=new TH1F("RiAllQ"  ,"Charge All"           ,4000 ,0  ,4000);// Q hists
 489   TH1F* pCerQ=new TH1F("RiCerQ"  ,"Charge Ckov"          ,4000 ,0  ,4000);
 490   TH1F* pMipQ=new TH1F("RiMipQ"  ,"Charge MIP"           ,4000 ,0  ,4000);
 491
 492   TH1F*    pS=new TH1F("RichCluSize"    ,"Cluster size;size"         ,100  ,0  ,100 );// size hists
 493   TH1F* pCerS=new TH1F("RichCluCerSize" ,"Ckov size;size"            ,100  ,0  ,100 );
 494   TH1F* pMipS=new TH1F("RichCluMipSize" ,"MIP size;size"             ,100  ,0  ,100 );
 495
 496   TH2F*    pM=new TH2F("RichCluMap"     ,"Cluster map;x [cm];y [cm]" ,1000 ,0  ,AliHMPIDDigit::SizePcX(),1000,0,AliHMPIDDigit::SizePcY()); // maps
 497   TH2F* pMipM=new TH2F("RichCluMipMap"  ,"MIP map;x [cm];y [cm]"     ,1000 ,0  ,AliHMPIDDigit::SizePcX(),1000,0,AliHMPIDDigit::SizePcY());
 498   TH2F* pCerM=new TH2F("RichCluCerMap"  ,"Ckov map;x [cm];y [cm]"    ,1000 ,0  ,AliHMPIDDigit::SizePcX(),1000,0,AliHMPIDDigit::SizePcY());
 499
 500
 501
 502   for(Int_t iEvt=0;iEvt<iNevt;iEvt++){
 503     pAL->GetEvent(iEvt);
 504     pRL->TreeD()->GetEntry(0);
 505     TClonesArray *pCluLst=new TClonesArray("AliHMPIDCluster");//tmp list of clusters for this event
 506
 507     for(Int_t iCh=0;iCh<7;iCh++) Dig2Clu(pRich->DigLst(iCh),pCluLst,kFALSE);//cluster finder for all chamber if any digits present
 508
 509     for(Int_t iClu=0;iClu<pCluLst->GetEntriesFast();iClu++){
 510       AliHMPIDCluster *pClu = (AliHMPIDCluster*)pCluLst->At(iClu);
 511       Int_t cfm=0; for(Int_t iDig=0;iDig<pClu->Size();iDig++)  cfm+=pClu->Dig(iDig)->Ch(); //collect ckov-fee-mip structure of current cluster ?????
 512       Int_t iNckov=cfm/1000000;      Int_t iNfee =cfm%1000000/1000;      Int_t iNmip =cfm%1000000%1000;
 513
 514                                              pQ   ->Fill(pClu->Q()) ; pS   ->Fill(pClu->Size()) ; pM    ->Fill(pClu->X(),pClu->Y()); //all clusters
 515       if(iNckov!=0 && iNfee==0 && iNmip==0) {pCerQ->Fill(pClu->Q()) ; pCerS->Fill(pClu->Size()) ; pCerM ->Fill(pClu->X(),pClu->Y());}//ckov only cluster
 516       if(iNckov==0 && iNfee==0 && iNmip!=0) {pMipQ->Fill(pClu->Q()) ; pMipS->Fill(pClu->Size()) ; pMipM ->Fill(pClu->X(),pClu->Y());}//mip only cluster
 517
 518     }//clusters loop
 519     pCluLst->Clear();delete pCluLst;
 520   }//events loop
 521
 522   pRL->UnloadDigits(); pAL->UnloadKinematics(); pAL->UnloadHeader();
 523   TCanvas *pC=new TCanvas("RichCluQA",Form("QA for cluster from %i events",iNevt),1000,900); pC->Divide(3,3);
 524   pC->cd(1);    pM->Draw();          pC->cd(2);    pQ->Draw();       pC->cd(3);    pS->Draw();
 525   pC->cd(4); pMipM->Draw();          pC->cd(5); pMipQ->Draw();       pC->cd(6); pMipS->Draw();
 526   pC->cd(7); pCerM->Draw();          pC->cd(8); pCerQ->Draw();       pC->cd(9); pCerS->Draw();
 527 }//CluQA()
 528
 529
 530
 531 void AliHMPID::OccupancyPrint(Int_t iEvtNreq)
 532 {
 533 //prints occupancy for each chamber in a given event
 534   Int_t iEvtNmin,iEvtNmax;
 535   if(iEvtNreq==-1){
 536     iEvtNmin=0;
 537     iEvtNmax=gAlice->GetEventsPerRun();
 538   } else {
 539     iEvtNmin=iEvtNreq;iEvtNmax=iEvtNreq+1;
 540   }
 541
 542   if(GetLoader()->GetRunLoader()->LoadHeader()) return;
 543   if(GetLoader()->GetRunLoader()->LoadKinematics()) return;
 544
 545 //  Info("Occupancy","for event %i",iEvtN);
 546   if(GetLoader()->LoadHits()) return;
 547   if(GetLoader()->LoadDigits()) return;
 548
 549
 550   for(Int_t iEvtN=iEvtNmin;iEvtN<iEvtNmax;iEvtN++){
 551     Int_t nDigCh[7]={0,0,0,0,0,0,0};
 552     Int_t iChHits[7]={0,0,0,0,0,0,0};
 553     Int_t nPrim[7]={0,0,0,0,0,0,0};
 554     Int_t nSec[7]={0,0,0,0,0,0,0};
 555     AliInfo(Form("events processed %i",iEvtN));
 556     if(GetLoader()->GetRunLoader()->GetEvent(iEvtN)) return;
 557     AliStack *pStack = GetLoader()->GetRunLoader()->Stack();
 558     for(Int_t iPrimN=0;iPrimN<GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop
 559       GetLoader()->TreeH()->GetEntry(iPrimN);
 560       for(Int_t iHitN=0;iHitN<Hits()->GetEntries();iHitN++){
 561         AliHMPIDHit *pHit = (AliHMPIDHit*)Hits()->At(iHitN);
 562         if(pHit->E()>0){
 563           iChHits[pHit->Ch()]++;
 564           if(pStack->Particle(pHit->GetTrack())->Rho()<0.01) nPrim[pHit->Ch()]++;else nSec[pHit->Ch()]++;
 565         }
 566       }
 567     }
 568
 569     GetLoader()->TreeD()->GetEntry(0);
 570     for(Int_t iCh=0;iCh<7;iCh++){
 571       for(Int_t iDig=0;iDig<DigLst(iCh)->GetEntries();iDig++){
 572         AliHMPIDDigit *pDig=(AliHMPIDDigit*)DigLst(iCh)->At(iDig);
 573         nDigCh[pDig->Ch()]++;
 574       }
 575       Printf("Occupancy for chamber %i = %4.2f %% and charged prim tracks %i and sec. tracks %i with total %i",
 576         iCh,Float_t(nDigCh[iCh])*100/AliHMPIDDigit::kPadAll,nPrim[iCh],nSec[iCh],iChHits[iCh]);
 577     }
 578
 579
 580   }//events loop
 581   GetLoader()->UnloadHits();
 582   GetLoader()->UnloadDigits();
 583   GetLoader()->GetRunLoader()->UnloadHeader();
 584   GetLoader()->GetRunLoader()->UnloadKinematics();
 585 }
 586
 587
 588 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 589 void AliHMPID::SummaryOfEvent(Int_t iEvtN) const
 590 {
 591 //prints a summary for a given event
 592   AliInfo(Form("Summary of event %i",iEvtN));
 593   GetLoader()->GetRunLoader()->GetEvent(iEvtN);
 594   if(GetLoader()->GetRunLoader()->LoadHeader()) return;
 595   if(GetLoader()->GetRunLoader()->LoadKinematics()) return;
 596   AliStack *pStack=GetLoader()->GetRunLoader()->Stack();
 597
 598   AliGenEventHeader* pGenHeader =  gAlice->GetHeader()->GenEventHeader();
 599   if(pGenHeader->InheritsFrom("AliGenHijingEventHeader")) {
 600     AliInfo(Form(" Hijing event with impact parameter b = %.2f (fm)",((AliGenHijingEventHeader*) pGenHeader)->ImpactParameter()));
 601   }
 602   Int_t nChargedPrimaries=0;
 603   for(Int_t i=0;i<pStack->GetNtrack();i++) {
 604     TParticle *pParticle = pStack->Particle(i);
 605     if(pParticle->IsPrimary()&&pParticle->GetPDG()->Charge()!=0) nChargedPrimaries++;
 606     }
 607   AliInfo(Form("Total number of         primaries %i",pStack->GetNprimary()));
 608   AliInfo(Form("Total number of charged primaries %i",nChargedPrimaries));
 609   AliInfo(Form("Total n. of tracks in stack(+sec) %i",pStack->GetNtrack()));
 610   GetLoader()->GetRunLoader()->UnloadHeader();
 611   GetLoader()->GetRunLoader()->UnloadKinematics();
 612 }
 613 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++