MUON/AliMUONproto.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /*
  16 $Log$
  17 Revision 1.1.2.1  2000/04/18 09:11:15  morsch
  18 Implementation of MUON Chamber Prototype Class
  19 Both read digits from raw data or use the Monte-Carlo.
  20 Rachid GUERNANE, IPN Lyon guernane@ipnl.in2p3.fr
  21
  22 */
  23
  24 /*
  25 Implementation of MUON Chamber Prototype Class
  26 Both read digits from raw data or use the Monte-Carlo.
  27 1-Feb-2000 Rachid GUERNANE, IPN Lyon guernane@ipnl.in2p3.fr
  28 */
  29
  30
  31 ////////////////////////////////////////////////
  32 //  Manager and hits classes for set:PROTO    //
  33 ////////////////////////////////////////////////
  34
  35 #include <TTUBE.h>
  36 #include <TBRIK.h>
  37 #include <TRotMatrix.h>
  38 #include <TNode.h>
  39 #include <TTree.h>
  40 #include <TRandom.h>
  41 #include <TObject.h>
  42 #include <TVector.h>
  43 #include <TObjArray.h>
  44 #include <TMinuit.h>
  45 #include <TParticle.h>
  46 #include <TROOT.h>
  47 #include <TFile.h>
  48 #include <TNtuple.h>
  49 #include <TCanvas.h>
  50 #include <TPad.h>
  51 #include <TDirectory.h>
  52 #include <TObjectTable.h>
  53 #include <AliPDG.h>
  54 #include <TArrayI.h>
  55 #include <AliDetector.h>
  56
  57 #include "AliMUONproto.h"
  58 #include "TTUBE.h"
  59 #include "AliMUONClusterFinder.h"
  60 #include "AliRun.h"
  61 #include "AliMC.h"
  62 #include "iostream.h"
  63 #include "AliCallf77.h"
  64 #include "AliConst.h"
  65 #include "chainalice2.h"
  66 #include "AliMUONSegResV1.h"
  67 #include "AliMUONSegResV11.h"
  68
  69 ClassImp(AliMUONproto)
  70
  71 //___________________________________________
  72 AliMUONproto::AliMUONproto()
  73         : AliMUON()
  74 {
  75     cout << "\n Calling AliMUONproto constructor..." << endl;
  76
  77         //
  78         //
  79         //
  80 }
  81
  82 //___________________________________________
  83 AliMUONproto::AliMUONproto(const char *name, const char *title)
  84         : AliMUON(name,title)
  85 //  AliMUON defines everything, but the chamber for NCH=1
  86 {
  87
  88 //  z-Positions of Chambers
  89     const Float_t zch[1] = {975.};
  90 //
  91 //  inner diameter
  92     const Float_t dmi[1] = {0.};
  93 //
  94 //  outer diameter
  95     const Float_t dma[1] = {29.2};
  96 //
  97 //
  98 //    Default Parameters for ALICE2 prototype
  99
 100 //
 101     (*fChambers)[0] = new AliMUONChamber();
 102     AliMUONChamber* chamber = (AliMUONChamber*) (*fChambers)[0];
 103     chamber->SetGid(0);
 104     chamber->SetZPOS(zch[0]);
 105 //
 106     chamber->InitGeo(zch[0]);
 107     chamber->SetRInner(dmi[0]/2);
 108     chamber->SetROuter(dma[0]/2);
 109
 110     for (int i = 0; i <= 99; i++) {
 111         fThreshold[i] = 0.;
 112     }
 113
 114 }
 115
 116 //___________________________________________
 117 void AliMUONproto::GetRawDigits(Int_t evnb, Int_t *lptr, Int_t ilen) {
 118
 119     Int_t ip = 0;
 120     Int_t equip = 0;
 121     Int_t nochnl;
 122     Int_t loop;
 123     Int_t val;
 124     Int_t itype;
 125     Int_t id;
 126     Int_t serial;
 127     Int_t equiplen;
 128
 129     Int_t digits[5];
 130
 131     AliMUON* MUON  = (AliMUON*)gAlice->GetModule("MUON");
 132     AliMUONSegmentationV1* seg = (AliMUONSegmentationV1*) Chamber(0).GetSegmentationModel(1);
 133
 134     Int_t tracks[10];
 135     Int_t charges[10];
 136
 137     for (Int_t i = 0; i < 10; i++) {
 138         tracks[i] = 0;
 139         charges[i] = 0;
 140     }
 141
 142     Int_t ich = 0;
 143
 144   nwtype:
 145
 146     itype = lptr[ip++];
 147     id = lptr[ip++];
 148     equiplen = lptr[ip++];
 149
 150     if (equiplen < 0 ) return;
 151
 152     if (itype != (int)(0XCACCA008)) {
 153         ip += equiplen;
 154         if (ip < ilen-2)  goto nwtype;
 155     }
 156     else {
 157         serial = id >> 16;
 158         equip = id & 0x1;
 159         if ((serial == 190) && (equip == 1)) {
 160             for (loop = 0; loop < equiplen; loop++) {
 161                 nochnl = (lptr[ip] & 0x7ff000 ) >> 12;
 162                 val = lptr[ip] & 0x3ff;
 163                 // fill digits from raw data according to cathode connexions
 164                 if (group[nochnl][2][1]!=0)
 165                     digits[0] = group[nochnl][2][1] - 12;
 166                 else if  (group[nochnl][1][1]!=0)
 167                     digits[0] = group[nochnl][1][1] - 12;
 168                 else
 169                     digits[0] = group[nochnl][0][1] - 12;
 170                 digits[1] = group[nochnl][0][0];
 171                 if (digits[0] != seg->Ix(digits[0], digits[1]))
 172                     printf("Pb pour ix=%d,iy=%d\n", digits[0], digits[1]);
 173                 digits[2] = val;
 174                 digits[3] = 0;
 175                 digits[4] = 0;
 176                 if (digits[2] >= fThreshold[nochnl])
 177                     MUON->AddDigits(ich, tracks, charges, digits);
 178
 179                 ip++;
 180             }
 181         }
 182         else
 183             ip += equiplen;
 184
 185         if (ip < ilen-2) goto nwtype;
 186     }
 187
 188     gAlice->TreeD()->Fill();
 189     MUON->ResetDigits();
 190
 191     gAlice->TreeD()->Fill();
 192     MUON->ResetDigits();
 193
 194     char hname[30];
 195     sprintf(hname, "TreeD%d", evnb);
 196     gAlice->TreeD()->Write(hname);
 197     // reset tree
 198     gAlice->TreeD()->Reset();
 199
 200 }
 201
 202 //___________________________________________
 203 void AliMUONproto::SetPadSize(Int_t id, Int_t isec, Float_t p1, Float_t p2)
 204 {
 205     Int_t i=2*(id-1);
 206     ((AliMUONChamber*)(*fChambers)[i])->SetPadSize(isec,p1,p2);
 207 }
 208
 209 //___________________________________________
 210 void AliMUONproto::SetChargeSlope(Int_t id, Float_t p1)
 211 {
 212     Int_t i=2*(id-1);
 213     ((AliMUONChamber*) (*fChambers)[i])->SetChargeSlope(p1);
 214 }
 215
 216 //___________________________________________
 217 void AliMUONproto::SetChargeSpread(Int_t id, Float_t p1, Float_t p2)
 218 {
 219     Int_t i=2*(id-1);
 220     ((AliMUONChamber*) (*fChambers)[i])->SetChargeSpread(p1,p2);
 221 }
 222
 223 //___________________________________________
 224 void AliMUONproto::SetMaxAdc(Int_t id, Float_t p1)
 225 {
 226     Int_t i=2*(id-1);
 227     ((AliMUONChamber*) (*fChambers)[i])->SetMaxAdc(p1);
 228 }
 229
 230 //___________________________________________
 231 void AliMUONproto::CreateGeometry()
 232 {
 233     Int_t *idtmed = fIdtmed->GetArray()-1099;
 234 //
 235     Float_t zpos;
 236     Float_t bpar[3];
 237     Int_t idGas=idtmed[1100];
 238
 239     AliMUONChamber *iChamber;
 240 //********************************************************************
 241 //                          Prototype ALICE2                        **
 242 //********************************************************************
 243     iChamber=(AliMUONChamber*) (*fChambers)[0];
 244
 245     zpos=iChamber->ZPosition();
 246 //
 247     const Float_t X_POS = 11*.975/2; //half x
 248     const Float_t Y_POS = 18*.55/2;  //half y
 249     const Float_t Z_POS = 0.325;
 250
 251     bpar[0] = X_POS;
 252     bpar[1] = Y_POS;
 253     bpar[2] = Z_POS;
 254
 255     gMC->Gsvolu("C01G", "BOX", idGas, bpar, 3);
 256     gMC->Gspos("C01G", 1, "ALIC", -bpar[0], bpar[1], zpos, 0, "ONLY");
 257 }
 258
 259 //___________________________________________
 260 void AliMUONproto::CreateMaterials()
 261 {
 262   // *** DEFINITION OF AVAILABLE MUON MATERIALS ***
 263   //
 264   //     Ar-CO2 gas
 265     Float_t ag1[3]   = { 39.95,12.01,16. };
 266     Float_t zg1[3]   = { 18.,6.,8. };
 267     Float_t wg1[3]   = { .8,.0667,.13333 };
 268     Float_t dg1      = .001821;
 269
 270     Float_t epsil  = .001; // Tracking precision,
 271     Float_t tmaxfd = -20.; // Maximum angle due to field deflection
 272     Float_t stmin  = -.8;
 273     Int_t ISXFLD   = gAlice->Field()->Integ();
 274     Float_t SXMGMX = gAlice->Field()->Max();
 275
 276     //
 277     // --- Define the various materials for GEANT ---
 278     AliMixture(22, "ArCO2 80%$", ag1, zg1, dg1, 3, wg1);
 279
 280     //
 281     //    Ar/CO2
 282     AliMedium(1, "ARG_CO2   ", 22, 1, ISXFLD, SXMGMX, tmaxfd, fMaxStepGas,
 283               fMaxDestepAlu, epsil, stmin);
 284     //    Air
 285         //AliMedium(1, "AIR_CH_US         ", 15, 1, ISXFLD, SXMGMX, tmaxfd, -1., -.3, epsil, stmin);
 286 }
 287
 288 //___________________________________________
 289 void AliMUONproto::Init()
 290 {
 291    printf("\n\n\n Start Init for Prototype ALICE2 - CPC chamber type\n\n\n");
 292
 293    //
 294    // Initialize Tracking Chambers
 295    //
 296
 297    ( (AliMUONChamber*) (*fChambers)[0])->Init();
 298
 299    //
 300    // Set the chamber (sensitive region) GEANT identifier
 301    AliMC* gMC = AliMC::GetMC();
 302    ((AliMUONChamber*)(*fChambers)[0])->SetGid(gMC->VolId("C01G"));
 303
 304    printf("\n\n\n Finished Init for Prototype ALICE2 - CPC chamber type\n\n\n");
 305
 306 }
 307
 308 //___________________________________________
 309 void AliMUONproto::StepManager()
 310 {
 311   Int_t          copy, id;
 312   static Int_t   idvol;
 313   static Int_t   vol[2];
 314   Int_t          ipart;
 315   TLorentzVector pos;
 316   TLorentzVector mom;
 317   Float_t        theta,phi;
 318   Float_t        destep, step;
 319
 320   static Float_t eloss, eloss2, xhit, yhit, tlength;
 321   const  Float_t big=1.e10;
 322
 323   //  modifs perso
 324   static Float_t hits[14];
 325
 326   TClonesArray &lhits = *fHits;
 327
 328   //
 329   // Set maximum step size for gas
 330   // numed=gMC->GetMedium();
 331   //
 332   // Only charged tracks
 333   if( !(gMC->TrackCharge()) ) return;
 334   //
 335   // Only gas gap inside chamber
 336   // Tag chambers and record hits when track enters
 337   idvol=-1;
 338   id=gMC->CurrentVolID(copy);
 339
 340     for (Int_t i=1; i<=NCH; i++) {
 341       if(id==((AliMUONChamber*)(*fChambers)[i-1])->GetGid()){
 342           vol[0]=i;
 343           idvol=i-1;
 344       }
 345     }
 346     if (idvol == -1) return;
 347   //
 348   // Get current particle id (ipart), track position (pos)  and momentum (mom)
 349   gMC->TrackPosition(pos);
 350   gMC->TrackMomentum(mom);
 351
 352   ipart  = gMC->TrackPid();
 353   //Int_t ipart1 = gMC->IdFromPDG(ipart);
 354   //printf("ich, ipart %d %d \n",vol[0],ipart1);
 355
 356   //
 357   // momentum loss and steplength in last step
 358   destep = gMC->Edep();
 359   step   = gMC->TrackStep();
 360
 361   //
 362   // record hits when track enters ...
 363   if( gMC->IsTrackEntering()) {
 364       gMC->SetMaxStep(fMaxStepGas);
 365       Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
 366       Double_t rt = TMath::Sqrt(tc);
 367       Double_t pmom = TMath::Sqrt(tc+mom[2]*mom[2]);
 368       Double_t tx=mom[0]/pmom;
 369       Double_t ty=mom[1]/pmom;
 370       Double_t tz=mom[2]/pmom;
 371       Double_t s=((AliMUONChamber*)(*fChambers)[idvol])
 372           ->GetResponseModel()
 373           ->Pitch()/tz;
 374       theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
 375       phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
 376       hits[0] = Float_t(ipart);         // Geant3 particle type
 377       hits[1] = pos[0]+s*tx;                 // X-position for hit
 378       hits[2] = pos[1]+s*ty;                 // Y-position for hit
 379       hits[3] = pos[2]+s*tz;                 // Z-position for hit
 380       hits[4] = theta;                  // theta angle of incidence
 381       hits[5] = phi;                    // phi angle of incidence
 382       hits[8] = (Float_t) fNclusters;   // first padhit
 383       hits[9] = -1;                     // last pad hit
 384
 385       // modifs perso
 386       hits[10] = mom[3]; // hit momentum P
 387       hits[11] = mom[0]; // Px/P
 388       hits[12] = mom[1]; // Py/P
 389       hits[13] = mom[2]; // Pz/P
 390       // fin modifs perso
 391
 392       // phi angle of incidence
 393       tlength = 0;
 394       eloss   = 0;
 395       eloss2  = 0;
 396       xhit    = pos[0];
 397       yhit    = pos[1];
 398       // Only if not trigger chamber
 399       if(idvol<10) {
 400           //
 401           //  Initialize hit position (cursor) in the segmentation model
 402           ((AliMUONChamber*) (*fChambers)[idvol])
 403               ->SigGenInit(pos[0], pos[1], pos[2]);
 404       } else {
 405           //geant3->Gpcxyz();
 406           //printf("In the Trigger Chamber #%d\n",idvol-9);
 407       }
 408   }
 409   eloss2+=destep;
 410
 411   //
 412   // Calculate the charge induced on a pad (disintegration) in case
 413   //
 414   // Mip left chamber ...
 415   if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
 416       gMC->SetMaxStep(big);
 417       eloss   += destep;
 418       tlength += step;
 419
 420       // Only if not trigger chamber
 421       if(idvol<10) {
 422           if (eloss > 0)
 423               MakePadHits(0.5*(xhit+pos[0]),0.5*(yhit+pos[1]),eloss,idvol);
 424       }
 425
 426       hits[6]=tlength;
 427       hits[7]=eloss2;
 428       if (fNclusters > (Int_t)hits[8]) {
 429           hits[8]= hits[8]+1;
 430           hits[9]= (Float_t) fNclusters;
 431       }
 432
 433       new(lhits[fNhits++])
 434           AliMUONHit(fIshunt,gAlice->CurrentTrack(),vol,hits);
 435       eloss = 0;
 436       //
 437       // Check additional signal generation conditions
 438       // defined by the segmentation
 439       // model (boundary crossing conditions)
 440   } else if
 441       (((AliMUONChamber*) (*fChambers)[idvol])
 442        ->SigGenCond(pos[0], pos[1], pos[2]))
 443   {
 444       ((AliMUONChamber*) (*fChambers)[idvol])
 445           ->SigGenInit(pos[0], pos[1], pos[2]);
 446 //      printf("\n-> MakePadHits, reason special %d",ipart);
 447       if (eloss > 0)
 448           MakePadHits(0.5*(xhit+pos[0]),0.5*(yhit+pos[1]),eloss,idvol);
 449       xhit     = pos[0];
 450       yhit     = pos[1];
 451       eloss    = destep;
 452       tlength += step ;
 453       //
 454       // nothing special  happened, add up energy loss
 455   } else {
 456       eloss   += destep;
 457       tlength += step ;
 458   }
 459 }
 460
 461 //___________________________________________
 462 void AliMUONproto::BuildGeometry()
 463 {
 464     TNode *Node;
 465     TNode *Top;
 466
 467     const int kColorMUON = kBlue;
 468     //
 469     Top=gAlice->GetGeometry()->GetNode("alice");
 470     //
 471     //
 472     //
 473     float dx, dy, dz, zpos;
 474
 475     const Float_t cz[1]={975.};
 476
 477     zpos=cz[0];
 478
 479     dx=11*.975/2;
 480     dy=18*.55/2;
 481     dz=0.325;
 482
 483     new TBRIK("C_MUON101","Mother volume for Proto.","void",dx*2,dy*2,dz*2);
 484     TBRIK* PROTO = new TBRIK("PROT", "Proto. Sens. Region","void",dx,dy,dz);
 485     Top->cd();
 486     Node = new TNode("MUON101","ChamberNode","C_MUON101",0,0,zpos,"");
 487     Node->SetLineColor(kColorMUON);
 488     Node->SetVisibility(0);
 489     fNodes->Add(Node);
 490     Node->cd();
 491     Node = new TNode("MUON201", "Proto. Sens. Region Node", PROTO, -dx, dy, dz);
 492     Node->SetLineColor(kColorMUON);
 493
 494 }
 495
 496 //_____________________________________________________________________________
 497 void AliMUONproto::FindClusters(Int_t nev, Int_t last_entry)
 498 {
 499
 500 //
 501 // Loop on chambers and on cathode planes
 502 //
 503   for (Int_t icat = 0; icat < 2; icat++) {
 504             gAlice->ResetDigits();
 505             gAlice->TreeD()->GetEvent(last_entry+icat);
 506
 507       for (Int_t ich = 0; ich < NTrackingCh; ich++) {
 508           AliMUONChamber* iChamber=(AliMUONChamber*) (*fChambers)[ich];
 509           TClonesArray *MUONdigits  = this->DigitsAddress(ich);
 510           if (MUONdigits == 0) continue;
 511           //
 512           // Get ready the current chamber stuff
 513           //
 514
 515           AliMUONResponse* response = iChamber->GetResponseModel();
 516           AliMUONSegmentation* seg = iChamber->GetSegmentationModel(icat+1);
 517           AliMUONClusterFinder* rec = iChamber->GetReconstructionModel();
 518
 519           if (seg) {
 520               rec->SetSegmentation(seg);
 521               rec->SetResponse(response);
 522               rec->SetDigits(MUONdigits);
 523               rec->SetChamber(ich);
 524               rec->FindRawClusters();
 525           }
 526
 527           TClonesArray *fRch;
 528           fRch=RawClustAddress(ich);
 529           fRch->Sort();
 530           // it seems to work
 531
 532       } // for ich
 533       // fill the tree
 534       gAlice->TreeR()->Fill();
 535
 536       ResetRawClusters();
 537   } // for icat
 538
 539   char hname[30];
 540   sprintf(hname,"TreeR%d",nev);
 541   gAlice->TreeR()->Write(hname);
 542   gAlice->TreeR()->Reset();
 543 }
 544
 545 //_____________________________________________________________________________
 546 void AliMUONproto::SetThreshold()
 547 {
 548
 549     ifstream inputFile("/home/alice/guernane/aliroot/pro/MUON/crped190.dat", ios::in);
 550
 551     if (inputFile.fail()) {
 552         cout << "Error opening file" << endl;
 553         exit(2);
 554     }
 555
 556     char buff[32];
 557     Int_t Serial;
 558     Int_t Ntrigger;
 559     Int_t Nchannel;
 560     Int_t i1;
 561     Int_t i2;
 562
 563     inputFile >> buff;
 564
 565     inputFile >> Serial;
 566     inputFile >> Ntrigger;
 567     inputFile >> Nchannel;
 568     inputFile >> i1;
 569     inputFile >> i2;
 570
 571     Float_t ped0[Nchannel];
 572     Float_t sig0[Nchannel];
 573     Float_t ped1[Nchannel];
 574     Float_t sig1[Nchannel];
 575     Int_t ichannel;
 576
 577     for (Int_t i = 0; i < Nchannel-1; i++) {
 578     ped0[i] = 0;
 579     sig0[i] = 0;
 580     ped1[i] = 0;
 581     sig1[i] = 0;
 582     }
 583
 584     for (Int_t i = 0; i < Nchannel-1; i++) {
 585         inputFile >> ichannel;
 586         inputFile >> ped0[i];
 587         inputFile >> sig0[i];
 588         inputFile >> ped1[i];
 589         inputFile >> sig1[i];
 590         fThreshold[i] = fNsigma*sig1[i];
 591     }
 592
 593     inputFile.close();
 594 }